{'fraction_non_alphanumeric': 0.029706835304680487, 'fraction_numerical': 0.010969278213686224, 'mean_word_length': 4.615519508987286, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fluctuations in protein aggregation: Design of preclinical screening for early diagnosis of neurodegenerative disease

Autocatalytic fibril nucleation has recently been proposed to be a determining factor for the spread of neurodegenerative diseases, but the same process could also be exploited to amplify minute quantities of protein aggregates in a diagnostic context. Recent advances in microfluidic technology allow analysis of protein aggregation in micron-scale samples potentially enabling such diagnostic approaches, but the theoretical foundations for the analysis and interpretation of such data are so far lacking. Here we study computationally the onset of protein aggregation in small volumes and show that the process is ruled by intrinsic fluctuations whose volume dependent distribution we also estimate theoretically. Based on these results, we develop a strategy to quantify in silico the statistical errors associated with the detection of aggregate containing samples. Our work opens a new perspective on the forecasting of protein aggregation in asymptomatic subjects.

Autocatalytic fibril nucleation has recently been proposed to be a determining factor for the spread of neurodegenerative diseases, but the same process could also be exploited to amplify minute quantities of protein aggregates in a diagnostic context. Recent advances in microfluidic technology allow analysis of protein aggregation in micron-scale samples potentially enabling such diagnostic approaches, but the theoretical foundations for the analysis and interpretation of such data are so far lacking. Here we study computationally the onset of protein aggregation in small volumes and show that the process is ruled by intrinsic fluctuations whose volume dependent distribution we also estimate theoretically. Based on these results, we develop a strategy to quantify in silico the statistical errors associated with the detection of aggregate containing samples. Our work opens a new perspective on the forecasting of protein aggregation in asymptomatic subjects.

I. INTRODUCTION
The presence of aberrant conformations of the amyloid β peptide and the protein α-synuclein is considered to be a key factor behind the development of Alzheimer's and Parkinson's diseases, respectively. The polymerization kinetics of these proteins has been shown to consist of nucleation and growth processes and to be strongly accelerated by the presence in solution of pre-existing fibrils [1,2], thereby circumventing the slow primary nucleation of aggregates. It was found that surfaces, such as lipid bilayers [3,4] and hydrophobic nanoparticles [5] can accelerate the nucleation process dramatically. Indeed, in the case of α-synuclein, it was found that in the absence of suitable surfaces, the primary nucleation rate is undetectably slow [2]. Under certain conditions, the surfaces of the aggregates themselves appear to be able to catalyse the formation of new fibrils, leading to autocatalytic behavior and exponential proliferation of the number of aggregates [2,6,7]. This so-called secondary nucleation process is likely to play an important role in the spreading of aggregate pathology in affected brains [8], as the transmission of a single aggregate into a healthy cell with a pool of soluble protein might be sufficient for the complete conversion of the soluble protein into aggregates.
An intriguing idea is to exploit this observation to screen biological samples based on the presence of very low concentrations of aggregates for pre-clinical diagnosis of neurodegenerative diseases. Indeed, this has been achieved in the case of the prion diseases in a methodology called prion misfolding cyclic amplification [9], which is based on the amplification of aggregates through repeated cycles of mechanically induced fragmentation and growth. Recently, the applicability of this approach to the detection of aggregates formed from the amyloid β peptide has been demonstrated [10]. Furthermore, the autocatalytic secondary nucleation of amyloid β fibrils has been exploited to demonstrate the presence of aggregates during the lag phase of aggregation [11].
However, none of these methods currently allow to easily determine the absolute number of aggregates in a given sample. One strategy to address this problem is to divide a given sample into a large number of sub volumes and determine for each of the sub-volumes whether it contains an aggregate or not. Due to advances in microfluidic technology and microdroplet fabrication [12], it is now possible to monitor protein aggregation in micron-scale samples [13], a technique that could be used to design microarrays targeted for protein polymerization assays. To be successful this program needs guidance from theory to quantify possible measuring errors due to false positive and negative detection. Current understanding of protein polymerization is based on mean-field reaction kinetics that have proved successful to describe key features of the aggregation process in macroscopic samples [7,14,15]. This theory is, however, designed to treat the system in the infinite volume limit, where the intrinsic stochasticity of the nucleation processes cannot manifest itself, so that its applicability to small volume samples is questionable. The importance of noise in protein aggregation was clearly illustrated in Ref. [16], who proposed and solved the master equation kinetics of a model for polymer elongation and fragmentation, obtaining good agreement with experiments on insulin aggregation [13].
Here we address the problem by numerical simulations of a three dimensional model of diffusion-limited aggregation of linear polymers [17], including explicitly auto-catalytic secondary nucleation processes [2,6,7]. A three dimensional model overcomes the limitations posed by both mean-field kinetics [7,14,15] and master equation approaches [16], which do not consider diffusion and spatial fluctuations. Most practical realizations of protein aggregation reactions are not diffusion limited, due to the slow nature of the aggregation steps, caused by significant free energy barriers [18]. This leads to the system being well mixed at all times and mean-field theories providing a good description. There are, however, cases both in vitro (e. g. when protein concentrations and ionic strengths are high, leading to gel formation [2]) and in vivo (due to the highly crowded interior of the cell), where a realistic description cannot be achieved without the explicit consideration of diffusion. Here we use our model to study fluctuations in the aggregation process induced by small volumes and to provide predictions for the reliability of a seed detection assay.

THREE DIMENSIONAL MODEL
Simulations are performed using a variant of the protein aggregation model described in Ref. [17] where individual protein molecules sit on a three dimensional cubic lattice. The model considers primary nucleation due to monomermonomer interaction, polymer elongation due to addition of monomers to the polymer endpoints and secondary nucleation processes in which the rate of monomer-monomer interaction is enhanced when the process occurs close to a polymer (see Fig. 1a for an illustration) [19]. In particular, monomers diffuse with rate k D and attach to neighboring monomers with rate k M (primary nucleation) but when a monomer is nearest neighbor to a site containing a polymer composed of at least n * monomers, then the nucleation rate increases to k S k M (secondary nucleation). We do not consider polymer fragmentation, since this term is mostly relevant for samples under strong mechanical action [7], and some of the most important amyloid-forming proteins have been shown to exhibit aggregation kinetics dominated by secondary nucleation under quiescent conditions [2,7]. A monomer can attach to a polymer with rate k H if it meets its endpoints. 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. [20] 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 [20]). Simulations start with a constant number of N monomers in a cubic system of size L = m 0 L 0 (with m 0 an integer) where L 0 is the typical monomer diameter, with periodic boundary conditions in all directions. We perform numerical simulations using the Gillespie Monte Carlo algorithm [21] and measure time in units of 1/k S and rates in units of k S . We explore the behavior of the model by varying independently both the monomer concentration ρ ≡ N/m 3 0 and the number of monomers N at fixed ρ, but also the rate constants. For the simulations results reported in the following, the rates describing polymer motion are chosen to be k E = k R = k K = 10 −2 which is smaller or equal than the diffusion rate of the monomers k D .
As expected, secondary nucleation efficiently decreases the half-time before rapid polymerization. We illustrate this by changing the critical polymer size n * needed to induce secondary nucleation. We observe that the lower n * the shorter the half-time (see Fig. 1(b)). Currently, no experimental data exists on the value of n * , but it can be expected to be of a similar magnitude as the smallest possible amyloid fibril, defined as the smallest structure for which monomer addition becomes independent of the size of the aggregate and an energetic downhill event.

