Protein accumulation in the endoplasmic reticulum as a non-equilibrium phase transition

Several neurological disorders are associated with the aggregation of aberrant proteins, often localized in intracellular organelles such as the endoplasmic reticulum. Here we study protein aggregation kinetics by mean-field reactions and three dimensional Monte carlo simulations of diffusion-limited aggregation of linear polymers in a confined space, representing the endoplasmic reticulum. By tuning the rates of protein production and degradation, we show that the system undergoes a non-equilibrium phase transition from a physiological phase with little or no polymer accumulation to a pathological phase characterized by persistent polymerization. A combination of external factors accumulating during the lifetime of a patient can thus slightly modify the phase transition control parameters, tipping the balance from a long symptomless lag phase to an accelerated pathological development. The model can be successfully used to interpret experimental data on amyloid-\b{eta} clearance from the central nervous system.

A number of neurological pathologies, including Alzheimer's disease (AD) and Parkinson's disease, spongiform encephalopathies and serpinopathies, are collectively identified as conformational diseases because they are all characterized by the aggregation and tissue deposition of aberrant conformations of proteins. In recent years, a large effort has been devoted to understanding the biological mechanisms underlying the biochemical properties of these proteins and the associated secretory pathways. While the key mechanisms triggering these diseases are still largely unknown, recent evidence points towards a pivotal role played by secretory pathways that is common to all these pathologies.
Serpinopathies result from point mutations in a1-antitrypsin and neuroserpin showing a delay in folding, with unstable intermediates being cleared by endoplasmic reticulum (ER)associated degradation (ERAD) [1][2][3] . The remaining proteins are either fully folded and secreted or retained as ordered polymers within the ER of the cell of synthesis. Mutations in neuroserpin result in the autosomal dominant inclusion body dementia (Familial Encephalopathy with Neuroserpin Inclusion Bodies (FENIB)) 4 , and mutants of neuroserpin in FENIB patients show accelerated rates of polymerization compared with wild type protein both at the protein level 5,6 and in cell models 7,8 . The current picture of these diseases shows that accumulation of misfolded proteins within the lumen of the ER activates ERAD pathway. When ERAD fails to remove sufficient mutants, the remaining proteins start to polymerize 7,9,10 .
In spongiform encephalopathies, recent evidence shows a new pathway that contributes to prion propagation owing to proteasomal dysfunction and ER stress, leading to an increase of prion protein in the secretory pathways 11 . The classical hystopathology of AD consists of extracellular plaques of the amyloid-b peptide and intracellular neurofibrillary tangles of hyperphoshorylated aggregates of microtubule-associated protein tau. Ab-degrading proteases involve a number of subcellular compartments, including lysosomes, endosomes and the ER 12 . Recent results show an overproduction of aggregation-prone amyloid b42, which might interest the ER 13 . This is confirmed by the analysis of the brains of AD patients displaying ER dysfunction 14,15 , particularly in the ERAD pathways 16 . Synucleinopathies are a group of neurodegenerative disorders, including Parkinson's disease, that are associated with proteins containing a-synuclein molecules. Experimental evidence suggests that the deposition of a-synuclein in intracellular inclusions together with ubiquitin lead to impaired functions in the ubiquitin-dependent proteasome system [17][18][19][20][21] . Finally, Huntington's disease, an inherited autosomal dominant disease, is caused by the expansion of CAG repeats within the HTT gene encoding huntintin, which leads to a protein prone to aggregation in the cytoplasm and alterations in the secretory pathways 22 .
The traditional theoretical framework to understand protein polymerization involved in these neurological disorders is based on either molecular dynamics simulations, which give an understanding of how individual proteins interact 23,24 , or on kinetic rate equations, which yield the growth of the polymers in a mean-field approximation [25][26][27][28][29][30][31] . Theoretical results are then usually compared with in vitro measurements of protein aggregation by dynamic light scattering or single molecule fluorescence under well defined conditions of temperature and concentration 30,[32][33][34] . The conditions are, however, not physiological because in the cell proteins undergo a well regulated cycle of synthesis by the ribosome and subsequent degradation through secretory pathways.
In this paper, we demonstrate that taking correctly into account the physiological conditions is important, because rates of synthesis and secretion control a non-equilibrium phase transition for the formation of protein aggregates in the ER. This can be seen in two complementary models. The first is a set of Monte Carlo simulations in a three dimensional (3D) confined space representing the ER, with synthesis and degradation explicitly accounted for. The second is a mean-field model of the same scenario, which offers analytical estimates of long-time protein aggregation. The results from these models illustrate that small variations in control parameters, which could be induced by a variety of biological factors, can lead to large variations in outcomes. This opens a new common perspective on possible diagnostics for these pathologies.

