A mathematical model of the dynamics of prion aggregates with chaperone-mediated fragmentation

Prions are proteins most commonly associated with fatal neurodegenerative diseases in mammals but are also responsible for a number of harmless heritable phenotypes in yeast. These states arise when a misfolded form of a protein appears and, rather than be removed by cellular quality control mechanisms, persists. The misfolded prion protein forms aggregates and is capable of converting normally folded protein to the misfolded state through direct interaction between the two forms. The dominant mathematical model for prion aggregate dynamics has been the nucleated polymerization model (NPM) which considers the dynamics of only the normal protein and the aggregates. However, for yeast prions the molecular chaperone Hsp104 is essential for prion propagation. Further, although mammals do not express Hsp104, experimental assays have shown Hsp104 also interacts with mammalian prion aggregates. In this study, we generalize the NPM to account for molecular chaperones and develop what we call the enzyme-limited nucleated polymerization model (ELNPM). We discuss existence, uniqueness and stability of solutions to our model and demonstrate that the NPM represents a quasi-steady-state reduction of our model. We validate the ELNPM by demonstrating agreement with experimental results on the yeast prion \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[$$\end{document}[PSI\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^{+}]$$\end{document}+] that could not be supported by the NPM. Finally, we demonstrate that, in contrast to the NPM, the ELNPM permits the coexistence of multiple prion strains. Electronic supplementary material The online version of this article (doi:10.1007/s00285-015-0921-0) contains supplementary material, which is available to authorized users.


Introduction
The central dogma of molecular biology stipulates that phenotypes, an organism's expressed states, are determined by genotypes, the vertically transmitted DNA (Crick 1958). However, the link between genotype and phenotype is not always this direct. Today we understand that a number of phenotypes are determined epigenetically, without a change to the nucleotide sequence of DNA (Goldberg et al. 2007). In 1965, a number of yeast phenotypes were found to violate the laws of Mendelian inheritance and were thus inconsistent with DNA-based transmission (Cox et al. 1988). Further experimental studies demonstrated that the phenotypic states were not the function of the underlying DNA but were the function of a misfolded (prion) protein (Wickner 1994). As such, the phenotypes were transmitted by the proteins themselves. This phenomenon of "protein only inheritance," also called the prion hypothesis, has over time gone from highly controversial to commonly accepted (Tuite and Serio 2010). Today, nearly a dozen proteins in yeast have been shown to be able to behave as prions (Liebman and Chernoff 2012). Of course, prions extend far beyond yeast. In mammals, prions are associated with a number of irreversible fatal neurological diseases such as Creutzfeldt-Jakob disease, fatal familial insomnia, chronic wasting disease and bovine spongiform encephalopathy. Mammalian prion diseases have varying modes of transmission and have been shown to be able to pass from one species to another. At present, all mammalian prion diseases are the result of a single protein, PrP (Aguzzi and Polymenidou 2004). In addition, prion diseases are closely related to other protein misfolding diseases such as Parkinson's, Huntington's, Alzheimer's and diseases (Brundin et al. 2010;Knowles et al. 2014;Brettschneider et al. 2015).
Although humans are vastly different from yeast, the dynamics of prion proteins in both hosts is quite similar. Both mammals and yeast have cellular machinery dedicated to identifying and removing misfolded proteins (Nelson et al. 2008)-prion proteins are capable of evading such protective mechanisms and transmitting their misfolded (prion) state to other normally folded proteins. Prion proteins aggregate into complexes which act as templates for initiating further misfolding of normally folded protein.
These aggregated complexes may also fragment into smaller units, each of which can template further misfolding (Sindi and Serio 2009;Tuite and Serio 2010). Finally, in order to spread the prion state to a colony or throughout a tissue, prion aggregates must be transmitted to other cells. In yeast colonies, prion aggregates are transmitted from mother to daughter cells during cell division (Tuite and Cox 2003;Derdowski et al. 2010). In mammalian prion diseases, PrP aggregates are thought to be transmitted extracellularly (Collinge 2001). For other mammalian neurodegenerative diseases, there is a growing body of evidence suggesting neuron to neuron propagation of the misfolded proteins (Brettschneider et al. 2015).
Many mathematical models have been developed to study the dynamics of prion aggregates primarily in the context of the mammalian host (Masel et al. 1999;Prüss et al. 2006;Greer et al. 2006;Calvez et al. 2010). Tanaka et al. (2006) applied these mathematical models to the [PSI + ] prion in yeast. However, experimental studies of [PSI + ] have shown that the molecular chaperone Hsp104 is essential for fragmentation, and recent studies have demonstrated that Hsp104 acts in a rate limiting fashion with respect to fragmentation (Satpute-Krishnan et al. 2007;Derdowski et al. 2010). We also note that fragmentation is important not only for efficient conversion of normal protein by providing more templates, but also to ensure there are sufficiently many templates to allow efficient transmission; thus an accurate model of fragmentation is essential to understanding the in vivo dynamics of prion aggregates. As such, to accurately model prion aggregates in yeast, the dynamics of Hsp104 and its interaction with prion aggregates must be considered.
Further, modeling Hsp104 will lend insight to more general prion and protein misfolding disorders. Although no chaperones are known to be involved in mammalian prion dynamics, prion amplification in mammals necessarily requires aggregate fragmentation (Masel et al. 1999). In addition, while mammals do not express Hsp104, recent in vitro work demonstrates that engineered mutants of Hsp104 suppress the toxicity of misfolded protein aggregates associated with mammalian neurodegenerative disorders (Jackrel and Shorter 2014). At present, no mathematical model exists which considers the dynamics of protein aggregates in the presence of a chaperone mediating fragmentation.
In this study we develop a mathematical model of prion dynamics where fragmentation requires the interaction of Hsp104 with aggregated proteins. In Sects. 2 and 3 we provide the mathematical background and analysis of our model, which we call the enzyme-limited nucleated polymerization model (ELNPM). In Sect. 4 we illustrate the necessity of including enzyme-limited fragmentation by demonstrating important experimental properties of [PSI + ] that are not described by previously published mathematical models. We also demonstrate that in contrast to models which consider only the prion aggregates, interactions with the enzyme Hsp104 permit stable co-existence of multiple prion strains. In Sect. 5 we provide a summary and concluding remarks.