MEAN-FIELD THEORY
The progress of reactions observed experimentally in bulk systems can be well approximated by a mean-field model [7,14,15], without fragmentation or depolymerization of polymers. Such models are in contrast to our three dimensional computational model, which describes also monomer diffusion and polymer motion due to reptation, kink motion and end-rotations, which are not treated by mean-field approximation. Despite this, it is still possible to fit polymerization curves resulting from three dimensional simulation results through mean-field theory with effective diffusion dependent parameters. The fact that both experimental and simulated polymerization curves are described by the same mean-field theory ensures that our model is appropriate to describe experiments. In the mean-field model, the evolution of the concentration f j of polymers of length j ≥ n c , where n c is the nucleation size, is given by [15] where dots indicate time derivatives and m(t) is the monomer concentration. The first term on the right-hand side represents an increase in the concentration of polymers of size n c due to polymer nucleation by aggregation of n c monomers with rate constant k n ; this is a generalized version of the dimer formation with rate constant k M in the 3d model. The second term represents an increase in the concentration of polymers of size j by attachment of a monomer to a polymer of size j − 1, with rate constant k + . The third term is the corresponding loss of concentration of polymers of size j when they attach to a monomer. These two terms are the mean-field equivalent of the endpoint attachment of monomers to polymers with rate constant k H in the 3d model. The final term represents secondary nucleation, which in the mean-field model is described as an increase in concentration of polymers of size n 2 (the secondary nucleus size) occurring at a rate proportional to the mass of polymers and with a rate constant k 2 . By conservation of mass, the evolution of the monomer concentration iṡ The evolution of the number concentration P (t) = j≥nc f j (t) and mass concentration M (t) = j≥nc jf j (t) can be found by summing over j in (1). After some algebra, one obtainṡ Analytical approximation [7,14,15] of the system of equations gives a solution that depends on two parameters, λ = 2k + k n m(0) nc and κ = 2k + k 2 m(0) n2+1 . We fit our data with the form given in Eq. 1 of Ref. [7], using a least squares method. Each curve is fitted independently. Diffusion plays an important role in the aggregation progress, shifting the aggregation curves as shown in Fig. 2a and Fig. 2b. For a considerable parameter range, however, the time evolution of the fractional polymer mass can be fitted by mean-field theory (lines in Fig. 2a and Fig. 2b) with effective parameters that now depend on the diffusion rate k D (see Fig. 2c and 2d). Similarly, mean-field theory describes the density dependence of the aggregation curves as shown in Fig. 2b.