Results
Protein polymerization in vitro. We consider a 3D model of linear polymer aggregation in a confined space simulating an intracellular organelle like the ER, a set of channel-like structures surrounding the cell nucleus (see Fig. 1 for an illustration of the model and the Methods section for a detailed description). To validate our model in a well-studied scenario, we first consider the kinetics of a typical in vitro experiment and perform 3D numerical simulations of linear polymer diffusion and aggregation at constant concentration and temperature (with k out ¼ k in ¼ 0 and periodic boundary conditions as discussed above). While the model is general for linear polymers, we focus here on neuroserpin, which has been shown to undergo processes of activation and latentization 33,34 .
We start the simulations from a fixed number N m of inactive monomers, randomly distributed in space, and then study the effect of different concentrations, measured by the dimensionless density r N m L 3 0 =L 3 . We quantify the aggregation kinetics by measuring the weighted polymer mass m w ¼ /i 2 S//iS, where i is the number of monomers in each polymer and the average is taken over different realizations of the process. This observable is accessible experimentally through dynamic light scattering 34 . Figure 2a shows that for k f ¼ 0 the weighted mass kinetics displays a crossover between a short-time regime growing as t 2 , owing to activation and dimerization ( Supplementary Fig. 1), and a long-time regime growing as t b with bC0.5, owing to polymerpolymer aggregation. The results are in agreement with experiments 34 . We fit the curves obtained at different densities to a crossover function (see Supplementary Information) depending on a crossover time t scaling as r À g (see the inset of Fig. 2b). The best fit yields g ¼ 0.149 ± 0.002 and allows rescaling of all the curves onto a single master curve. This result confirms that the concentration only sets the timescale of the kinetics. In Fig. 2c, we study the effect of polymer fragmentation showing that if k f 40 the polymerization process is blocked after a characteristic time that depends on k f . The role of latentization is illustrated in Fig. 2d, which shows that the long-time growth is not affected if k L 40. Latentization induces a plateau in the crossover region, a feature that is also observed in experiments (see the inset of Fig. 2d and ref. 34). Mean-field theory can also be used to study polymerization kinetics in vitro, and yields qualitative agreement with the 3D model, albeit with quantitative differences, as shown in Fig. 3. However, the two models agree on essential features.
Polymerization in vivo displays a non-equilibrium phase transition. To describe in vivo conditions, we introduce to our models non-zero rates for protein synthesis (k in ) and degradation (k out ). These effects are in competition: protein synthesis allows greater polymer growth via a flux of monomers that combine into larger polymers; however this growth can be balanced by polymer degradation. To study this quantitatively, we turn first to the mean-field model. Here we characterize polymer aggregation by the mean length m p of polymers of size iZ2, which is given by the ratio of the mass and the number of polymers, P i!2 in i = P i!2 n i . As discussed above, when protein synthesis and degradation are turned off, as in experiments performed in vitro, the mean-field model predicts a finite steady state mean polymer size 26 . However, if protein synthesis is turned on, with rate k in , then if k out is small, the mean polymer size increases as Ct b , with b ¼ 0.5. On the other hand, if k out is greater than a critical value k c out , which depends on the rates of aggregation and fragmentation, a non-growing steady state is again obtained ( Fig. 4a and Supplementary Fig. 2 and 3). As shown in Fig. 4b, we find a sharp transition between growing and stationary phases as quantified by the prefactor of the square-root growth of the mean polymer size C, which scales to zero as ðk out À k c out Þ y , with y ¼ 1/2, when the transition is approached. Figure 5 shows the phase diagram estimated by the mean-field model as a function of the rates k in , k out , k f and k p . The transition region is governed by the products k in k p and k out k f . Intuitively, k in k p controls the rate of polymer growth, by the introduction of new monomers and their aggregation into larger polymers, while k out k f controls the rate of polymer shortening by the removal of material and fragmentation of large polymers. When k f is large relative to k in k p , the dependence on parameters is further reduced to a dependence on the ratio k in k p /k out k f . Indeed, in the limit c ¼ 1, it can be shown that if a steady state exists, the relationship ( This relation determines the location of the critical region as discussed in detail in the Supplementary Information. We have confirmed the existence of a polymerization phase transition controlled by protein production and degradation by numerical simulations in a 3D system representing the ER, thus mimicking more closely physiological conditions in the cell. As shown in Fig. 4c, 3D simulations show growth curves analogous to the one obtained in the mean-field model, namely m p BCt 1/2 (Fig. 4a). In Fig. 4d, we report the prefactor C estimated from the fits in Fig. 4c and observe transition curve similar to the one observed in the mean-field curve (Fig. 4b), with a typical broadening expected for finite size systems.
Polymerization in vivo displays critical-point fluctuations. To confirm that we are in presence of a critical phase transition, we study the temporal fluctuations of the average polymer length in 3D simulations. As shown in Fig. 6a for a single realization of the process, the polymer length undergoes an intermittent dynamics with bursty activity that is reminiscent of the crackling noise observed in materials close to non-equilibrium phase transitions 35 . We characterize the statistical properties of the signal by measuring the distribution of pulse durations T and sizes s, defined as the area under a pulse, for different values of k out . Figure 6b shows that for k out 4k c out , bursts display a power law distribution with a cutoff that diverges at the transition. In the growing phase, the distribution develops a peak at large T typical of system-wide spanning events. We fit all the subcritical distributions simultaneously by the scaling function obtaining a ¼ 1.94 ± 0.02, d ¼ 2.6 ± 0.2 and k c out ¼ 1:82 AE 0:02 as best fit parameters. Using these values we can collapse the distributions onto a universal scaling function as expected for critical phenomena (see Fig. 6c). A similar analysis for the size distribution yields a power law exponent k ¼ 1.75 ± 0.01, cutoff exponent 1/s ¼ 2.5±0.2 and k c out ¼ 1:89 AE 0:02 (see Supplementary Fig. 4). We notice that these value differ significantly from the predictions of mean-field theory, which Polymerization is system-size dependent. Three dimensional simulations of protein aggregation display qualitative agreement with mean-field theory in conditions that should apply to experiments in vitro and in vivo. Quantitative differences between the two models are mainly owing to the absence of geometrical features in the mean-field model. Investigating the role of spatial confinement is important because large local fluctuations can drive polymerization in small systems. To this end, we use 3D simulations to estimate the size-dependence of some effective mean-field parameters such as the polymerization and degradation rates, as a function of k out and H. The effective degradation rate k out depends on the sample size H in three dimensions because proteins can exit only through the boundaries, a feature that is absent in the mean-field model (see Fig. 7). The effective polymerization rate k p displays a peak as a function of k out , in correspondence with the phase transition (Fig. 8). Furthermore, for smaller system sizes the polymerization rate is larger.
Comparison with experimental data on amyloid-b clearance from the brain. Although quantitative experimental data for protein aggregation in the ER are not yet available, the kinetic mechanism we study is general and can be applied also to other situations such as Ab aggregation and clearance in the brain 12 , which has been subject to extensive experimental study 37,38 . In particular, we use our mean-field model to simulate an amino acid labelling experiment like those reported by Bateman et al. 37 .
In those experiments, labelled leucine, an amino acid composing Ab, was injected into the bloodstream of healthy subjects 37 . After a 9 h infusion, the concentration of labelled Ab in cerebral spinal fluid was monitored for a further 27 h.
To understand experimental results, we simulate the scenario using our mean-field model. We first obtain a steady state by simulating protein aggregation in the steady phase for 50 years. We then follow the labelling experiment of ref. 37 so that over 36 h a fraction of injected monomers are labelled. Labelled monomers evolve following the same dynamics as unlabelled monomers, with full mixing. We calculate the labelled/unlabelled mass fraction in the cerebral spinal fluid, as reported in Fig. 8a. (c) Allowing polymer fragmentation limits the growth of polymers. Here the activation rate is k A ¼ 1, the polymerization rate k p ¼ 1, and there is no latentization (k L ¼ 0), and the growth in unweighted mean size of polymers of length iZ2, m p has been plotted. (d) Latentization slows polymer growth, as seen here for the weighted mean polymer size m w . Activation is not incorporated here, and no fragmentation occurs (k f ¼ 0), and k p ¼ 1. We plot against time rescaled by k p . When fragmentation is allowed, the system always attains a steady state in which all mass consists of latent monomers. Analogous plots for the 3D model are shown in Fig. 2. From this, we determine the fractional clearance rate (FCR) of labelled proteins, as described in the Supplementary Methods (see Supplementary Fig. 5). Measured FCRs depend on many parameters (see Supplementary Table 1) but here we focus on the role of protein degradation, controlled by k out . A single FCR measurement does not allow a determination of all parameters, but our results indicate that a series of labelling experiments can reveal changes in degradation rates over time and the approach to the physiological-pathological phase transition. We take the steadystate system and decrease k out slightly, and allow it to evolve a further 3 years to attain a new steady state, after which we perform another labelling test. As reported in Fig. 9b, a reduced k out yields a reduction in FCR, as well as the maximum labelled/ unlabelled ratio. Further decreases in k out cause measured FCR to continue to decrease, but as the critical point is approached the lifetime of labelled proteins diverges and near the critical point no clearance is observed in experimental timescales. Further experiments by the same group indeed reveal that FCR decreases for AD patients with mild symptoms 38 , in agreement with the predictions of our model.

Discussion
Conformational diseases, such as Alzheimer's and Parkinson's diseases, spongiform encephalopathies and serpinopathies, are all associated to aberrant protein aggregation in which the secretory activity of intracellular organelles plays a critical role. In this paper, we have shown by 3D numerical simulations and meanfield calculations that protein aggregation undergoes a nonequilibrium phase transition controlled by the rates of protein synthesis and degradation. Our theoretical analysis compares well with experimental results both in vitro and in vivo. We can describe with great accuracy the time-dependent polymerization of neuroserpin in vitro 34 and experimental measurements of Ab clearance from the central nervous system in healthy subjects (or with mild AD symptoms) 37,38 . We predict that Ab clearance measured experimentally in AD patients should decrease as the disease progresses and vanish as the non-equilibrium phase transition is approached. At this point, protein polymerization would not stop. We believe that combining quantitative clinical  Phase diagram of the mean-field model. The manifold in parameter space separates those parameters for which a steady state is obtained (the stationary phase) and those that lead to continual growth of polymers (the growing phase). The model used here has a polymer degradation rate that scales with polymer size i as i À 3 , up to a cutoff c ¼ 5, but similar phase diagrams are obtained for different size dependences and cutoffs and in much of the phase space the location of the transition is predicted well by a model with c ¼ 1, as discussed in the Supplementary Information. information with realistic theoretical modelling could be useful to improve diagnostic tools for these pathologies.
In our model, the phase transition is associated with a breakdown of homeostasis in intracellular organelles, controlled by the rates k in and k out . This result naturally leads to two fundamental questions: (i) which specific biological processes tune these control parameters and (ii) what are the biological consequences of the transition?
A possible answer to the first question might involve lipid metabolism, which has been recently shown to play a role in all these neurodegenerative diseases [39][40][41][42] . In particular, ref. 42 shows that in FENIB the inhibition of HMGCoA reductase, the limiting enzyme of the cholesterol biosynthetic pathway, has a critical role in the clearance of mutant neuroserpin from the ER. Numerous studies have also reported that the modification of cholesterol content can affect amyloid precursor protein processing, which is needed for neuronal activity 40 . Interestingly, a recent paper shows that E693D (Osaka) mutation in amyloid precursor protein promotes intracellular accumulation of Ab, reducing its excretion 43 . This paper also shows the importance of Ab trafficking for intracellular cholesterol transport and efflux and that the Osaka mutation potentiates cholesterol-dependent exacerbation of intracellular Ab toxicity by disturbing amyloidmediated cholesterol efflux from the cell 43 . These observations can be translated in our model considering that the alteration of lipid metabolism (that is, the level of cholesterol) should lead to a reduction of the parameter k out eventually triggering the phase transition.
As for the second question, our model allows the interpretation of the emergence of conformational diseases in the framework of phase transitions and critical phenomena. This implies that a minimal change in the control parameters can lead to drastic changes in the system when we cross the transition line leading to the disease. The continuous nature of the transition implies, however, that large fluctuations are expected as we approach it. - In perspective this might provide the basis for early detection strategies before the pathology has taken place.
It is important to remark that the phase transition we observe does not occur when the concentration is held constant as usually assumed in experiments in vitro. In this case, the response to the variation of control parameters is smooth: small changes of the rates correspond to small changes in the polymerization kinetics. This might be one of the reason behind the difficulty in extrapolating in vitro results to the in vivo situation. One possibility to overcome this problem would be to use microfluidic devices where the influx and outflux of proteins are externally controlled 44 . Using a similar experimental device, it would be possible to reproduce more faithfully protein aggregation in the cell.

Methods
3D model. In the model, individual proteins are modelled as monomers sitting on a 3D square lattice. Monomers diffuse with rate k D and attach to neighbouring monomers or polymer endpoints with rate k H . Polymers move collectively by reptation with a length-dependent rate k R /i 2 , where i is the number of monomers in the polymer (see ref. 45 p. 89), and locally by end rotations with rate k E and kink moves with rate k K (for a review of lattice polymer models see ref. 45). A polymer can attach to another polymer with rate k H if their endpoints meet, and can fragment by breaking an internal bond with rate k f (Fig. 1c). Inspired by experimental results on neuroserpin polymerization 34 , we allow for polymerization after at least one of the monomers has been activated with rate k A . Active monomers can also become latent with rate k L and after that they do not aggregate (Fig. 1d).
We consider two types of boundary and initial conditions for the polymerization kinetics. (i) To simulate experiments in vitro, we start with a constant number of inactive monomers in a cubic system of size L ¼ 60L 0 , where L 0 is the typical monomer diameter, with periodic boundary conditions in all directions. (ii) To simulate polymerization in the ER, we consider a system of size (L Â L Â H) with L ¼ 100L 0 and H ¼ 25L 0 with periodic boundary conditions along x and y and closed boundary conditions along z. This choice is justified by the structure of the ER in which a channel of small width (that is, H) is bounded by two extended membrane sheets. Monomers enter the system from both closed boundaries with rate k in and monomers and polymers can exit from the same boundaries with rate k out /i 3 . The fact that the exit rate decreases with the polymer length is suggested by experiments showing that proteasome degradation is slower for larger aggregates 46,47 , although no specific measurements exist for the ER. We have checked that different functional dependencies of the exit rate on the number of monomers in the polymer yield similar results.
We perform numerical simulations using Gillespie Monte carlo algorithm 48 . In case (i) we measure time in units of 1/k D , setting varying r, k L and k f . In case (ii) we measure time in units of 1/k in and set Mean-field model. In the mean-field model, the evolution of the populations n i of polymers of size i is described by a set of coupled nonlinear differential equations, which are based on the assumption that polymer aggregation and fragmentation have no spatial dependence. We have used the model introduced by Blatz and Tobolsky 26 in which polymers aggregate with rate controlled by parameter k p and fragment with rate controlled by k f . We add two ingredients to this model. The first is production of monomers, which enters as a term þ k in added to _ n i , the rate of change of monomer population. The second is polymer degradation, which occurs at a rate that decreases with polymer size. This is described by a set of terms À k out n i i À 3 . The population of polymers of size i, iZ1, is given by the differential equation where f(i) ¼ i À 3 for irc and 0 otherwise. Rate equations for the moments of the size distribution and equations for _ n i , irc are solved numerically using standard techniques. The cutoff c is introduced only for numerical convenience, because otherwise it would be impossible to solve the equations in closed form, but we have observed rapid convergence to a c independent solution when f(i) decreases faster than i À 2 . The limit c ¼ 1 may have biological significance since it implies that only monomers can exit the system, which implies that polymers have to be fragmented before being degraded. A full discussion of these equations is given in the Supplementary Methods.