Mathematical models of prion aggregate fragmentation
We develop our model of enzyme-mediated fragmentation by considering the key biochemical processes involved in the dynamics of prions. We first discuss the dynamics included in previous mathematical formulations and then detail the additional features necessary to depict interactions between enzymes and aggregates. Finally, we demonstrate that through a series of assumptions, consistent with the yeast prion [PSI + ], our system of infinite ordinary differential equations can be analyzed with a 5-dimensional system of differential equations which approximates the full system dynamics.

Prion aggregate dynamics
While the biochemical processes depicted differs between prion model formulations (Masel et al. 1999;Prüss et al. 2006;Greer et al. 2006;Calvez et al. 2010), in all cases aggregates change size through conversion and fragmentation. That is, a prion aggregate increases in length by actively converting and incorporating normal (healthy) protein monomers. Typically, aggregates are assumed to be linear fibrils and, as such, conversion of normal protein can only take place on one of the fibril ends. Aggregates may also fragment into two smaller aggregates, each of which now act as a template to convert additional protein. It is often assumed that any aggregates smaller than the minimum stable size, n 0 , immediately disassociate into healthy prion monomers (see Fig. 1).
Such models are referred to as nucleated polymerization models (NPM); mathematical formulations of the NPM were first introduced and subsequently validated by Nowak et al. (1998) and Masel et al. (1999). This model is so-named due to the assumption that there is a minimum "stable" size of a prion aggregate (a nucleus). The spontaneous formation of such an initial nucleus (or seed, as it is also called) is the time-limiting step in prion disease initialization, but once seeded, the disease progresses primarily by the processes of conversion, fragmentation and transmission.
The NPM equations are derived from the Law of Mass Action applied to a minimal set of kinetic rate equations. Masel et al. (1999) give them as where s(t) denotes the concentration of healthy protein and u m (t) the density of aggregates of size m. Many authors (Greer et al. 2006;Prüss et al. 2006;Engler et al. 2006) have studied the more analytically-tractable equations that come from a continuous relaxation of aggregate sizes: Fig. 1 Nucleated polymerization model: conversion and fragmentation (n 0 = 2). Conversion of healthy protein (circles) lengthens the aggregate (squares), which may in turn fragment. If a daughter fragment is smaller than the stable nucleus size (n 0 ), it is immediately disassociated into healthy protein monomers Though it is known that this latter system converges weakly to the former in the limit of large average aggregate size under very general assumptions (Doumic et al. 2009;Doumic and Gabriel 2010), we choose to generalize the discrete model for simplicity.
In adapting this prion model to yeast, we identify the kinetic parameters as representing the following physical quantities: • α s the basal rate of transcription of Sup35, • μ s the decay (or dilution) rate of Sup35, • n 0 the minimum stable aggregate size (aggregates of size smaller than n 0 immediately disassociate into soluble Sup35), • μ 0 the decay (or dilution) rate of aggregated protein, • β the rate of conversion of healthy protein (from the end of a prion filament), and • γ (m − 1) the fragmentation rate of a prion aggregate of size m.

