A Two-Step Biopolymer Nucleation Model Shows a Nonequilibrium Critical Point

Biopolymer self-assembly pathways are central to biological activity, but are complicated by the ability of the monomeric subunits of biopolymers to adopt different conformational states. As a result, biopolymer nucleation often involves a two-step mechanism where the monomers first condense to form a metastable intermediate, and this then converts to a stable polymer by conformational rearrangement of its constituent monomers. While existing mathematical models neglect the dynamics by which intermediates convert to stable polymers, experiments and simulations show that these dynamics frequently occur on comparable timescales to condensation of intermediates and growth of mature polymers, and thus cannot be ignored. Moreover, nucleation intermediates are responsible for cell toxicity in pathologies such as Alzheimer's, Parkinson's, and prion diseases. Due to the relationship between conformation and biological function, the slow conversion dynamics of these species will strongly affect their toxicity. In this study, we present a modified Oosawa model which explicitly accounts for simultaneous assembly and conversion. To describe the conversion dynamics, we propose an experimentally motivated initiation-propagation (IP) mechanism in which the stable phase arises locally within the intermediate, and then spreads through additional conversion events induced by nearest-neighbor interactions, analogous to one-dimensional Glauber dynamics. Our mathematical analysis shows that the competing timescales of assembly and conversion result in a nonequilibrium critical point, separating a regime where intermediates are kinetically unstable from one where conformationally mixed intermediates can accumulate. Our work provides the first general model of two-step biopolymer nucleation, which can be used to quantitatively predict the concentration and composition of biologically crucial intermediates.