FLUCTUATIONS IN PROTEIN AGGREGATION
Having confirmed that our computational model faithfully reproduces polymerization kinetics in macroscopic samples, we now turn to the main focus of the paper, the study of sample-to-sample fluctuations in small volumes, a feature that can not be studied with mean-field kinetics. When the sample volume is reduced, we observe increasing fluctuations in the aggregation kinetics as shown in Figs. 3a and 3b. These results are summarized in Fig 3c showing the complementary cumulative distributions of half-times for different monomer numbers N and constant monomer concentration where P (x) is the probability density function and t 1/2 is defined as the half-time of the polymerization curve (i.e. the time at which M/M 0 = 1/2). The steepness of the aggregation curves in Figs. 3a and 3b suggests that, for k S k M , fluctuations are mostly ruled by the time of the first primary nucleation event, t 0 , whose complementary cumulative distribution S 0 (t 0 ) can be estimated analytically as where f M is the average number of possible primary nucleation events per unit monomer. We estimate that f M = 3ρ using a Poisson approximation, as we show in details in the following section. Note that Eq. 6 displays a size dependence that is reminiscent of extreme value statistics S 0 (x, N ) = exp(−N F (x)) where F (x) is a function that does not depend on N [22,23]. If k S k M , the half-time t 1/2 differs from the nucleation time by a weakly fluctuating time τ . This comes from the observation that, once the first primary nucleation event has happened, the polymerization follows rapidly, thanks to fast growth and secondary nucleation. This yields a weakly fluctuating delay τ (N, ρ) = t 1/2 − t 0 , which in general depends on the number density ρ and on the number of monomers N . The distribution and average values of τ are reported in Fig. 7. The average value of τ decreases with ρ and displays only a smaller dependence on N . The distribution of τ is always peaked around its average but while at small values of ρ the peak shifts with N while the standard deviation remains constant, for higher values of ρ only the standard deviation depends on N and the peak position does not change. Since the fluctuations in τ are much smaller than the fluctuations of t 0 we can safely assume that t 1/2 t 0 + τ for t 0 ≥ 0, so the complementary cumulative distribution takes the form The prediction of Eqs. 6 and 7 are in agreement with numerical simulations results for S 0 (t 0 ) and S(t 1/2 ), respectively (see Figs. 3c and 4). In particular, the behavior of S 0 (t 0 ) is obtained without any fitting parameters, while S(t

THEORETICAL DERIVATION OF THE HALF-TIME DISTRIBUTION
In this section, we provide a detailed derivation of Eqs. 6 and 7 in the limit of relatively large diffusion when the system is well mixed. To this end, we consider a cubic lattice composed by m 3 0 nodes, in which N monomers are placed randomly at time t = 0. As illustrated in Fig. 8, each monomer i sits near l (i) neighboring monomers and 6 − l (i) neighboring empty sites, where l (i) is in general a fluctuating time dependent quantity. In the model, each monomer i can either diffuse into one of the 6 − l (i) empty sites or form a dimer with one of the l (i) neighboring monomers. Therefore, at any given time t the number of possible diffusion events in the system is n D (t) = i (6 − l (i) ) and the number of possible aggregation events is n M (t) = 1/2 i l (i) , where the factor 1/2 is needed to correct for the double counting of the number of monomer pairs. We can compute the time of first aggregation of N monomers using Poisson statistics, considering for simplicity the case in which the number of possible aggregation events n M would not depend on time. In this case, the probability of having an aggregation event within ∆t is n M k M ∆t. We can then divide the time interval t 0 in n elementary time subintervals ∆t = t0 n , the rate of aggregation event at t 0 , i.e. the probability per unit time to have the first dimer formed after a time interval t 0 has elapsed, is given by the following expression:P As stressed previously, n M and n D are in principle fluctuating quantities and therefore Eq.8 is not valid. Yet, as shown in Fig. 9: n M and n D are are both i) stationary, ii) ergodic, iii) weakly fluctuating and iv) linearly dependent on N , on average. Hence, the probabilityP 0 (t 0 ) for a monomer to form a dimer at t 0 can reasonably be approximated by its ensemble average where we have replaced n M by its average value n M and defined f M ≡ n M N . From Eq.(9), we easily obtain the complementary cumulative distribution function recovering Eq. 6. To conclude our calculation, we still need to evaluate f M . To this end, we perform a discrete enumeration of the possible configurations of a single monomer, in the spirit of cluster expansions for percolation models. In particular, the six relevant configurations for a single monomer in contact with other monomers are reported in Fig. 8. The weight p l of a configuration in which a monomer has l occupied neighbors is assumed to be given by the binomial distribution This single particle picture suggests that the average number of primary nucleation events per monomer f M , corresponds to the average number of nearest neighbors l , divided by a factor 2 since any nucleation event encompasses 2 particles. With a similar reasoning, we estimate f D = 6 − l , i.e. the average number of empty directions. Then, from the binomial distribution (11) we get f M = 6ρ 2 and f D = 6(1 − ρ). Using these values in Eq. 6 and Eq. (7), we obtain agreement with numerical simulations as illustrated in Figs. 5 and 6 (panels (a) and (b) respectively) for different values of ρ.
Finally, we calculate the averages of the first aggregation time and the half-time as t 0 = ∞ 0 dt 0 t 0 P 0 (t 0 ) and t 1/2 = ∞ 0 dt 1/2 t 1/2 P (t 1/2 ), where P (t 1/2 ) = − dS(t 1/2 ) dt 1/2 . The expression for t 0 is given by In Fig. 4b we show the perfect agreement of the theoretical estimate given by Eq. (12) with the numerical values of the average time for the first primary nucleation event as a function of 1/N , for several densities. Notice that no fitting parameters are involved. The average time of t 1/2 follows from t 1/2 = t 0 + τ : the inset of Fig. 3c confirms the agreement between the theoretical estimates (dashed lines) and the numerical values (symbols).