Enzyme-mediated fragmentation
We now draw attention to an underlying assumption in these equations that we will modify: for the mammalian prion PrP, it was assumed that the fragmentation rate is an intrinsic function of the aggregate size itself. However, with yeast prion systems, it has been demonstrated that fragmentation requires the additional presence of heat-shock protein 104 (Hsp104) (Satpute-Krishnan et al. 2007). Its under-or overexpression can eliminate prion aggregates entirely (Chernoff et al. 1995). Additionally, over-expressing Sup35 results in a translational shift in the aggregate-size density (Derdowski et al. 2010)-the model equations of Masel et al. (1999) do not admit such behavior (see Sect. 4.3), suggesting the possibility of a rate-limited fragmentation mechanism rooted in the Hsp104 interactions.
Though yeast prion systems have been studied with the NPM (Tanaka et al. 2006), we believe the impact of Hsp104 to be nonnegligible and explicitly consider the Hsp104 concentration and dynamics. We assume a prion aggregate of size i has i − 1 sites to which a hexamer (the active unit) of Hsp104 can bind and subsequently fragmentwe denote such an aggregate with j bound hexamers as X i, j . We note that there are actually i−1 j unique configurations of bound Hsp104 that X i, j could refer to, but the proposed kinetic equations will best be described by the amount of Hsp104 bound, not their configuration. We use standard terminology for the enzyme kinetics (parameters k on and k off ), and for simplicity, we do not model the formation of Hsp104 hexamers from monomers explicitly. We additionally define α h and μ h for Hsp104 similarly as α s and μ s for Sup35. Lastly, while we write the dilution rates μ s , μ h , and μ 0 separately, we will take them as all equal to the rate of growth of the yeast cell in our numerical experimentation. We now propose our generalization in the form of the following kinetic relations (and illustrate in Fig. 2).

Enzyme Kinetics
Hsp104 Fragmentation (with unspecified probability κ(m, n; i, j)) The key unknown in our model is the density of configurations of bound Hsp104, which is incorporated into the fragmentation kernel κ (m, n; i, j). Typically, all fragmentation sites in an aggregate are taken to be equally likely (Masel et al. 1999;Prüss et al. 2006;Greer et al. 2006): Furthermore, we require total Sup35 and Hsp104 to be conserved across fragmentation events, which corresponds to We claim that is the natural choice, which follows from taking each i−1 j configuration of X i j to be equally likely, which in turn corresponds to enzyme binding acting on a faster timescale than conversion and fragmentation (refer to the Supplemental Materials for the argument).