I. INTRODUCTION
Biopolymer formation is essential to life. Uncontrolled self-association of biological molecules to form non-covalent polymers is implicated in diseases such as Alzheimer's [1], Parkinson's [2], and sickle-cell anemia [3]. In many cases, the monomeric subunits of biopolymers adopt distinct conformational states. As a result, nucleation often proceeds via a two-step mechanism, in which the monomers first condense into a metastable intermediate, and the stable polymer develops by conformational conversion of the monomers within this phase [4][5][6][7][8][9][10]. Intermediates have diverse morphologies ( Fig. 1(a)), depending on the monomer interactions which drive condensation. Polymeric intermediates are commonly observed in experimental and computational studies of amyloid formation [4,[11][12][13], and convert to stable polymers via reorganization of their secondary or tertiary structure. In addition, diverse biopolymers including actin and prions are known to nucleate in spheroidal intermediates [5-7, 10, 14-16]. Biopolymer nucleation intermediates are implicated in cell death and resulting In order to accurately predict the rate and mechanism of two-step biopolymer nucleation, theoretical models must realistically represent the process by which intermediates convert to stable polymers. Existing analytical models treat conversion as an instantaneous process [19][20][21][22][23], neglecting the dynamics by which the stable polymeric phase spreads through the intermediate. However, there is no general reason why the ordering transition should be sufficiently rapid, or cooperative, for this to be valid. In contrast, experimental and computational studies often point to progressive ordering on slow timescales [7,15,[24][25][26][27], indicating that the stable phase arises locally and spreads through the intermediate. This means that condensation and ordering can occur concurrently [15,[24][25][26]. In addition, simulations support an autocatalytic mechanism by which conformationally ordered monomers promote conversion of their neighbors [24,28,29]. This raises the possibility that condensation and ordering may be kinetically coupled, with larger intermediates converting at a higher rate [20,22]; however, the kinetic prerequisites for such a mechanism have not been rigorously determined. Thus, in order to describe the time-evolution of intermediate populations and elucidate the role of autocatalysis in two-step nucleation, there is an urgent requirement for a model which explicitly considers the conversion of intermediates by an autocatalytic initiation-propagation (IP) mechanism.
In this paper, we address these issues. We develop a minimal model of two-step nucleation ( Fig. 1(b)-(f)) in which polymeric intermediates undergo an autocatalytic conversion process analogous to one-dimensional (1D) Glauber dynamics [30]. Through a combination of analytical calculations and stochastic simulations [31], we show that this simple addition to Oosawa's classical model of nucleated polymerization [32,33] produces rich, nonlinear dynamics. In autocatalytic systems, a nonequilibrium critical point arises, separating a regime where intermediates do not accumulate from one where the two-step nature of nucleation becomes apparent. In the vicinity of this critical point, nucleation dynamics resemble a second-order phase transition similar to the two-dimensional (2D) Ising model [34,35]. It is clear from our derivation that this phase transition does not require a particular dimensionality, and we expect our results to be applicable to morphologically diverse intermediates. Thus, we propose that models of explicit conformational conversion have the capacity to reconcile apparently contradictory interpretations of biopolymer nucleation, and resolve important unanswered questions regarding the origin and nature of metastable intermediates.

II. OUTLINE OF THE MODEL
In this section, we outline the mathematical basis of our two-step biopolymer nucleation model. This model combines the established formalism for nucleated biopolymer self-assembly [32,33] with a proposed initiation-propagation (IP) mechanism of conformational conversion, in which the stable phase arises locally within the intermediate and spreads autocatalytically by nearest-neighbor interactions (Fig. 1). In this respect, our IP mechanism resembles 1D Glauber dynamics [30]. Although the proposed dynamics affect all polymerized monomers, we are particularly interested in their effects on macroscopic, observable nucleated polymerization. Therefore, we concentrate on the aspects of these dynamics which are most relevant to the biopolymer nucleation and growth, and leave a more comprehensive treatment of the internal composition of developing biopolymers as a subject for future work.

A. Microscopic processes
We consider nucleated polymerization of a free monomer population maintained at a fixed concentration m 0 , far from equilibrium. This situation arises physiologically, when there is a monomer source, and is also valid for experimental systems on timescales shorter than the typical timescale τ inf for significant monomer depletion to occur, which we discuss in greater detail in Sec. II. B. Each biopolymer is treated as a 1D lattice of monomers Schematic of the proposed two-step biopolymer nucleation model. (a) Biopolymer nucleation pathways frequently involve diverse metastable intermediates, which progress to stable polymers by conversion of monomer subunits. (b-f) Schematic representation of our model, accounting for explicit growth and conversion of linear intermediates as described in Sec. II. A. Color scheme for monomers: green, unassembled; red, metastable (type-1); blue, stable (type-2).
with length j. A polymerized monomer is in either the metastable (type-1) or stable (type-2) state ( Fig. 1(b)-(f)). Polymerization is initiated at a rate v n = k n m j0 0 , and each polymer initially contains j 0 ≥ 2 monomers in the type-1 state ( Fig. 1(b)). Polymerized monomers irreversibly convert to the type-2 state at a rate of k c , or k c + k p if adjacent to a type-2 monomer ( Fig. 1(e) & (f)). Thus, k p represents the excess rate due to autocatalysis. For the purpose of this study, it is not necessary to consider cases where a type-1 monomer is bounded by two type-2 monomers, as such cases do not affect the macroscopic kinetics of nucleated polymerization.
Each polymer can elongate by reversible monomer addition at a single 'active' end ( Fig. 1(c) & (d)). Strong bias towards elongation at a single end is often observed experimentally [36][37][38], and is caused by the fact that most biopolymers exhibit structural differences between their ends, which result in contrasting physicochemical properties. Although inclusion of this bias simplifies our analysis, it is not expected to qualitatively alter the results. Therefore, we fully expect our theory to generalize to more unusual cases [39] where both ends elongate at similar rates. It is also important to note that the forward and reverse elongation rates must exhibit an equal bias, in order to satisfy detailed balance when m 0 is at the solubility limit. From this point onwards, we simply refer to the forward process as elongation, and the reverse process as dissociation.
During elongation, the conformational state of the monomer being incorporated into the polymer is templated by that of the monomer already at the active end [5,40] (Fig. 1(c) & (d)). Therefore, if the monomer at the active end is in the type-1 state, then elongation proceeds via addition of type-1 monomers with a forward rate k + 1 m 0 and a reverse rate k − 1 . Similarly, if the monomer at the active end is in the type-2 state, then type-2 monomers are added with forward and reverse rates k + 2 m 0 and k − 2 , respectively. Thus, polymers exhibit one of two distinct growth modes, type-1 and type-2, corresponding to the conformation of the monomer at the active end and the resulting mechanism of elongation. Growth mode switching occurs when the monomer at the active end converts, causing the emergence of an elongating type-2 phase. Experimental results typically support k + 2 k + 1 [4,41], meaning that the switching rate determines the lifetime of intermediates with a high type-1 structure content, as well as the observed macroscopic kinetics of nucleated polymerization.

B. Characteristic timescales
Before deriving master equations to predict the macroscopic kinetics of nucleated polymerization, it is informative to consider the key conformational transitions that a nascent polymer undergoes, and the timescales on which they occur (Fig. 2). There are three key events that lead to emergence of a type-2 growth mode: (i) initiation of the polymer at a time t 0 , with all constituent monomers initially in the type-1 state; (ii) non-autocatalytic conversion of the first polymerized monomer at time t c > t 0 , causing the formerly 'pure' intermediate to become a mixed intermediate; and (iii) growth mode switching at time t s ≥ t c , which coincides with conversion of the monomer at the active end ( Fig. 2(a)). If the first monomer to convert is situated at the active end, then t s = t c ; otherwise, additional conversion events are required for growth mode switching to occur, so that t s > t c . If k p k c , these additional conversion events will be primarily non-autocatalytic and uncorrelated with one another; if k p k c , conversion of the first monomer will stimulate autocatalytic conversion of its neighbours, causing the type-2 state to propagate from its point of origin to the active end. We also identify three important time intervals: (i) the lifetime of a pure intermediate, τ c = t c − t 0 ; (ii) the propagation time τ p = t s − t c , which is the time for a mixed intermediate to adopt a type-2 growth mode; and (iii) the total switching time τ s = t s − t 0 = τ c + τ p , which is the total time taken for an intermediate to acquire a type-2 growth mode ( Fig. 2(a)). In addition, we introduce a single domain wall propagation time τ p , which is distinct from τ p and will be discussed in greater detail in Sec. II. C.
We also consider two timescales related to the overall self-assembly of populations of mature polymers: (i) the characteristic self-assembly time τ char , and (ii) the inflection time τ inf (Fig. 2(b)). Both are common experimen- Characteristic timescales of two-step biopolymer nucleation. (a) Schematic showing an example of how the length (j; black line) and AEM domain size (x; red line) of a single polymer can change over time. The three key events (Sec. II. B.) of initiation (t = t0), conversion of the first monomer (t = tc), and growth mode switching (t = ts) are annotated with diagrams representing the change in conformational state of the polymer, with the following color scheme for constituent monomers: green, unassembled; red, metastable (type-1); blue, stable (type-2). (b) Schematic showing the macroscopic self-assembly kinetics of a typical polymer population (Sec. II. B.), and relevant characteristic times. Curves represent the proportions of monomers which are unassembled (m(t); green), or incorporated into polymers with a type-1 (M1(t); red) or type-2 (M2(t); blue) growth mode. tal observables which are derived from phenomenological analysis of biopolymer self-assembly curves, and can be used to derive mechanistic conclusions [42,43]. Specifically, the time-dependent mass of mature polymers follows a convex profile under circumstances where m 0 is approximately constant. If monomer depletion can occur, the self-assembly curve is sigmoidal and approaches a constant limit due to exhaustion of the free monomer [42] (Fig. 2(b)). For the purpose of this study, we identify mature polymers as species with a type-2 growth mode, and denote the effective concentration of monomer incorporated into such species as M 2 (t). The self-assembly time τ char is the time required for M 2 (t) to reach an arbitrary, often small proportion of the initial free monomer concentration m 0 , such that M 2 (τ char ) = Cm 0 and C is a constant in the interval (0, 1) [42] (Fig. 2(b)). While a number of different values of C are used, for small C these different characteristic times typically produce similar scaling behaviors. The inflection time τ inf is the time at which M 2 (t) has its inflection point, such that ∂ 2 t M 2 (τ inf ) = 0 [43] (Fig. 2(b)). Because inflection is a consequence of monomer depletion, τ inf provides a general timescale on which monomer depletion occurs. Therefore, the assumption that m 0 is approximately constant implies that τ s τ inf , meaning that M 2 (t) is locally convex.
C. Growth mode switching as an absorbing boundary problem The ability of type-2 monomers to autocatalytically induce conversion of their neighbors allows the type-2 state to propagate to the active end from a distant site of origin. Thus, the rate of growth mode switching depends not only on the conformational state of the monomer at the active end, but also the length distribution of the domain of type-1 monomers adjacent to the active end, which we term the 'active end metastable' (AEM) domain. When the size x of this domain becomes zero, growth mode switching occurs ( Fig. 2(a)). We are primarily interested in far-from-equilibrium cases where m 0 k − 2 /k + 2 , so that growth mode switching has a negligible reverse process and can be treated as an absorbing boundary problem. Initially, the polymer is formed with all monomers in the type-1 state, so that x = j for t 0 ≤ t < t c . Following the first conversion event at t = t c , the polymer acquires a mixed composition, so that x < j for t ≥ t c . For t > t c , the AEM domain can grow due to elongation at the active end ( Fig. 1(c)), or contract due to either dissociation of monomers from the active end ( Fig. 1(d)) or conversion of monomers within the AEM domain ( Fig. 1(e) & (f)). Growth mode switching occurs when the AEM domain size hits an absorbing boundary at x = 0, at the time t s .
The fact that both autocatalytic and non-autocatalytic conversion events may occur during propagation leads to a distinction between the total propagation time τ p , and the single domain wall propagation time τ p (Fig. 2(a)). While autocatalytic conversion is restricted to the nearest neighbors of type-2 monomers, non-autocatalytic conversion is not. Thus, non-autocatalytic conversion of type-1 monomers within the AEM domain may cause a new domain wall to form, superseding the existing AEM domain wall. While τ p denotes the total time between the onset of propagation and growth mode switching, τ p denotes the time taken for a single domain wall to successfully propagate to the active end, without being superseded. Therefore, τ p only applies to propagation of the last AEM domain wall to form before growth mode switching. If the original domain wall successfully propagates to the active end, then τ p = τ p ; otherwise, τ p < τ p .

D. General master equations
We consider the multivariate probability mass function (PMF) p j,x (t; j 0 , t 0 ), which gives the probability that a polymer will have total length j, AEM domain size x, and a type-1 growth mode at time t. The master equation for p j,x (t) has the following form: where δ j,j0 and δ(t − t 0 ) are the Kronecker and Dirac deltas, respectively, and θ(z) is a step function such that θ(z) = 1 if z ≥ 0 and θ(z) = 0 if z < 0. As with previous models [43][44][45], j 0 is taken to be a minimum polymer length below which rapid disassembly almost always occurs; therefore, an absorbing boundary exists such that p j0−1,x (t) = 0. In addition, to account for irreversible growth mode switching, we impose the absorbing boundary condition p j,0 (t) = 0. While the first term in Eq. (1) accounts for initiation at t = t 0 , the terms proportional to k c , k p , k + 1 m 0 , and k − 1 account for non-autocatalytic conversion, autocatalytic conversion, elongation, and dissociation, respectively.
The terms accounting for non-autocatalytic conversion introduce a jump process to Eq. (1), which complicates efforts to satisfy the boundary conditions. Nonetheless, analytical solutions can be obtained in specific cases. In the case where x = j, corresponding to the length distribution of intermediates with a purely type-1 composition, the summation disappears and the problem becomes significantly more tractable. We use this special case to obtain τ c in Sec. III. B. For mixed intermediates, the time-dependence of AEM domain size can be solved by summing over the polymer length. Let φ x (t; x c , t c ) represent the probability that a mixed intermediate has an AEM domain of size x and a type-1 growth mode at time t, given the boundary condition that x = x c when t = t c . We can determine φ x (t) by summing p j,x (t) over j: where j 0 ∨ x + 1 denotes the maximum of j 0 and x + 1.
The corresponding master equation is as follows: As with related theories [42][43][44], we consider the case where k − 1 is small, allowing losses due to destabilization of polymers with length j < j 0 to be neglected. The absorbing boundary at x = 0 corresponds to growth mode switching, and is retained. Eq. (3) is the discrete analog of a jump-diffusion equation: where µ = k + 1 m 0 − k − 1 − k p and D = (k + 1 m 0 + k − 1 + k p )/2 are the drift and diffusion coefficients for AEM domain size, with D ≥ 0 and −2D ≤ µ ≤ 2D. Because Eq. (3) can be written in an equivalent form to Eq. (4), µ and D are also informative in a discrete context. We define an effective Péclet number x c |µ|/D for propagation of the type-2 state to the active end. Propagation is primarily diffusive when x c |µ|/D < 1, and advective when x c |µ|/D > 1. In the advective case, the AEM domain tends to contract if µ < 0, and expand if µ > 0. Thus, µ is the expected rate of AEM domain expansion due to the competing effects of type-1 growth and type-2 propagation, and this competition provides the basis of the critical behavior we observe in this paper.

III. CRITICAL PHENOMENA IN GROWTH MODE SWITCHING
In this section, we investigate how the competing effects of propagation and elongation on AEM domain size alter the lifetime of intermediates. Using a combination of stochastic simulations and analytical theory, we uncover the existence of a nonequilibrium critical point resulting from transcritical bifurcation of an order parameter P e around the point µ = 0. This critical point shares many common features with equilibrium secondorder phase transitions, including a jump in ∂ µ P e , critical slowing, and diverging susceptibility to autocatalytic effects. These effects will be shown in the following section to significantly affect the macroscopic kinetics of two-step biopolymer self-assembly.

A. Existence of a nonequilibrium critical point
To evaluate how the competing effects of propagation and elongation on AEM domain size alter the lifetime of intermediates, we calculate the mean switching time τ s . For this purpose, we use a rejection-free kinetic Monte Carlo (rfKMC; see Supplemental Material at [31]) algorithm to determine hitting times at the boundary x = 0 in Eq. (1). The rfKMC algorithm converges to the exact solution for the dynamics as the number of polymers N → ∞ [46][47][48], and its application to biochemical problems including biopolymer self-assembly is well-established [47,49,50]. In this case, the rfKMC algorithm is preferable to deterministic methods as it samples broad polymer size distributions more efficiently. For the purpose of the simulations, we set k − 1 = 0 so that there is a single absorbing boundary, although our subsequent mathematical analysis (Sec. III. C. & D.; [31]) is more generally applicable to small, non-zero k − 1 . In addition, while we have carried out simulations for diverse j 0 , the results are qualitatively similar in all cases [31]; therefore, we focus on the j 0 = 2 case as a representative example. To explore the system's limiting behavior in a 2D space, it is convenient to normalize both the timescales and rate constants with respect to k c . Therefore, we use relative time k c t and sample points in (k p /k c , k + 1 m 0 /k c ) space by varying a pair of bounded, mechanistically significant parameters, µ/D and ζ = k p /(k c + k p ). While µ/D defines the tendency for expansion or contraction of the AEM domain (Sec. II. D.), ζ is the proportion of the conversion rate adjacent to a domain wall which is attributable to autocatalytic effects. Thus, ζ provides a measure of the importance of autocatalysis in conformational conversion of the assembled monomers.
In Fig. 3(a), we show how the normalized mean growth mode switching time k c τ s varies as a function of the model parameters (µ/D, ζ). We use these data to plot a nonequilibrium phase diagram, which we divide into three different dynamical regimes ( [31]; Fig. 3(c)): (A) random, (B) cooperative, and (C) competitive. In the random switching regime, elongation and autocatalysis do not occur, so switching depends entirely on nonautocatalytic conversion of the monomer at the active end. In the cooperative regime (ζ → 1 when µ/D 0) (Fig. 3), autocatalytic conversion events allow the type-2 state to spread from a distant site of origin to the active end. Thus, non-autocatalytic conversion events occurring anywhere in the intermediate can indirectly lead to growth mode switching, causing the instantaneous switching rate to approach k c j. This means that mass accumulation and autocatalytic conversion cooperate to bring about a change in growth mode. In the competitive regime, which arises for µ/D 0, the effects of elongation on AEM domain size outweigh those of autocatalytic conversion. Therefore, the AEM domain undergoes net expansion, which inhibits propagation of the type-2 state to the active end. This causes an increase in k c τ s , allowing accumulation of large, metastable intermediates with a mixed composition. In this regime, mass accumulation and autocatalytic conversion can be said to compete with one another. Although regimes A-C are continuous with one another, we observe a singular point at (µ/D, ζ) = (0, 1), where the second derivative of τ s with respect to µ/D diverges in a manner resembling a second-order phase transition ( Fig. 3

B. Lifetime of pure intermediates
The switching time τ s must reflect a similar critical point in the lifetime of pure intermediates τ c , or the propagation time of mixed intermediates τ p , since τ s = τ c + τ p . The fact that the critical point occurs on the µ = 0 line, along which the competing effects of autocatalytic conversion and elongation on AEM domain size are balanced, indicates that an effect on propagation is responsible for the observed criticality. To ascertain whether τ c also plays a role in the critical point, we solve Eq. (1) in the case x = j to obtain a moment-generating function (MGF) for τ c . As our simulations were carried out in the k − 1 = 0 case, we apply the same condition when determining τ c . The MGF of the first-conversion time is given by By taking the Laplace transform of Eq. (1) in the case x = j and rearranging, we obtain the following result [31]: where γ is a lower incomplete gamma function. Thus, Eq. (6) decreases monotonically with k + 1 m 0 /k c , and has the limits τ c → (k c j 0 ) −1 as k + 1 m 0 /k c → 0, and τ c → 0 as k + 1 m 0 /k c → ∞. This behavior reflects the fact that elongation increases the number of monomers within the polymer which can non-autocatalytically convert, reducing the time until the first monomer converts. Thus, if growth does not occur (k + 1 m 0 /k c → 0), then the total conversion rate is simply k c j 0 . Conversely, if growth is rapid, then the length goes to infinity and the time until the first conversion event becomes negligible. Importantly, Eq. (6) does not predict any form of criticality, and the fact that it does not depend on k c means that it cannot be expressed as a function of µ/D and ζ alone. This prediction is confirmed by analysis of k c τ c values from our rfKMC trajectories, which we present in Fig. 4(a). Given that k + , the phase behavior of k c τ c in Fig. 4(a) is accounted for by Eq. (6) alone, with no critical point. Thus, changes in τ p rather than τ c must be responsible for the observed criticality. The critical point in τ s must be caused by an abrupt change in the manner in which the type-2 state propagates to the active end. This deduction is supported by analysis of k c τ p values from our rfKMC trajectories, which we present in Fig. 4(b). While k c τ c is unaffected by the critical point, k c τ p exhibits second-order critical behavior similar to k c τ s . Specifically, mixed intermediates are kinetically unstable (k c τ p = 0) along a line from (µ/D, ζ) = (−2, 1) to (0, 1), corresponding to the idealized case in the autocatalytic limit where propagation of the AEM domain wall to the active end is instantaneous. Because the value of k c τ p is constant along this line, it follows that ∂ τ p /∂(µ/D) = 0 along this line as well. At the critical point at (µ/D, ζ) = (0, 1), we observe a jump in ∂ τ p /∂(µ/D) = 0, allowing k c τ p to take nonzero values to the right of this point. To understand this result, let us consider the propagation of a single AEM domain wall from its site of origin to the active end. For this purpose, we solve Eq. (3) in the case where k c x k p , so that the terms proportional to k c become negligible. This describes the propagation of a single domain wall when the probability of additional non-autocatalytic conversion events is negligible (ζ → 1), and x remains finite. We obtain φ x (t; x c , t c ) using lattice Fourier transforms and the image method (see Supplemental Material at [31]): for x ≥ 0 and t ≥ t c , where I x (t) is a modified Bessel function of the first kind and ν = 4D 2 − µ 2 . The switching rate, which is equal to the hitting rate at the boundary x = 0, is given by , which is the probability that a given AEM domain wall does not reach the active end as t → ∞. In the ζ → 1 limit, this can only occur if x goes to infinity. In this case, additional nonautocatalytic conversion events become non-negligible, allowing propagation to be interrupted by formation of a new domain wall closer to the active end. Thus, P e (x c ) is the probability that the AEM domain does not disappear without additional non-autocatalytic conversion events occurring. From final value theorem, we find that where r ± = ± s 2 + 4Ds + µ 2 . As shown in Fig. 5(b), Eq. (8) delineates a parabola-like curve, which opens to the right and whose upper (r − , unphysical) and lower (r + , physical) branches meet at the point (s * ,ĥ * ), where s * = ν − 2D ≤ 0. When µ = 0, both branches intercept theĥ axis; when µ = 0, (s * ,ĥ * ) touches theĥ axis and the branches exchange intercepts, causing an abrupt change in ∂ µ P e (x c ): In the ζ → 1 limit, autocatalytic conversion and growth occur on much faster timescales than nonautocatalytic conversion. Therefore, the mean switching time is asymptotically k c τ s ∼ lim j→∞ [31]. Thus, Eq. (9) explains the critical point in k c τ s , and the agreement between our calculated scaling law and the simulated result is shown in Fig. 5(a). The emergence of the critical point in the ζ → 1 limit (Fig. 5(a)) is remarkably similar to the behavior of the 2D Ising model in the case of vanishing magnetic field H [35]. Therefore, we propose that the appearance of nonzero P e (x c ) for µ > 0 represents a second-order phase transition, where P e (x c ) is an order parameter and (µ, ζ) play an analogous role to (T, H) in the Ising model. This transition is associated with the emergence of a regime where mixed intermediates have a finite lifetime, and their continued formation allows a nonequilibrium steady-state population to persist.
D. Diverging propagation times and susceptibility at the critical point As further evidence for this critical point, we observe that the moments of the propagation time diverge when (µ/D, ζ) = (0, 1) (Fig. 5(c)). As discussed in Sec. II. C., τ p is the time taken for a single domain wall to successfully propagate to the active end, without being interrupted by formation of a new domain wall closer to the end ( Fig. 2(a)). In the ζ → 1 limit, the n th -order moment is given by (τ p ) n = (−1) n lim s→0 + (∂ n sĥ )/ĥ. The upper and lower branches ofĥ meet at s * = 0 when µ = 0, so the derivatives ∂ n sĥ diverge and the moments behave correspondingly. More generally, which strongly resembles the diverging magnetic susceptibility χ H = ∂ H M near the Curie point. Because τ p is the time interval in which non-autocatalytic conversion events can interrupt propagation, a connection to a susceptibility of the form χ ζ ∝ ∂ ζ P e is highly intuitive.
To test this connection, we determine analytic bounds for the susceptibility in the ζ → 1 limit. We define the susceptibility as χ ζ = −∂ ζ P e = k p ∂ kc P e , where the sign of our definition ensures a positive value, but does not affect the presence or absence of a divergence. This is essentially the effect of ζ on the probability that propagation will be interrupted by a non-autocatalytic conversion event in the AEM domain. By differentiating Eq.
(3) with respect to k c , we obtain a master equation which can be used to determine ∂ kc φ 1 (t, k c ; x c , t c ) [31]. In turn, this can be used to assess the impact of k c on the hitting rate and P e . We solve this master equation to obtain upper and lower bounds for ∂ kc φ 1 (t, k c ; x c , t c ): Thus, there is a divergence in χ ζ at the critical point which has a scaling law between |µ| −1 and |µ| −3 and is closely related to the divergence of τ p and Var(τ p ).

E. Relationship to a transcritical bifurcation
The exchange of intercepts by the branches ofĥ(s; x c ) closely resembles a transcritical bifurcation (Fig. 5(b)), raising the question of whether P e (x c ) can be written as a fixed point of the corresponding normal form. We begin by noting that the escape probability can be expressed has not yet propagated to the active end, under conditions where k c is sufficiently small that additional nonautocatalytic conversion events are negligible. Because P e (x c ) = 1 − [1 − P e (1)] xc , a bifurcation in P e (1) must cause a similar behavior in P e (x c ) [31]; therefore, we focus on the case where x c = 1. Using the conservation Either by applying final value theorem and the power rule of limits [51] to Eq. (12), or by taking the inverse Laplace transform and evaluating the convolution that corresponds to theΦ(s; 1) 2 term, it is possible to obtain the following result in the t → ∞ limit [31]: This is the normal form of a transcritical bifurcation. The physical and non-physical branches of Eq. (8) correspond to the regions of Eq. (13) where ∂ t Φ(t; 1, t c ) decreases or increases with Φ(t; 1, t c ), respectively, and the exchange of intercepts between these branches corresponds to an exchange of stability between the fixed points of Φ(t; 1, t c ). Solutions to Eq. (13) provide insights into the diverging mean and variance of τ p at the critical point, as the hitting rate can be derived from Φ(t; 1, t c ) according to the formula h(t; 1, [31]. When µ = 0, the solution is a logistic function, so that the hitting rate is exponentially bounded. However, when µ = 0, the solution is ∼ (k + 1 m 0 t) −1 , meaning that the hitting time distribution becomes heavy-tailed and its moments diverge. This is an example of critical slowing, and is caused by the disappearance of the linear term in ∂ t Φ(t; 1, t c ). In equilibrium systems, a relationship between critical slowing and diverging susceptibility can be established by considering the Landau potential V (Φ), which satisfies the law ∂ t Φ = −∂ Φ V (Φ). Both the relaxation time following a small variation, and the susceptibility of Φ to a linearly coupled external field, are inversely proportional to ∂ 2 Φ V (Φ). At the critical point, disappearance of the linear term in ∂ t Φ means that ∂ 2 Φ V (Φ) = 0, so that both the relaxation time and susceptibility diverge. In our system, the connection between critical slowing and the susceptibility is more direct, as the mean propagation time defines a lower bound for the susceptibility. Nonetheless, the principle is similar, in that the exchange of stability by the fixed points causes disappearance of the linear term in Eq. (13), allowing the susceptibility to another parameter to diverge.
It is possible to obtain an effective potential of the form V (Φ) = −µΦ 2 /2 + k + 1 m 0 Φ 3 /3 by integrating Eq. (13) with respect to Φ. However, this result must be interpreted with caution. Thermodynamic equilibrium relies on microscopic reversibility, which does not exist for our parameter Φ(t; 1, t c ). Due to the irreversible nature of growth mode switching, Φ(t; 1, t c ) can only decrease; thus, the regions of Eq. (13) for which ∂ t Φ(t; 1, t c ) > 0 are non-physical, and correspond to the parts ofĥ(s; 1) situated at s < 0 ( Fig. 5(b)). Instead, we propose that the shared feature between our critical point and those observed in equilibrium systems is the presence of a bifurcation. While our bifurcation is the result of irreversible dynamics, bifurcations observed in equilibrium systems are caused by large numbers of microscopically reversible processes. Nonetheless, because these processes result in similar dynamical laws, our model has many features in common with equilibrium phase transitions.

IV. EFFECTS OF CRITICALITY ON THE MACROSCOPIC SELF-ASSEMBLY KINETICS
Autocatalytic biopolymer systems with high ζ will exhibit a nonlinear response to µ near the critical point. In this section, we evaluate the effects of this criticality on the macroscopic self-assembly kinetics. We begin by estimating likely ζ values based on typical biological free energy scales. We find that stabilization of conformational transition states by even modest free energy differences (< 10k B T ) can result in highly autocatalytic conformational conversion. Thus, different biopolymers are expected to exhibit the full spectrum of autocatalytic behaviors predicted by our model. In systems with high ζ, the critical point has a profound effect on the macroscopic kinetics, inducing rapid changes in the size and abundance of steady-state populations of intermediates, and the concentration-dependent scaling behavior of the formation of mature polymers. These effects are not predicted by existing local equilibrium models [52][53][54], and indicate that nonequilibrium criticality may explain widespread experimental phenomena [55,56] which onestep nucleation models have not been able to address.

A. Estimation of ζ from Kramers theory
Critical behavior arises in the autocatalytic (ζ → 1) limit ( Fig. 3(b)). Although k c is non-zero and k p is finite, meaning that a value of ζ = 1 is not strictly possible, ζ can be asymptotically close to this limit. Realistic ζ values can be estimated from reaction rate theory. Let us consider the conversion rate of a single type-1 monomer in the AEM domain, situated n monomers from the active end. We represent the conformation of the monomer with a single continuous degree of freedom, q n ; however, we note that our analysis can be generalized to cases where monomers have additional slow degrees of freedom without changing the conclusions. The free energy G(q n ; q n−1 , q n+1 ) of a polymerized monomer depends on its own conformation q n , and a term accounting for interactions with the adjacent monomers q n−1 and q n+1 . The monomer at position n has free energy minima corresponding to the type-1 and type-2 states, which are situated at q n = p 1 and q n = p 2 . These minima are separated by a barrier at q n = p ‡ , which corresponds to a conformational transition state. The fact that conversion has a negligible reverse pro-cess entails that G(p 1 ; q n−1 , q n+1 ) G(p 2 ; q n−1 , q n+1 ), and the effectively two-state nature of our model entails that G (p 1 ; q n−1 , q n+1 ) and G (p 2 ; q n−1 , q n+1 ) are both large. Therefore, the conformation q n of all nonconverting monomers n = n is assumed to be close to p 1 or p 2 . Monomers situated closer to the active end (n < n) are part of the AEM domain, and must be in the type-1 state (q n ≈ p 1 ); monomers situated further from the active end (n > n) may be in either the type-1 or type-2 state (q n ≈ p 1 or q n ≈ p 2 ). Although it is possible to write a Langer equation [57,58] for the conversion rate which accounts for continuous conformational variation in the monomers at positions n = n, the results are hard to interpret. The main reason for this is that such a model predicts dynamics based on the j-dimensional free energy function n G(q n ; q n−1 , q j+1 ), which has 2 j−1 relevant saddle points. Instead, we simplify the analysis by neglecting variation in q n about the points p 1 or p 2 . We write the free energy of the monomer at position n as a pair of functions corresponding to the two possible conformations of the monomer at position (n + 1): The monomer conformation exhibits ovderdamped Langevin dynamics, so that the corresponding conformational probability density functions ρ 1 (q n , t) and ρ 2 (q n , t) satisfy the following Fokker-Planck equation: where i = 1, 2 is the conformation of the monomer at position (n + 1), β = 1/k B T , and D i (q n ) is the conformational diffusion coefficient. The conversion rate can be estimated from G i (q n ) and D i (q n ) using Kramers theory [58,59]. The non-autocatalytic conversion rate is given by the passage rate at the barrier (q n , q n+1 ) = (p ‡ , p 1 ): where . Similarly, the autocatalytic conversion rate is given by the passage rate at the saddle point (q n , q n+1 ) = (p ‡ , p 2 ): where ∆G ‡ 2 = G 2 (p ‡ ) − G 2 (p 1 ). Therefore, the definition ζ = k p /(k c + k p ) gives the result: where ∆∆G ‡ = ∆G ‡ 1 − ∆G ‡ 2 is the reduction in the free energy barrier of conversion due to autocatalytic effects. It is not immediately obvious whether the preexponential ratio would have a value greater or less than 1. While ordering of the (n + 1) monomer might be expected to reduce conformational diffusion, the simultaneous loss of complementarity between the monomers' structures might exert the opposite effect. The effects on G (p 1 ) and G (p ‡ ) are also hard to gauge, and may be biopolymer-specific. In the absence of further information, we set the value of the pre-exponential ratio to 1, so that ζ = 1 − exp (−β∆∆G ‡ ), and we examine the relationship between ζ and ∆∆G ‡ . Eq. (18) predicts that even high ζ values will be associated with relatively modest free energy differences; for example, ζ = 0.999 ⇔ k p ≈ 10 3 k c implies ∆∆G ‡ = 6.9 k B T , which is well within the range of typical biological free energy scales at T = 310 K. For comparison, formation of a single hydrogen bond in protein has a free energy change of ∼ 2 k B T [60,61], and estimated free energy barriers for protein folding are mostly in the 2 − 20 k B T range [62], accounting for the effects of large numbers of competing, non-additive interactions. In addition, a recent experimental study of biopolymer nucleation by the Aβ peptide indicated that autocatalytic interactions between polymers reduce the nucleation free energy barrier by ∆∆G ‡ = 19 k B T [63]. Although this form of autocatalysis is different from the form that we consider here, as it depends on interactions between polymers rather than within a polymer, the two processes may share a common physical basis. It is interesting to note that an equivalent ∆∆G ‡ value in our model would give ζ ≈ 1 − 10 −8 ⇔ k p ≈ 10 8 k c , resulting in highly autocatalytic behavior. Thus, we believe that values of ζ close to the ζ = 1 limit are highly plausible, meaning that biopolymer systems will frequently exhibit high levels of autocatalysis and resulting critical behavior.

B. Nonlinear response to perturbations in the autocatalytic limit
We now investigate the response of the macroscopic self-assembly kinetics to variations in µ/D, when ζ is high. To do so, we calculate the time-evolution of the principal moments of the polymer length distribution. Let f i,j (t) represent the concentration of polymers with growth mode i = 1, 2 and length j at time t. In parallel with the notation used by [43,44], we define the following principal moments: The zeroth-order moment P i (t) is the total concentration of polymers with growth mode i, and acts as an effective normalization constant for f i,j (t). The first-order moment M i (t) is the effective concentration of monomers incorporated into polymers with growth mode i, or the 'polymer mass-concentration'. The mean polymer length L i (t) can be obtained by the normalization: Let us begin by examining the effect of the critical point on M 1 (t) and L 1 (t). The length distribution f 1,j (t) of the population of type-1 species contains contributions from individual polymers initiated at all t 0 < t. If nucleated polymerization begins when t = 0, then where v n = k n m j0 0 is the nucleation rate. Thus, In Sec. III. A., we obtained numerical solutions for p j,x (t; t 0 ) using an rfKMC algorithm. We will now use these solutions to predict M 1 (t) and L 1 (t), by applying Eq. (23-24). As we are specifically interested in the response of the system to varying µ/D at high ζ, we will focus on cases where µ/D varies through the range [−2, 2), while ζ is fixed at the value ζ = 0.999 ⇔ k p ≈ 10 3 k c . This value is associated with a comparatively small ∆∆G ‡ , and so is physically very plausible (Sec. IV. A.). In Fig. 6, we present the effect of µ/D on the timeevolution of αM 1 (t) and L 1 (t), where α = k c /v n is a normalization constant facilitating the use of relative time k c t. The integral ∞ 0 p j,x (t; t 0 ) converges, allowing both αM 1 (t) and L 1 (t) to approach steady-state values αM ss = α lim t→∞ M 1 (t) (Fig. 6(a) & (b)) and L ss = lim t→∞ L 1 (t) (Fig. 6(b)) in the t → ∞ limit. These nonequilibrium steady states reflect the concentrations and size of intermediates likely to be observed in physiological contexts, if the free monomer concentration and chemical environment are stable on timescales τ steady , where τ steady is the approximate timescale taken to reach the steady state. In experimental contexts, the pre-steady-state behavior is likely to dominate in the self-assembly lag phase if τ inf τ steady , and the steadystate behavior is likely to dominate if τ inf τ steady . We will evaluate the effects of the critical point on both the pre-steady-state kinetics and the steady-state kinetics. The pre-steady-state kinetics appear linear when plotted on double-logarithmic axes ( Fig. 6(a)), indicating that they are described by a power law of the form αM 1 (t) ∝ t d . Power laws are often observed in biopolymer self-assembly kinetics where t is small, and d can be interpreted as the number of sequential processes contributing to biopolymer mass accumulation; thus, nucleation without polymerization gives d = 1, and nucleated polymerization gives d = 2. To assess whether the presteady-state kinetics are affected by the critical point, we calculated d = ∂ ln M 1 (t)/∂ ln t| kct=0.01 for varying µ/D (Fig. 6(b)). Although d increases with µ/D, reflecting the increasing importance of polymerization in mass accumulation, there is no clear response to the critical point. This is unsurprising, given that the critical point results from changes in growth mode switching, which is negligible in the pre-steady-state.
In contrast, the steady-state kinetics are strongly affected by the critical point. We begin by examining the effect of the critical point on τ steady . For the purpose of this analysis, we define τ steady such that αM 1 (τ steady ) = 0.9αM ss . As shown in Fig. 6(a), τ steady initially decreases with increasing µ/D, as both elongation and propagation cooperate to accelerate growth mode switching. However, this trend reverses close to the critical point, causing a delay in the approach to the steady-state and allowing greater accumulation of polymer mass. This delay corresponds to the emergence of a population of large, conformationally mixed intermediates. As a result, both αM ss and L ss exhibit a highly nonlinear response to µ/D around the critical point ( Fig. 6(a) & (b)). In physiological and experimental contexts, this means that small variations in the free monomer concentration or solvent composition have the potential to induce large changes in the populations of biopolymer nucleation intermediates. In addition, the total rate of growth mode switching is strongly dependent on the size and nature of the intermediate populations, so the increase in τ steady will strongly affect the accumulation of mature polymers with a type-2 growth mode.

C. Low concentration-dependence is a dynamical signature of a nonequilibrium critical point
As discussed in Sec. II. B., a common experimental descriptor of the macroscopic self-assembly kinetics is the characteristic time τ char at which M 2 (τ char ) = Cm 0 , where C is an arbitrary constant ( Fig. 2(b)). The concentration-dependence of τ char can be described using a scaling exponent γ, which has the following definition: In cases where C is sufficiently small (0 < C 0.5), mechanistic conclusions can be drawn from the value of γ [42]. For example, the Oosawa model predicts γ = n c /2 [33], where n c is the nucleation order, and the nucleationconversion-polymerization (NCP) model, which assumes instantaneous growth mode switching, predicts γ = n c /3 [21]. To illustrate the impact of the abrupt change in τ steady on the macroscopic self-assembly kinetics, we predict γ values based on our simulations. We begin by observing that: Our rfKMC simulations were carried out in the k − 1 → 0 limit, meaning that the following conservation principle applies: To simplify our analysis, let us consider the case where k + 2 m 0 M ss /τ s , so that the predominant contribution to the type-2 polymer mass comes from elongation, rather than conversion of intermediates. Therefore, In the pre-steady-state (t τ steady ), P 1 (t) is approximately linear, so both P 2 (t) and M 2 (t) will exhibit higher-order scaling with time. In the steady-state (t τ steady ), P 1 (t) is approximately constant, meaning that P 2 (t) and M 2 (t) will exhibit linear and quadratic scaling with time, respectively. It is more convenient to express concentrations in relative units than it is to choose arbitrary values of the rate constants. Therefore, we calculate α M 2 (t) using Eq. (23,(27)(28), where α = (k + 1 /k p ) j0+1 (k 2 c /k n k + 2 ) is a normalization constant, and set C = 0.01/α . From the definition of µ and D, we find that Thus, under circumstances where the rate constants are kept fixed, variation in µ/D can be attributed to changes in the free monomer concentration m 0 . This means that we are able to determine γ by calculating τ char for varying, finely spaced µ/D values, estimating the logarithmic derivative with respect to µ/D from the corresponding finite difference equation, and then multiplying by the right-hand side of Eq. (29). In Fig. 6(c), we present the results obtained for various ζ. At low ζ, γ monotonically transitions between the limits lim µ→−2 γ = 1 and lim µ→2 γ = 2/3, which agree with the predictions of [33] and [21] discussed above. However, a minimum in γ appears at high ζ, which starts at the point µ/D = 2 and moves towards µ/D = 0 as ζ → 1. This depression occurs because µ/D scales with m 0 , so that increasing m 0 causes a reduction in growth mode switching, inhibiting further acceleration of type-2 mass accumulation. In highly autocatalytic systems (k p > 10 3 k c ), the effect is severe enough to cause a temporary reversal in concentrationdependence. As the minimum γ decreases, the drop in γ becomes increasingly sharp about the point µ/D = 0. Both the low concentration-dependence (γ < 0) and sharp drop in γ predicted by our model are commonly observed in experimental studies, but cannot be explained by existing models [52][53][54]. Furthermore, while existing models invoke local equilibria to explain low γ [52][53][54], our model predicts a low γ purely as a result of the dynamics. We thus propose that a low concentrationdependence in biopolymer self-assembly may represent the dynamical signature of nonequilibrium critical point.

V. DISCUSSION
Biopolymer self-assembly is a fundamental mechanism of biological organization, and plays a key role in disease. Experimental and computational studies have revealed that monomers often populate metastable conformational states during the process of self-assembling to form a biopolymer [40,[64][65][66]. As a result, nucleation of a biopolymer often proceeds via one or more metastable intermediates, whose constituent monomers are conformationally distinct from those within the mature polymer. In addition, monomers within the same intermediate often exhibit considerable conformational heterogeneity. This complexity has challenged our ability to extract detailed mechanistic information at the molecular level. Current mathematical models either neglect the conformational heterogeneity of the assembled monomers entirely, or assume that monomers within the same intermediate rapidly assume the same conformational state, so that intermediates with a mixed composition are not observed. However, mixed intermediates are commonly observed in experiments and simulations, indicating that mathematical models which neglect these dynamics are only applicable under highly simplified experimental conditions [64]. Moreover, this conformational diversity is responsible for many of the most important biological behaviors of biopolymers [67][68][69]. One notable example is the role of intermediates in the assembly of amyloid fibres from protein monomers. The conformational heterogeneity of amyloid self-assembly intermediates allows them to interact with a diverse range of other biomolecules [70]. This disrupts cellular function and is causative in the majority of degenerative conditions of aging, including Alzheimer's and Parkinson's diseases [16,18,[71][72][73][74].
In this paper, we address this issue by introducing a minimally complex process by which the self-assembled monomers are able to exhibit conformational heterogeneity, and change conformation over time. This simple addition causes our model to exhibit rich, critical dynamics. As a result, the macroscopic self-assembly kinetics display features such as a highly nonlinear dependence of steady-state intermediate populations on the self-assembly conditions, and loss or even inversion of the concentration-dependence of the self-assembly rate of the stable phase. These phenomena are commonly observed in experimental studies [55,56], but are not explained by existing models. Thus, for the first time, the initiationpropagation (IP) model proposed here allows quantitative predictions to be made regarding the concentration of conformationally heterogeneous intermediates.

A. Competing dynamics in morphologically diverse intermediates
While polymeric intermediates are commonly observed in biopolymer self-assembly, a diverse range of other intermediate morphologies are also observed, including small prenucleation clusters with an apparently random organization of constituent monomers [75,76], large spheroidal intermediates [52,77,78], and species possessing fractal geometry which is evident at a macroscopic scale [79]. In addition to applying to polymeric intermediates, we expect the conclusions of our study to generalize to these other diverse geometries. To make this generalization, we note that propagation involves incor-poration of monomers into a growing stable phase by autocatalytic conversion of condensed monomers in the intermediate phase. Thus, the dimensionality of propagation is defined by that of the stable phase. In an ndimensional intermediate, propagation will manifest as linear growth of the stable phase through the intermediate, with the growth mode switching when an autocatalytic end reaches the solvent. This form of propagation is well-supported by experimental observations of spheroidal intermediates [6,14,15], which are often large enough to be observed by light microscopy techniques. In these studies, the stable polymeric phase appears to originate within the intermediate by conformational rearrangement of a subset of the monomers, and then emerge from the surface of the intermediate. When the active end reaches the solvent, the dramatic change in the physicochemical environment, and the availability and conformational state of nearby monomers, would be expected to induce a marked change in the growth behavior of the polymeric phase analogous to the switch from propagation (k c + k p ) to type-2 growth (k + 2 m 0 ) behavior in our polymer model.
The fact that polymerization is likely to manifest as linear growth-like behavior in multidimensional intermediates means that we expect IP models involving such intermediates to reduce to pseudo-1D absorbing boundary problems. A transition from a regime where P e = 0 to one where P e > 0 is a general feature of such problems [80], so we expect similar transitions for diverse intermediate geometries. It is also important to note that our critical point is conserved in the continuous limit of x, ie. the case where the intermediate is so large as to behave like a continuous phase, rather than a discrete collection of monomers. In this limit, we obtain the classical results for first-passage of Brownian motion with drift [81], and a jump in ∂ µ P e (x c ) occurs about the point µ = 0 [31]. Besides the discrete polymer example considered in this paper, IP mechanisms with a number of other intermediate geometries are likely to be mathematically tractable. For example, in the case of spherical intermediates with a continuous-like composition, the problem considered in Sec. II. can be recast as a continuous advection-diffusion (Eq. (4)) along a 1D trajectory towards a spherical absorbing boundary. In the simplified case of radial propagation, we have performed a quick calculation to estimate the rate at which large spherical intermediates give rise to a growing stable polymers [31]. While this rate is proportional to the volume of the intermediates in the cooperative regime, it is proportional to their surface area in the competitive regime. This behavior closely resembles the way in which initiation sites for successful propagation events are restricted to the ends of polymeric intermediates in the competitive switching regime. Thus, a more general principle may exist, in which the competition between growth and propagation dynamics creates two distinct regimes. In one, a polymerizing stable phase may originate from anywhere within the volume of the intermediate; in the other, the sites where such a phase may originate are restricted to the boundary of the intermediate. The fact that a critical point separating these regimes appears to be a widespread feature of two-step biopolymer nucleation pathways underlines the more general principle that the prerequisites for critical behavior are relaxed in far-from-equilibrium systems.
In addition to being compatible with a wide range of intermediate morphologies, our IP mechanism also predicts significant morphological diversity. The incorporation of a single alternative monomer conformation means that polymers of a given length may exhibit a large number of different sequences of type-1 and type-2 monomers. Because the conformation of a monomer affects its mechanical properties, we would expect these different sequences to result in diverse morphologies. In this study, we have primarily focused on the mechanism of growth mode switching, and have ignored aspects of the conformational conversion dynamics which are not relevant to this problem. However, in future studies it may be possible to rigorously predict changing polymer morphologies by obtaining analytical solutions to describe the progressive conversion of monomers throughout the polymer. As an example of how polymer composition might affect morphology, we observe that monomers in the type-1 and type-2 states are likely to exhibit differing degrees of conformational ordering. It is most often assumed that the intermediate is less ordered than the stable species, so that conformational conversion can be regarded as an ordering transition. In this case, the conformational degeneracy of the type-1 state would be expected to result in greater flexibility for type-1 domains. As the monomers within a polymer converted to the type-2 state, we would expect to see a progressive increase in the persistence length l p . However, in cases where the persistence length of the type-1 and type-2 monomers are markedly different, presence of a small number of type-1 monomers would be expected to strongly affect the overall morphology of the polymer. Thus, polymers with a small but significant type-1 content might exhibit radically different morphologies from mature polymers, as is often observed experimentally [66,82,83].
It would also be possible to introduce further conformational diversity to our model, by introducing additional metastable conformations for the polymerized monomers. While analysis of such a model would not be trivial, additional conformational freedom and the presence of other slow dynamics might be expected to yield even richer behaviors than those presented here. However, it is useful to note that structural data from solidstate nuclear magnetic resonance (NMR) studies suggest that monomers within amyloid fibrils occupy a limited range of stable or metastable conformations [84,85], while high-resolution microscopy data have revealed a diverse range of morphologies [82,86]. Therefore, even in amyloid systems where the free monomer is often highly conformationally degenerate, the range of accessible conformations seems to remain limited within the polymer. Thus, a simple model in which polymerized monomers occupy one of two conformations may be sufficient to account for much of the observed morphological diversity.

B. Explicit conversion dynamics are essential to a general two-step biopolymer nucleation model
The inclusion of explicit conformational conversion dynamics also allows our IP model to exhibit a generality not seen in other two-step nucleation models. Propagation of the stable phase is not assumed to be instantaneous, and can occur on timescales similar to either polymerization of the intermediate or non-autocatalytic conversion (Fig. 2). As a result, the dynamical behavior of our model depends on two key parameters, µ/D and ζ. The first of these, µ/D, describes the competing effects of propagation and polymerization on the time-dependence of growth mode switching. When µ/D < 0, the stable phase propagates through the intermediate faster than the intermediate phase can polymerize, causing accelerated switching. When µ/D > 0, polymerization is rapid and reduces the likelihood of successful propagation. The second parameter, ζ, compares the timescale for propagation due to autocatalytic conversion with the timescale for non-autocatalytic conversion. When ζ = 0, nearestneighbor interactions between monomers do not promote conformational conversion, so that the rate of growth mode switching is independent of the length of the polymer. When ζ = 1, nearest-neighbor interactions strongly promote conversion, so that the timescale for propagation is much faster than that of non-autocatalytic conversion, and the rate of growth mode switching is accelerated.
Existing models which neglect the conversion dynamics are not able to consider the variation accounted for by µ/D and ζ, and thus make various assumptions regarding the rate at which the new growth mode emerges. While certain models [20,22] have correctly supposed a connection between intermediate size and the rate of emergence of a growing stable phase (eg. k c τ s ≈ 1/j), the lack of competing dynamics means they are restricted to a limited range of points close to the top left (ζ → 1, µ/D < 0) of our phase diagram (Fig. 3(a)). In contrast, other models have assumed that there is no relationship between intermediate size and the switching rate [19,21] (ie. k c τ s ∝ 1). Justifications for this behavior include irrelevance of size due to the lack of a growth process, or a scenario in which monomer conversion is highly cooperative. While the former is plausible in specific cases where growth of intermediates is slow, the latter is unlikely to be widely observed, as intermediates are typically either low-dimensional or have highly degenerate monomer interactions. However, a sufficient level of cooperativity may occur in exceptional cases where monomers within intermediates exhibit a limited variety of interactions with a large number of neighbors. In our model, lengthindependence arises in two scenarios: when conversion is entirely non-autocatalytic so that adjacent monomers convert independently (ζ = 0); and when growth is suf-ficiently rapid that growth mode switching is inhibited (µ/D → 2). Thus, the behavior assumed by [19,21] will be observed at points on our phase diagram corresponding to these limits ( Fig. 3(a)). It is important to emphasize that we expect length-independence to be a somewhat unusual phenomenon. Even stabilization of the conformational transition state by a small free energy (∼ 1 k B T ) would result in a significant ζ value, and in the vast majority of our phase space some level of reduction in k c τ s is observed (Fig. 3(a)).
We have seen that existing models of two-step biopolymer nucleation describe a limited range of possible scenarios, which are restricted to points close to three lines (ζ → 1 given µ/D < 0; ζ = 0; and µ/D → 2) on the limits of our phase diagram (Fig. 3(a)). As a result, these models exhibit markedly different behaviors which appear difficult to reconcile. By explicitly considering the dynamics by which conversion of monomers within a polymeric intermediate leads to emergence of a new growth mode, we are able to explore a much broader region of parameter space, and show that these apparently contrasting scenarios are in fact different limits of the same underlying reaction scheme. In addition, we find that these previously unexplored regions contain rich phase behavior which is not anticipated based on the behavior observed in the limits. Therefore, our model has the capacity to explain a wide variety of 'poorly behaved' experimental or physiological systems which are not described by existing models, such as those with mixed intermediates or low concentration-dependence (Sec. IV. B. & C.). Moreover, by classifying existing two-step nucleation models in the context of our phase diagram ( Fig. 3(a)), we have shown that these models occupy disparate regions of phase space which can only be united by considering the competitive dynamics accounted for by µ/D and ζ. Therefore, while our model still does not describe some of the more complex scenarios that have been suggested [23], we have shown that explicit conversion dynamics are essential features of a general model of two-step biopolymer nucleation. We thus believe that our study represents a crucial step towards such a model.

C. Biological importance of the nonequilibrium critical point
In addition to generalizing existing models, the inclusion of an explicit conformational conversion process allows our model to predict a nonequilibrium critical point. While critical points have previously been described in biopolymer systems, the possibility of nonequilibrium criticality has largely been overlooked. Instead, the most commonly suggested critical point is the critical micellar/oligomer concentration (CMC) [52,77,87], which is the concentration at which the free monomer and intermediate are in local equilibrium, such that m 0 = k − 1 /k + 1 . CMC models typically assume that formation of intermediates causes rapid depletion of the free monomer, so that k − 1 /k + 1 represents an upper limit on the free monomer concentration; as a result, saturation behavior occurs. In this paper, we consider a different scenario. We focus primarily on the case where the free monomer depletes slowly, as occurs in physiological circumstances where there is a monomer source, or in many experimental situations [88][89][90][91]. Thus, while our basic model (Fig. 1) has the required processes to produce a CMCtype effect, we are primarily interested in the case where the free monomer is supersaturated with respect to the intermediate (m 0 > k − 1 /k − 1 ), and thus out of local equilibrium. Under these circumstances, a second critical point emerges at m 0 = (k − 1 + k p )/k + 1 (Fig. 3 & 5), at a higher concentration than the CMC. On our nonequilibrium phase diagram (Fig. 3a), the CMC would manifest as a vertical line (ie. constant µ/D) situated to the left of the nonequilibrium critical point. The difference in concentration between these two points will be given by k p /k + 1 . Although we set k − 1 = 0 in our simulations to focus primarily on the nonequilibrium critical point, our mathematical analysis is valid for more general k − 1 subject to the condition that destabilization of intermediates is negligible. Thus, the CMC and the nonequilibrium critical point described in this paper are compatible so long as the concentration difference between the two is sufficiently large. This situation is highly relevant to nonequilibrium scenarios where the dissociation of monomers from the intermediates is relatively slow, as occurs in Alzheimer's and Parkinson's diseases [89][90][91]. In cases where k p k − 1 , the CMC and nonequilibrium critical point will become close to one another on the phase diagram, and the presence of two non-negligible absorbing boundaries in Eq. (1) will further complicate the dynamics of self-assembly and conversion. A rigorous analysis of this effect is both mathematically involved and outside the scope of this paper, and so has been left as a subject for future work.
It is interesting to note that in cases where physiological and experimental systems exhibit slow variation in free monomer levels over time, the phase behavior of biopolymer self-assembly will vary accordingly. Experimental systems are likely to exhibit a monotonic depletion of the free monomer. In this case, µ/D will decrease over time, causing the system to pass the nonequilibrium critical point en route to the CMC. The formation of stable polymers due to growth mode switching will then cause further monomer depletion, resulting in either destabilization or switching of the remaining intermediates. Thus, accelerated growth mode switching will occur when the free monomer concentration passes the nonequilibrium critical point at m 0 = (k − 1 + k p )/k + 1 , with destabilization of the remaining intermediates occurring shortly after. This sort of behavior has been observed for the amyloid-β (Aβ) peptide. Under conditions of high supersaturation, Chimon et al. [7] described the formation of large, spheroidal intermediates by Aβ, similar to those discussed in Sec. V. A. Over time, these intermediates progressively acquired secondary structure content similar to the stable polymeric phase, but growing, stable polymers only appeared following depletion of the free monomer. This behavior is strongly indicative of a transition from the competitive switching regime to the cooperative regime, as described by our model.
The biological implications of the nonequilibrium critical point extend beyond the biopolymer self-assembly process itself. The rapid change in steady-state behavior when µ/D ≈ 0 (Fig. 6) means that small changes in free monomer concentration or the physicochemical environment will have a pronounced effect on the population of intermediates. This provides a way for biological systems to rapidly switch between producing stable polymers and metastable intermediates with distinct biological activity. This sort of switch-like regulation is typical of biological systems [68,92], and may for example reflect the importance of controlling the population of potentially toxic species in the assembly of amyloids [93]. Thus, while previous work is restricted to specific sets of conditions, our model encompasses a much wider context. This makes it possible to predict the response of the system to a wide range of stimuli, including changes in internal conditions during the course of a reaction or external stimuli such as mutations and temperature changes, and to interpret the possible role that particular variations in self-assembly conditions play in healthy functioning and disease.
Biological systems are typically maintained far-fromequilibrium, with free monomer concentrations often highly supersaturated in relation to their respective polymers [88]. However, existing models of two-step biopolymer nucleation neglect the dynamics which arise under these nonequilibrium conditions. The requirements for criticality are relaxed under nonequilibrium conditions, and many of these neglected dynamics have the potential to produce complex, biologically relevant behaviors, such as we observe in our IP model. As discussed above, these behaviors may explain experimental observations which are not predicted by existing models, and rationalize the role of physiological conditions in human disease. As a result, we propose that the role of nonequilibrium critical phenomena in physiological biopolymer self-assembly has been severely underestimated.

D. Relationship to autocatalytic secondary nucleation
The model developed here can be applied to any twostep biopolymer nucleation process. A fundamental aspect of our model is its ability to predict mixed intermediates consisting of self-assembled monomers in two distinct conformations, resulting in diverse intermediate morphologies. Conversion of a self-assembled monomer from the metastable type-1 state to the stable type-2 state occurs either non-autocatalytically, with rate k c , or autocatalytically with rate k c +k p . Thus, the stable phase is able to initiate locally within the intermediate, and then propagate autocatalytically to the active end. This type of intramolecular structural propagation has been extensively documented in the biochemical literature in models of multi-subunit proteins [7,15,24,26,27,94], and is responsible for the nonequilibrium criticality that our model exhibits.
A different type of autocatalysis which is often discussed in models of biopolymer nucleation is secondary nucleation, the process by which the surface of an existing stable biopolymer catalyzes the formation of additional biopolymers by heterogeneous nucleation [95][96][97][98]. Thus, while the proposed IP mechanism can be considered a form of intra-polymer autocatalysis, secondary nucleation is an example of inter-polymer autocatalysis. Secondary nucleation has been used to explain features such as exponential (M 2 (t) ∝ exp(κt)) rather than quadratic (M 2 (t) ∝ t 2 ) accumulation of the stable phase, and the formation of significant concentrations of toxic intermediates during the 'growth phase' (ie. when t ≈ τ inf ) [97]. Given the prominence of secondary nucleation models in the protein polymerization literature, it is pertinent to question how they are mechanistically related to our IP model. At the microscopic level, the process of secondary nucleation results from the alignment of free monomers to the surface of the mature biopolymer [99]. In the context of our model, this alignment could promote initiation of type-1 polymers which then detach from the mature polymer and independently undergo growth mode switching. This would manifest as an increase in the initiation rate v n . There is also the possibility that interactions between conformationally ordered monomers in the stable phase and disordered monomers in the nascent intermediate promote ordering of the latter. This would create a more tightly coupled process in which conversion and propagation are accelerated, resulting in an increase in k c and k p . Thus, it is possible that the IP mechanism and secondary nucleation share a common mechanistic basis, and are simply different examples of a common tendency for autocatalytic conformational change induced by interactions between monomers. All these effects are easy to incorporate into the framework of our IP model, and we anticipate that future studies which combine the IP mechanism with inter-polymer autocatalysis may be able to gain deeper insights into the microscopic mechanism of secondary nucleation.

E. Competing dynamics can explain low concentration-dependences
One of the key experimental validators for biopolymer assembly models is their ability to accurately describe the concentration-dependence of the process and relate this to theory [5,23,32,33,40,55,56,97,100,101]. The Oosawa model [32,33] was first to ascribe a physical meaning to the value of γ (Eq. (25)), the scaling exponent for the concentration-dependence of accumulation of the stable polymer phase (Sec. IV. C.). The Oosawa model predicts that γ = n c /2, where n c is the effective order of nucleation and often represents the minimum size for a polymer to be stable or metastable (ie. j 0 ) [42]. While n c is often interpreted as reflecting the number of monomers within the critical nucleus (n * = n c − 1) [102], which is the least stable species during the self-assembly pathway, this explanation has been challenged by the widespread observation of γ < 1 [55]. Such low values of γ would imply n * < 1, which is physically implausible [56]. In addition, many biopolymer systems with higher γ still exhibit much lower concentration-dependences than would be expected, based on biophysical predictions of their n * [20,23,103]. Some authors have proposed that saturation effects, such as the presence of a saturable catalytic surface which induces heterogeneous nucleation, may explain the low concentration-dependences [53,54,104]. Indeed, γ < 1 is generally observed under conditions where γ appears to saturate with increasing monomer concentration. However, under some experimental conditions the concentration-dependence drops further to γ < 0.5 [55]. This is not only incompatible with the Oosawa model, but it also occurs under circumstances where secondary nucleation remains active, meaning it is incompatible with saturation models. To our knowledge, the nucleation-conversion-polymerization (NCP) [21] model is the only model in the literature which produces low γ without invoking a saturation effect. In the NCP model, low γ values are explained in the context of a 'cascade nucleation system', where sequential conversion of N −1 pre-fibrillar intermediates results in a concentration-dependence of n c /(N + 1) ≤ γ ≤ n c /2. In this model, although the number of intermediate types can be infinite, their conversion occurs as a single-step process without considering the underlying dynamics, and the intermediates do not exhibit a growth process. As a result, the conversion mechanism of an NCP intermediate corresponds to the scenario predicted at the point (µ/D, ζ) = (−2, 0) on our nonequilibrium phase diagram ( Fig. 3(a)).
The IP mechanism proposed here generalizes the conversion mechanism considered in the NCP model to explicitly consider the dynamics by which intermediates grow and individual monomers convert. As a result, our model makes predictions regarding the pre-steady and steady-state scaling behaviors of intermediates and mature polymers which can easily be generalized to the case of multiple intermediates in the same manner as the NCP model. However, the inclusion of additional dynamics allows our model to explore a broader region of parameter space, leading to our observation of sharp changes in γ near the critical point ( Fig. 5(b) & (c)). This results in local depression of the concentration-dependence, which is more pronounced and occurs more abruptly at higher ζ.
High ζ values correspond to systems in which autocatalysis plays a dominant role in conformational conversion, and can be achieved on comparatively modest free energy scales (Sec. IV. A.). This reduction in concentrationdependence allows our far-from-equilibrium model to violate constraints on the concentration-dependence of ex-isting models over a broad range of µ/D values, and this effect is likely to be further enhanced when saturation effects are also active. Under many experimental conditions, a value of γ < 1 is thus expected, and is simply a result of the competing dynamics. In this study, the competing processes are polymerization of the metastable intermediate, and propagation of the stable conformational state through this intermediate. However, models which explore other nonequilibrium dynamics are also likely to produce such an effect. Therefore, our work highlights the potential that far-from-equilibrium selfassembly models with competing dynamics have to resolve outstanding issues relating to low concentrationdependences.
F. Interpretation of the nucleus size Fig. 6(c) shows how γ varies with µ/D for different ζ values. While both µ/D and ζ are likely to depend on the physicochemical environment under which selfassembly occurs, and the underlying properties of the biopolymer system, µ/D exhibits an additional dependence on the free monomer concentration m 0 , which can easily be varied in experimental contexts. It is useful to note that different ζ values result in highly characteristic changes in concentration-dependence (Fig 6(c)), so that the ζ value intrinsic to a particular biopolymer system can be determined by experimentallly varying µ/D and measuring the changing concentration-dependence of assembly of the stable polymer. Thus, while it is possible to determine n c from the values of γ observed in the limits where µ/D → ±2, our model predictions provide an easier alternative. Instead, we anticipate that it will be possible to determine n c and ζ from the characteristic curve shape with which γ varies at intermediate values of µ/D. This may be preferable as it eliminates the need to explore extreme values of µ/D, which is challenging due to lack of sensitivity of experimental methods at low concentrations, and the likely presence of additional competing processes at high concentrations. This approach can also be generalized to cases similar to the NCP model [21] in which multiple intermediates are present. In this scenario, we would typically expect the Oosawa behavior (γ = n c /2) to be recovered in the µ/D → −2 limit, and the NCP-type behavior (γ = n c /(N + 1)) to be recovered in the µ/D → 2 limit, as we have observed in Sec. IV. C. The inclusion of extra intermediates would also be expected to introduce additional critical points, which could be experimentally observable.
It is also important to note that, while our current model identifies n c as being equal to the minimum size of a metastable intermediate j 0 , future models which consider additional complexities of the conversion dynamics may arrive at different interpretations for n c . For example, a previous coarse-grained simulation study [22] which investigated a scenario similar to our cooperative switching regime identified n c with the minimum size at which conversion can initiate within an intermediate, rather than the absolute minimum size of an intermediate. While we take these sizes to be the same in our model, future work which treats these quantities as different may provide analytical justification for this result.

VI. CONCLUSIONS
We have developed a minimal model of two-step biopolymer nucleation which accounts for conformational conversion of metastable intermediates by an initiationpropagation (IP) mechanism. In our framework, the stable polymer phase initiates locally within the intermediate, and propagates by autocatalytic nearest-neighbor interactions between the assembled monomers. This mechanism is well-supported by experimental studies showing progressive ordering of intermediates on slow timescales [7,15,[24][25][26][27], and computational studies suggesting an autocatalytic mechanism for structural conversion of the monomers [20,24,28,29]. The inclusion of conformational ordering on timescales comparable to condensation of the intermediates causes our model to exhibit rich dynamics and a nonequilibrium critical point, which are not predicted by models with an implicit ordering step [19][20][21][22][23]. Furthermore, by explicitly considering the dynamics of conformational conversion, our framework unifies existing models, and shows that their apparently unconnected behaviors occur as different limits of a single underlying mechanism. The IP mechanism has the potential to explain a wide variety of experimental behaviors which are not predicted by existing models, such as the formation of large, conformationally mixed intermediates [6,7,14,15], and nucleation of the stable polymer phase at a rate independent of the free monomer concentration [55]. In our model, these phenomena arise naturally as a result of the competing dynamics of selfassembly and conversion, resolving the apparent para-dox that concentration-independent nucleation rates suggested critical nuclei with a non-physical size [56].
Biopolymer nucleation intermediates are responsible for toxicity in a huge range of human diseases [16,18,[71][72][73][74], and the nonequilibrium critical point predicted by our model may explain the extreme sensitivity of many polymerization disorders to small changes in self-assembly conditions [105][106][107]. As a result, our model provides a deeper understanding of the origin of these toxic intermediates, and suggests promising therapeutic strategies to treat diseases caused by their formation. Moreover, future theoretical studies which build on the results presented here are likely to yield analytical solutions for the size distribution and composition of nucleation intermediates, which will provide a powerful tool in therapeutic development. Biopolymers are also increasingly exploited by materials scientists [108][109][110], and the development of a general model of their nucleation will allow us to regulate their formation and prevent unwarranted side reactions. In this paper, we have shown that explicit conformational conversion dynamics, such as those addressed by our initiation-propagation mechanism, are essential to such a theory. We thus believe that our work represents a crucial step towards a general theory of two-step biopolymer nucleation, with important implications for both experiments and human disease.