STATISTICAL ANALYSIS OF SEED DETECTION TESTS
While the fluctuations we observe are intrinsic to the random nature of nucleation events, the ones usually encountered in bulk experiments are likely due to contamination or differences in initial conditions [3,24,25]. In those bulk systems (µl and larger), the number of protein molecules involved in the aggregation process is extremely large, even at low concentrations, so that we can exclude intrinsic kinetic fluctuations. For instance, a volume of 100µl at a concentration of 1µM still contains 10 14 monomers, leading to a large number (hundreds to thousands) of nucleation events per second for a realistic value of the nucleation rate [7,26]. However, if the relevant volumes are made small enough (pico-to nanolitres), the stochastic nature of primary nucleation can be directly observed. This has been exploited by aggregation experiments performed inside single microdroplets, where individual nucleation events could be observed, due to their amplification by secondary processes [13]. In these experiments, the average half-time is found to scale with volume in a similar manner to what is shown in the inset of Fig. 3, thus in agreement with our simulations.
We are now in a position to use our model to design a test in silico to detect the presence of pre-formed polymers, that act as seeds and nucleation sites for the secondary nucleation process, and that are thus amplified. As illustrated in Fig. 10a, the test considers a set of small volume samples containing protein solutions with a given concentration at time t 0 = 0. The aim of the test is to detect the samples containing at least one seed (case B in Fig. 10a). In an ideal experiment, the size of the microdroplets would be adjusted such that most droplets contain no seeds, some contain one seed, and the proportion of droplets containing more than one seed is negligible. In practice, these conditions can be easily adjusted experimentally by progressively reducing droplet volumes until only a small fraction of them display aggregates. After a fixed time t 1 ∼ τ , one can observe which samples contain macroscopic, detectable amounts of aggregates, enabling exact quantification of the number of aggregates present in the initial sample. Ideally the test should be positive only when at least one seed was initially present, but given the large fluctuations intrinsic to the nucleation processes we demonstrated above, as well as the competition with de novo nucleation, there is a chance for false tests. In particular, a false positive test occurs when an unseeded sample is found to contain aggregates, while a false negative test corresponds to the case in which a seeded sample does not produce detectable amounts of aggregates within the time scale of the experiment.
In Fig. 10b we report the complementary cumulative distribution of aggregation half-times t 1/2 as a function of the primary nucleation rate k M for samples with or without seeds, in this case a single pre-formed trimer. For small values of k M , seeded and unseeded samples yield distinct results, as also illustrated by the average half-times reported in Fig. 10c. As the value of k M increases, however, the distributions become closer in the two cases. In Fig. 10d, we quantify the fraction of false positive and false negative tests for two different testing times (e.g. t 1 = 500 and t 1 = 600). As expected, for large k M errors are very likely and the test would not be reliable. For intermediate values of k M , one can try to adjust t 1 to reduce possible errors with a caveat: decreasing t 1 reduces false positive errors, but at the same time increases false negatives (Fig. 10d). It is, however, possible to optimize t 1 so that both types of errors are minimized. In an experimental realisation of such a setup, the most important system parameter that needs to be optimized for any given protein is the ratio of secondary nucleation rate to primary nucleation rate. Due to potentially different dependencies of these two rates on the monomer concentration [27], pH [2] and potentially other factors, such as temperature, salt concentration etc., it is possible to fine tune this ratio and adjust it to a value that allows for an easy discrimination between droplets that do contain a seed aggregate and those that do not.