Enzyme-limited nucleated polymerization model
With the biochemical processes defined, we are now able to formally derive our ELNPM. We define s(t) as the concentration of soluble Sup35, η(t) as the concentration of aggregates, z(t) as the concentration of bound Sup35 (z(t) ≥ n 0 η(t)), h(t) as the concentration of unbound Hsp104, and z b (t) as the concentration of bound Hsp104. Using [S] to denote the concentration of chemical species S, these definitions correspond to and apply the Law of Mass Action to our proposed kinetic equations. We obtain We note that u i j /η defines a probability mass function for the joint random variable thus the sums may be interpreted as an expected value. In the interest of obtaining a simple set of equations, we approximate and The error of this approximation is on the order of the inverse square of the average aggregate size, a size which is typically large by assumption Greer et al. 2006;Doumic et al. 2009). (This analysis is provided in the Supplemental Materials.) Making these substitutions, we obtain the enzyme-limited, nucleated polymerization model (ELNPM): and marginal density equations (13) We note that these equations may alternatively be derived by prescribing the bino- , then finding the unique p(t) that preserves Hsp104 conservation across fragmentation events.
We further observe that the systems (s, {u m }) and (s, η, z) are identical to that of the original NPM (Masel et al. 1999), except the constant fragmentation rate γ has been replaced by the time-varying quantity γ p(t). This provides the interpretation of p(t) as the measure of how effectively the system is fragmenting at a given time. Furthermore, this formulation suggests a quasi-steady-state interpretation of our approximation, since now every aggregate is always bound with Hsp104 proportionally as p(t); that is, the enzyme binding reaches equilibrium before any conversion or fragmentation events occur. We finally note that the first 2 terms in Eq. (12e) reflect a Michaelis-Menten simplification of the enzyme kinetics (Segel and Slemrod 1989); however, since Hsp104 off-binding results in a change in the amount of binding substrate with probability γ /(γ + k off ), we may view the last term as the correction to preserve Hsp104 conservation.
Before detailed analysis of the ELNPM, we briefly examine the qualitative form of the aggregate size distribution. In Fig. 3 we plot a typical equilibrium solution to Eq. (13), where we have defined x m = u n 0 −1+m /η to be the corresponding probability mass function over the natural numbers. As expected, given the asymptotic similarity  2006), appropriately modified to match the steady-state, effective fragmentation rate with the paper's constant rate of our system to (Masel et al. 1999) (in that γ p(t) presumably converges to a fixed γ ), the equilibrium size density is of the same distribution.

Analysis of the ELNPM
We first prove a few results on existence and uniqueness for this system and then provide a non-dimensionalized, transformed system we will use to study stability. We analytically demonstrate the stability of the disease-free state and derive conditions which will ensure aggregate persistence.

Existence and uniqueness of solutions
Theorem 1 Trajectories of Eq. (12) remain invariant under a bounded, "feasible" subset of the non-negative cone R 5 + .
Proof With the exception of the η z−η term in p , the derivatives are polynomial in the dependent variables, which yields continuous partial derivatives. Considering now this last term, z ≥ n 0 η in our feasible region so we have 0 ≤ η z−η ≤ 1 n 0 −1 ≤ 1. This term's partial derivatives are also continuous in this region; let q(η, z) = η/(z − η). Then, z/η = 1 + 1/q and and Thus, q, q η , q z are continuous in our region, even as η → 0 + . This establishes existence and uniqueness.

Non-dimensionalized equations
We now reduce the ELNPM equations to non-dimensional form, where we have replaced z(t) by the displacement of the average aggregate size from the minimum size ω(t) = z(t)/η(t) − n 0 ≥ 0. We have scaled time by γ , s(t) and η(t) by α s /μ s and h(t) by α h /μ h , and used the following non-dimensional constants All subsequent analysis will done with respect to these non-dimensional equations. We note that by construction k −1 > 1 and by assumption, B > 1.

Asymptotic behavior of ELNPM
With the nondimensional equations established, we next consider the asymptotic behavior of the ELNPM. We call any trajectory satisfying lim t→∞ η(t) = 0 diseasefree; otherwise, we call the prion state persistent. Prüss et al. (2006) observed that an appropriate transformation could reduce the NPM equations to the standard SEIS model of mathematical epidemiology-a model which is governed entirely by a single parameter R 0 (the basic reproductive number). If R 0 < 1, the only equilibrium is disease-free and is globally stable. If R 0 > 1, a unique endemic equilibrium appears and exchanges stability with the disease-free state; that is, the endemic equilibrium is asymptotically globally stable (Li et al. 1999). Applying the transformation from Prüss et al. (2006), our R 0 will vary in time through its dependence on p(t): It is convenient to think of the quantity R 0 ( p(t)) = R t as the effective reproductive number of the disease system, where typically R t < R 0 . Though p(t) appears to always converge to a fixed steady-state value, it does so in a non-trivial way making it difficult to provide a Lyapunov analysis. Instead, we provide a local analysis of the disease-free state and offer numerical evidence in support of R 0 ( p disease-free ) determining global stability.

Disease-free steady state
At a disease-free equilibrium, we have η = 0 and subsequently, s = h = 1. Thus, we need only study solutions to This system has five solutions in general, though we shall show there is only a single solution inside our feasible region.
We now establish local stability criteria of this root; let p 0 and ω 0 be the unique solution to Eq. (19).
Proof The eigenvalues of the localized Jacobian will satisfy Of the five roots, the first 2 are clearly negative. The quadratic factor will also admit 2 stable roots: we see that its quadratic and linear coefficients are strictly positive and now show that the constant term is as well. If k 1 > 0, then k 1 +k −1 −2 p 0 > 2(1− p 0 ) > 0.

Endemic steady state
The local instability of the disease-free equilibrium yields a persistent disease state; however, the numerical experimentation in Fig. 4 suggests further that there is an attractive, endemic equilibrium. Generally speaking, steady-state solutions to our system will be solutions of a quintic polynomial in five variables, thus preventing closed-form descriptions of such states. However, one can parameterize these values in terms of a fixed (but unknown) valuep corresponding to p =p: wherep satisfies We note that this quadratic is uniquely invertible within our feasible region, which yields the recursive relation We draw attention to the similarity between Eq. (27) and what would be the Michaelis-Menten value forp. Our feasible region requires s, h, η, ω > 0;s andω are always positive, whilẽ η > 0 ⇒s < 1 ⇒ R 0 (p) > 1. This also gives ush > 0, since we typically Thus, we will have endemic equilibria when Eq. (26) has solutions satisfying R 0 (p) > 1. Based on considerable numerical studies, we conjecture that this solution uniquely exists and is globally asymptotically stable when R 0 ( p 0 ) > 1.

Discussion
Formulating the ELNPM allows us to consider aspects of prion aggregate dynamics that cannot be explained by prior mathematical approaches that neglected the role of the Hsp104 chaperone in fragmentation. We first demonstrate that the NPM is a limiting case of the ELNPM and comment on implications of the ELNPM to the larger question of the appearance of prion strains in a population. We then demonstrate that the ELNPM is the first model capable of supporting two experimentally observed phenomena. First, the ELNPM is the first model capable of reproducing shifts in the aggregate densities associated with increases in synthesis of Sup35. Finally, we demonstrate that the binding kinetics of Hsp104 in the ELNPM allows the possibility of multiple co-existing prion strains. Intriguingly, the co-existence of multiple strains is thought to be crucial towards understanding the transmission of prion diseases between species (Chien et al. 2004).