CONCLUSIONS
In conclusion, we study protein polymerization in a three dimensional computational model and elucidate the role of protein diffusion in the polymerization process. Most theoretical studies of protein aggregation neglect completely the role of diffusion and any other spatial effect. When the polymer diffusion and elongation rates are large enough we recover the standard polymerization curves that can also be obtained from mean-field analytical treatments and that can be used to fit, for example, kinetic data of amyloid β aggregation [7]. It would be interesting to explore if for small diffusion rates and small densities mean-field kinetics would eventually fail to describe the results, but this is a challenging computational task. At low densities, diffusion could play an important since a critical timescale would be set by the time needed by two monomers to meet before aggregating. This time-scale can be estimated considering the time for a monomer to cover a distance x D ∼ ρ −1/3 , yielding t D ∼ x 2 D /D ∼ D −1 ρ −2/3 . This timescale is not relevant for our simulations since at the relatively high densities we study a considerable fraction of monomers are close to at least another monomer (see n M in Fig. 9a). Consequently, the distribution of the first aggregation time does not depend on the diffusion rate k D (see Eq. 6). The half-time distribution, however, depends on diffusion even in this regime (Fig. 2).
Our simulations show intrinsic sample-to-sample fluctuations that become very large in the limit of small volumes and low aggregation rates. We show that the corresponding half-time distributions are described by Poisson statistics and display size dependence. As a consequence of this, the average half-times scale as the inverse of the sample volume, in agreement with insulin aggregation experiments performed in microdroplets [13] and with calculations based on a master-equation approach [16]. We use this result to design and validate in silico a pre-clinical screening test based on a subdivision of the macroscopic sample volume that will ultimately allow the determination of the exact number of aggregates that was initially present. This is the first step to develop microarray-based in vitro tests for early diagnosis of neurodegenerative diseases.