The NPM is a limiting case of ELNPM
Our moment-closed ELNPM model, Eq. (9), is nearly identical to the original NPM model but with a time-varying, effective fragmentation rate γ p(t) instead of the constant γ . When p(t) →p, the dynamics of ELNPM will mirror that of the NPM with γ replaced by γp. As such, it is convenient to think of the NPM as a quasi-steady-state approximation of the full enzyme kinetics we have considered in our model. We informally used this observation in Sect. 3.3 to motivate (but not prove) global stability based on known results of the NPM system.
We plot in Fig. 5 p(t) over the same parameters in Fig. 4 and note that-since yeast has a doubling time of roughly 90 min (Hartwell and Unger 1977) p(t) will not reach its asymptotic value for a few cell-divisions. As such, even though the NPM represents  Fig. 4. The transient fragmentation efficiency may be higher (for small R 0 ) or lower (for large R 0 ) than the asymptotic efficiency a quasi-steady-state approximation it may not represent the aggregate dynamics during the early cell divisions following the introduction of prion aggregates.

Transient fragmentation efficiency may impact prion stability
Although the NPM may be viewed as an asymptotic simplification of the ELNPM, we remark that the differences in transient behavior may provide insight into the underlying stochastic dynamics that arise when a single prion seed is introduced into a healthy yeast colony. The ELNPM predicts an initial fragmentation rate that can be larger or smaller than the asymptotic rate-this is because the availability of enzyme (Hsp104) is much larger than the availability of substrate (binding sites) in this initial configuration (see Fig. 5).
Since aggregate amplification is essential to spreading a prion disease these transient fragmentation rates may impact prion stability. For example, a higher transient fragmentation rate for a prion strain with low R 0 would represent a barrier to successful "seeding" of the prion state than would otherwise be predicted by a constant fragmentation rate. This provides a plausible mechanism for the removal of an initial prion aggregate appearing in a population and therefore the low frequency of spontaneous appearance of the prion state.

Hsp104 acts as a rate limiter for fragmentation
We noted in Sect. 2 that the NPM does not admit translational shifts in the aggregate density as a function of increasing synthesis of the prion protein (α s ), we now formally demonstrate this is the case. Let us revisit the quantity u m (t), the density of Both the scaling and translation are affected by α s . Initial kinetic parameters are as in Fig. 3 aggregates of size m. We write x m = u m+n 0 −1 /η-this quantity defines a probability mass function that is independent of the amount of aggregated protein. In our rescaled variables, we have Davis and Sindi (2015) gave a closed-form for x m at steady-state: where ζ = μ γp + n 0 . The size distribution's dependence on α s can only occur through its relationship with the steady-state fragmentation efficiencyp. This is fixed in the NPM, thus the size distribution will not change in response to changes in α s . This is in contradiction to the experimental results described by Derdowski et al. (2010). Since our model does allowp to vary as a function of the kinetic parameters, we are able to numerically investigate qualitative shifts in the distribution. We demonstrate these shifts in Fig. 6, which are in qualitative agreement with the experiments of Derdowski et al. (2010). As such, experimental evidence supports that fragmentation can not be purely a function of the number of available fragmentation sites and must depend on the amount of Hsp104 in the system.

Prion extinction and Hsp104 expression levels
Beyond translational shifts in aggregate size distribution, the [PSI + ] prion phenotype in yeast has been shown to be very sensitive to the amount of Hsp104 in the system Fig. 7 Hsp104 production is up-regulated or Hsp104 is deactivated after 3 h, both by a factor of 10 4 . The new system in either case is unable to stably support the presence of prion aggregates with our engineered parameters (Higurashi et al. 2008;Doyle and Wickner 2009;Shorter and Lindquist 2008). Our mathematical formulation correctly captures the dependency of all prions to the underexpression of Hsp104. In contrast, and in agreement with recent experimental studies (Klaips et al. 2014), over-expression of Hsp104 does not necessarily drive prions to extinction.
First, sufficiently high concentrations of guanidine hydrochloride GdnHCl have been shown to severely disrupt the fragmentation process by inactivating Hsp104 (Ferreira et al. 2001;Byrne et al. 2007). Since fragmentation is halted the total number of aggregates will not change and aggregates will eventually be lost through dilution in the population due to cell division. Over time the population will be cured of the prion disease as the fraction of cells with aggregates approaches zero.
Quantitatively, we treat the inactivation of Hsp104 as letting k on → 0. In the limit, we'll obtain p = −k 1 p + p 2 − Bsp which corresponds to the elimination of prion aggregates. We demonstrate this in Fig. 7.
While earlier studies seemed to indicate over-expression of Hsp104 could cause loss of the prion state (Chernoff et al. 1995), more recent experimental evidence indicates that Hsp104 over-expression is not sufficient to drive prions to extinction (Klaips et al. 2014). In our formulation over-expressing Hsp104 drives p → 1. This is readily demonstrated by assuming k 1 h = 1/ 1 and rescaling time by this quantity; then p = 1 − p + O( ). However, as is experimentally, this alone is not mathematically sufficient to cure the prion state.
Consider an endemic state with R 0 (p) = B/p Depending on the other kinetic parameters, R 0 may be either increasing or decreasing with respect Fig. 8 The reproductive number R 0 as a function of α h . Prion strains will only be driven to extinction by Hsp104 over-expression if lim α h →∞ R 0 (α h ) < 1 top-prion extinction would only occur if R 0 (1) < 1. Specifically, sgn(R 0 ( p)) = sgn(A 2 0 − (n 0 − 1)n 0 p 2 ). If A 0 > n 0 − 1, R 0 is always increasing and extinction is impossible under our model. This is consistent with the belief that there are other (unmodeled) factors more likely to contribute to prion phenotype loss (Palmer et al. 2011;Klaips et al. 2014); nonetheless, we do provide an engineered parameter set (described in Table 1) that demonstrates prion extinction by maximizing fragmentation efficiency in Fig. 7. We additionally plot the dependence of R 0 on α h in this parameter set as well as the original set we've used in Fig. 8.

Stability of the coendemic prion strains
Up until this point we have considered the aggregates of a single prion species; however, a prion protein is capable of adopting a host of misfolded confirmations each of which is capable of the biochemical processes of conversion of normal protein and fragmentation (Sindi and Serio 2009;Tuite and Cox 2003;Tuite and Serio 2010). That is, [PSI + ] or PrP Sc does not refer to a single prion phenotype, but many related ones, each characterized by different pathology (implying different kinetic parameter values). These distinct prion states are referred to as prion "strains." Biologists have observed multiple coexisting strains (Strbuncelj 2009;Polymenidou et al. 2005), but there has been limited mathematical modeling of multiple prion strains. Previously, Tanaka et al. (2006) considered the NPM, under the continuous relaxation of aggregate sizes, with n 0 = 1 and demonstrated that if two strains were present then, asymptotically, one strain would dominate and drive the other to extinction. The outcome was determined by the strain which had maximized βγ (which is proportional to the reproductive number in the case of continuous-size, n 0 = 1 NPM). Since level curves of βγ represent a set of measure 0 in parameter space, realistically this prevents asymptotic prion strain coexistence. By coexistence, we mean that there exists i = j such that lim t→∞ η i (t), η j (t) > 0 when η i (0), η j (0) > 0 (where η i is the concentration of aggregates of strain i).
We now generalize the ELNPM to include multiple prion strains, each capable of converting the same normal protein. Because aggregation is based on conversion to a particular prion strain conformation, we consider an aggregate as consisting of misfolded protein of a single strain. We write the equations for k strains with similar constants as before, but scale time by Γ = k i=1 γ i and write Γ i = γ i /Γ : Because our mathematical formulation also requires the molecular chaperone Hsp104, this opens up the possibility for an alternative pathway to prion strain coexistence-rather than out-compete solely on conversion β and fragmentation γ , a second strain may be more efficient at sequestering Hsp104 (k on /k off ). Increasing this ratio improves the strain's own fragmentation efficiency, as well as decreases the other strain's efficiency by decreasing the amount of available Hsp104.
It is helpful to think of prion strain competition and coexistence in terms of the reproductive numbers described in Sect. 3.3. Intuitively, the strain with the highest reproductive number will dominate and drive others to extinction. As already stated, this is exactly what Tanaka et al. (2006) observed-however, recall that with the NPM, the reproductive number is fixed, so there will not be any dependency on the kind of or number of strains present. Our model, however, has an effective reproductive number dependent on the current fragmentation efficiency. In terms of strain-specific constants, this number is given by The fragmentation efficiency of strain i, p i , is dependent on the current concentration of soluble Sup35 and free Hsp104 (Eq. (31e)), which in turn depend on all of the strains' concentrations. As such, the reproductive numbers of the strains are coupled to one another. These nonlinear, secondary interactions make analytic determinations of coexistence difficult. However, we are able to numerically demonstrate coexistence of prion strains (Fig. 9). Remarkably, Fig. 9b demonstrates that strain coexistence is possible for parameters lying in a non-zero area of parameter space. We note that, in contrast to prior models, this type of coexistence is biologically feasible because is it robust to small perturbations in parameter space. Thus, our numerical studies demonstrate that strains with different reproductive numbers (in isolation) can coexist. Further, at the coendemic state each strain's "cooperative" reproductive number is different from their isolated value, but equal to each of the other strains' cooperative numbers.
We choose two specific parameter sets from Fig. 9b and plot of their steady-state, cooperative size densities in Fig. 10a and concentration of aggregated Sup35 in each strain over time in Fig. 10b.

Conclusion
In this work we successfully developed a mathematical formulation of aggregate dynamics where fragmentation occurs through the molecular chaperone Hsp104. We demonstrate that, under certain restrictions, our model reduces to a numerically tractable form which we call the ELNPM. By including chaperone-mediated fragmentation, this work represents an important step towards a more complete understanding of prion and protein misfolding in vivo.
We derived a unique disease-free steady-state of the ELNPM and analyzed its stability. We demonstrated that the ELNPM supports experimentally observed results such as shifts in the aggregate-size distribution with increasing Sup35 synthesis and response to over-and under-expression of Hsp104. Additionally, it represents a first step towards quantifying prion strain coexistence.
While the ELNPM successfully describes the effects of varying amounts of Sup35 and Hsp104 in the system, we note that there are factors common to enzyme-substrate kinetics that were not included in our model. First, in many biochemical systems there is evidence of cooperation between binding sites (Nelson et al. 2008). Since there is no evidence of interaction between binding sites for Hsp104, we have modeled binding events as purely a function of the free enzyme and available binding sites.
Second, by assuming that k off was large we were able to assume that for an aggregate consisting of i Sup35 monomers with j sites bound by Hsp104, all possible configurations of bound Hsp104 are equally likely (see the Supplemental Material). Since under normal expression Hsp104 is observed to be only minimally bound to [PSI + ] aggregates (Klaips et al. 2014), we interpret this to indicate that k off must indeed be large relative to k on . Third, we considered a generalization of the uniform fragmentation kernel which corresponds to equality in fragmentation at all binding sites. Together, these three assumptions allowed the use of analytical approaches previously employed in the analysis of the NPM to demonstrate existence, uniqueness and asymptotic stability of the disease-free steady-state.
Beyond Hsp104, other enzymes have been identified as important players in the dynamics of prion aggregate fragmentation (Inoue et al. 2004;Shorter and Lindquist 2008). As such, our mathematical formulation may be taken as the representing collective impact of enzymes on fragmentation. However, compared to the other enzymes, Hsp104 occurs in the lowest molecular number (Ghaemmaghami et al. 2003) and is thus likely to represent a rate limiting step. In addition, we have evaluated our model by comparison to experimental results on the [PSI + ] prion which has shown to have greatest sensitivity to Hsp104 (Higurashi et al. 2008). Lastly, again note that the form of Hsp104 is a hexamer (Doyle and Wickner 2009)-we have assumed the kinetics of hexamer formation are not relevant to the aggregate dynamics.
In addition to including additional biological complexities, in future studies we plan to investigate global asymptotic stability and explore the conditions underlying prion strain coexistence.