A unifying model for the propagation of prion proteins in yeast brings insight into the [PSI+] prion

The use of yeast systems to study the propagation of prions and amyloids has emerged as a crucial aspect of the global endeavor to understand those mechanisms. Yeast prion systems are intrinsically multi-scale: the molecular chemical processes are indeed coupled to the cellular processes of cell growth and division to influence phenotypical traits, observable at the scale of colonies. We introduce a novel modeling framework to tackle this difficulty using impulsive differential equations. We apply this approach to the [PSI+] yeast prion, which associated with the misconformation and aggregation of Sup35. We build a model that reproduces and unifies previously conflicting experimental observations on [PSI+] and thus sheds light onto characteristics of the intracellular molecular processes driving aggregate replication. In particular our model uncovers a kinetic barrier for aggregate replication at low densities, meaning the change between prion or prion-free phenotype is a bi-stable transition. This result is based on the study of prion curing experiments, as well as the phenomenon of colony sectoring, a phenotype which is often ignored in experimental assays and has never been modeled. Furthermore, our results provide further insight into the effect of guanidine hydrochloride (GdnHCl) on Sup35 aggregates. To qualitatively reproduce the GdnHCl curing experiment, aggregate replication must not be completely inhibited, which suggests the existence of a mechanism different than Hsp104-mediated fragmentation. Those results are promising for further development of the [PSI+] model, but also for extending the use of this novel framework to other yeast prion or amyloid systems. Author summary In the study of yeast prions, mathematical modeling is a powerful tool, in particular when it comes to facing the difficulties of multi-scale systems. In this study, we introduce a mathematical framework for investigating this problem in a unifying way. We focus on the yeast prion [PSI+] and present a simple molecular scheme for prion replication and a model of yeast budding. In order to qualitatively reproduce experiments, we need to introduce a non-linear mechanism in the molecular rates. This transforms the intracellular system into a bi-stable switch and allows for curing to occur, which is a crucial phenomenon for the study of yeast prions. To the best of our knowledge, no model in the literature includes such a mechanism, at least not explicitly. We also describe the GdnHCl curing experiment, and the propagon counting procedure. Reproducing this result requires challenging hypotheses that are commonly accepted, and our interpretation gives a new perspective on the concept of propagon. This study may be considered as a good example of how mathematical modeling can bring valuable insight into biological concepts and observations.


Author summary
In the study of yeast prions, mathematical modeling is a powerful tool, in particular when it comes to facing the difficulties of multi-scale systems. In this study, we introduce a mathematical framework for investigating this problem in a unifying way. We focus on the yeast prion [P SI + ] and present a simple molecular scheme for prion replication and a model of yeast budding. In order to qualitatively reproduce experiments, we need to introduce a non-linear mechanism in the molecular rates. This transforms the intracellular system into a bi-stable switch and allows for curing to occur, which is a crucial phenomenon for the study of yeast prions. To the best of our knowledge, no model in the literature includes such a mechanism, at least not explicitly. We also describe the GdnHCl curing experiment, and the propagon counting procedure.

Introduction 1
Amyloids are self-perpetuating protein aggregates, involved in many neurodegenerative 2 diseases in mammals such as Alzheimer's Disease, Parkinsons's Disease, 3 Creutzfeldt-Jakob Disease, Huntington's Disease [1]. Prions are a particular case of 4 amyloids, where a conformational change is infectious and transmissible through 5 self-catalyzed aggregation [2]. Understanding the molecular mechanisms associated with 6 amyloid growth and prion replication is crucial in the effort against neurodegenerative 7 diseases. The yeast Saccharomyces cerevisiae has emerged as a powerful model system 8 in the study of these processes [3]. Native yeast prions share many specificities with 9 mammalian prions [4], and the protein quality control machinery involved in their 10 propagation is highly conserved between mammals and yeast. Furthermore, yeast 11 models are used to screen for anti-amyloid drugs through the artificial expression of 12 mammalian proteins [5]. Many kinetic models of aggregate growth and nucleation have 13 been proposed [6], but their validation using data from yeast colonies is challenging. 14 Indeed, the propagation of yeast prions is a multi-scale system where molecular 15 processes are coupled to cellular processes such as cell growth and division to produce 16 effects that are observable at the phenotypical scale, i.e. at the scale of colonies. An 17 illustration is proposed in Fig 1. Taking into account and relating those scales calls for 18 specific modeling approaches, but to our knowledge no mathematical modeling work 19 describes the whole system rigorously and in a controlled way. We propose a novel 20 approach based on the use of impulsive differential equations. It is important to 21 emphasize that this framework is versatile and can be used to study the evolution of any 22 chemically active intracellular components through mass-action kinetics inside growing 23 yeast colonies. 24 As a first application of this approach, we focus here on the yeast prion [P SI + ]. It 25 has been studied in detail [7], and experimental results on [P SI + ] have become essential 26 groundwork for the development of prion biology [8,9]. The [P SI + ] phenotype is 27 associated with the oligomerization and loss-of-function of the Sup35 protein, an 28 enzyme release factor responsible that promotes stop-codon recognition. When 29 aggregates are formed, a change of color of the yeast colonies is observed (from red to 30 white) if they encode a stop-codon in either the ADE1 or ADE2 genes [10]. See Fig 1C 31 for an illustration of the change in phenotype. As is the case with most native yeast 32 prions, the [P SI + ] phenotype is reversible. It can be cured using different 33 physico-chemical treatments [11,12] but also by introducing point mutations in Sup35 34 or its co-chaperones [13,14]. Curing experiments are widely used to derive information 35 on the molecular properties of Sup35 aggregates, depending on the curing treatment 36 applied or the strain studied. The most exploited one is Guanidine Hydrochloride 37 (GdnHCl) treatment, because it allows to infer prion seed numbers in [P SI + ] 38 cells [15,16]. We reproduce some of those curing experiments using a model for 39 aggregate replication at the molecular scale combined with a model for cell growth and 40 division, as described in Figure 2. By coupling those aspects, we build a full model of 41 prion propagation inside yeast colonies, which is detailed in Methods. Even by keeping 42 things as simple as possible, our model sheds light onto two specific characteristics of 43 the [P SI + ] prion that were not uncovered in the past, as we detail in Results, the 44 replication of aggregates at low densities follows non-linear effects, with a kinetic barrier 45 preventing the expansion of the population of aggregates below a certain threshold. A second implication of our modeling work is that the effect of GdnHCl is not a complete 47 block of aggregate fragmentation, as was suggested before [11], rather a significant 48 reduction of the replication efficiency.

50
We develop a multi-scale model of prion replication and propagation in yeast to 51 investigate the relationship between molecular aggregation mechanisms and colony level 52 phenotypical properties. The model is presented in Fig 2 and Fig 1C. This property is exploited to track the dynamics of curing during 71 continuous treatment with a de-stabilizing agent [11,12,15]. In order to reproduce experimental results involving [P SI + ] curing, this is an essential aspect that must be 73 captured by the model.

74
A multi-stable molecular model 75 In terms of mathematical modeling, the fact that different outcomes are possible in the 76 same conditions is interpreted as multi-stability. Different stable solutions are able to 77 attract the trajectories of the system, and the outcome depends on the initial condition. 78 Our model introduces this property by including a non-linearity in the replication rate 79 of aggregates, with a sigmoidal dependency on the concentration of aggregates 80 portrayed by the function f in the model equations (see Fig 2). This shape of function 81 is frequently used to model cooperative reactions [20], it appears naturally from 82 ligand-receptor interaction equations. Without knowing the specific kinetic scheme, we 83 empirically choose a Hill function. With this non-linearity, the colony grown from a cell 84 will exhibit one of three possible behaviors depending only on its initial condition.
Daughter growth rate λ 0.7 ‡ µM.hr −1 Sup35 monomer basal production rate ρ 10 hr −1 Maximal aggregate replication rate K 1 µM Replication efficiency threshold n 5 -Replication efficiency order Simulations presented in the text use those values unless specified otherwise. See Figure 2 and Methods for a description of the equations. *: [17]. †: [18]. ‡: [19] of aggregates. de-stabilizing agents. In Fig 4, we see how reducing the proportion of aggregated Sup35 124 will bring cells into the region corresponding to a sectored phenotype. Furthermore, we 125 also know that "weak" strains are more likely to exhibit sectors [9]. Fig 4 (B) shows the 126 phenotype map for a strain that has a slightly reduced aggregate replication rate. We 127 see immediately that the region leading to sectored phenotypes is enlarged, and also 128 closer to the periodic solutions. This means that such a strain is more likely to be  Table 1. Panel (A) shows a strong strain with a maximal replication rate ρ = 10hr −1 and panel (B) is a weak strain with ρ = 0.87hr −1 sectoring is often dismissed by counting mosaic colonies as either full [P SI + ] or half 131 and half [3,11]. Perhaps some valuable information could be obtained by studying the 132 phenomenon in more detail. Among the various agents that de-stabilize [P SI + ], GdnHCl is one of the most studied. 136 It acts by impairing the chaperone Hsp104, which fragments Sup35 aggregates, thus 137 preventing the creation of new aggregates [11]. [P SI + ] colonies growing under GdnHCl 138 treatment progressively lose the prion phenotype, by dilution of the aggregates between 139 dividing cells. This property was used to infer the number of prion "seeds" in individual 140 cells [15] by fitting an exponential model to the curing curve. This approach was 141 supported by [3], with the discovery that during GdnHCl treatment, the number of 142 [P SI + ] revertants (cells that will grow into a [P SI + ] colony when transferred back onto 143 regular medium) reaches a plateau. This number was interpreted as the number of  Propagon counts and statistics on those counts are often used to characterize 147 strains [23], to quantify the bias of transmission between mother and daughter cells [17], 148 to study the effect of different point mutations [13]. It needs to be emphasized that this 149 interpretation is based on two hypotheses. The first hypothesis is that the threshold  [11,17,24]. Furthermore, [25] used fluorescent tagging to track aggregates during 168 GdnHCl treatment, and report a decrease in cells with foci slower than the halving 169 predicted by the segregation hypothesis.

170
Reproducing the propagon count experiment with our model requires very specific 171 kinetic parameters for the replication reaction. Those conditions are found numerically 172 by studying the mother-only lineage in the model. Indeed when this lineage is attracted 173 to a [P SI + ] solution but each of its daughters becomes [psi − ], we are assured that the 174 total number of [P SI + ] cells in the colony reaches a plateau. Once again this is the 175 consequence of bi-stability in the molecular model as well as the asymmetric division 176 favoring aggregate retention by mother cells. The fact that those conditions correspond 177 to a very narrow parameter range is worth emphasizing and discussing further.  Table 1). Cells are scored as [P SI + ] as soon as they contain a concentration of aggregated Sup35 higher than 0.5µM . In previous modeling work, one approach was to study aggregates as discrete 182 entities [15,26,27]. When doing so, the kinetic barrier hypothesis is implicit and the 183 threshold lies at the difference between 0 and 1 aggregate. However, it is not a proper 184 modeling approach when it comes to studying chemical reactions, because mass-action 185 kinetics cannot be applied to numbers in isolation but rather to densities. Our model 186 implements chemical modeling inside a multi-scale phenomenon in a rigorous way and 187 provides a powerful tool for generating insights into existing experimental outcomes: 188 colony sectoring and propagon counting assays. 189 Our work suggests that the multi-stability in the intracellular scale is essential to Hsp104. However, even if fragmentation is limited by the presence of Hsp104 the model 198 does not exhibit bi-stability [28]. We would like to emphasize that this type of modeling 199 (global models) is universally used in the prion biology, in particular with the 200 widespread use of the Nucleated Polymerization Model (NPM) [29,30]. The 201 introduction of bi-stable systems was suggested in earlier work [31,32], in the case of 202 mammals and more recently to study interaction between prion strains [33]. The benefit 203 of the framework we suggest here is that it could allow to use data to validate a precise 204 kinetic scheme of aggregate replication. In particular, the sectoring phenomenon should 205 be studied with a new perspective, especially during curing experiments.

206
What is the effect of GdnHCl?

207
Our results question all previous assumptions made about the effect of GdnHCl. With 208 our model, the propagon experiment is reproduced only if the chemical replication rates 209 are chosen very carefully. This is potentially explained by two reasons. The first is that 210 our model might be too simple to capture the possibilities of the full biological system. 211 This would also explain why sectoring is reduced to such a narrow region in the 212 phenotype map Fig 4.

213
A second explanation is that our model may be missing an essential chemical 214 process, which remains unaffected when GdnHCl is present. There is precedent for this 215 assumption, because there is evidence that aggregates are still chemically active under 216 GdnHCl treatment. First GdnHCl does not prevent aggregates from growing by 217 polymerizing newly synthesized Sup35 [17,34]. Furthermore, there is evidence for an 218 action of Hsp104 that is not affected by GdnHCl [35,36]. This brings up the 219 controversial question of the roles of Hsp104 in the propagation of [P SI + ] and other 220 yeast prions. The only conclusion our results bring is that GdnHCl curing is not 221 explained by an exponential dilution model as suggested previously [15,18], because of 222 the plateau of [P SI + ] cells. Our study is a first step in the design of a more elaborate 223 kinetic model that includes the effect of Hsp104.

225
We introduced a novel modeling tool to the field of yeast prions, with the major benefit 226 of relating different scales in a controlled and rigorous way. From the molecular 227 mechanisms, with a kinetic scheme built using mass-action kinetics, to the phenotypical 228 traits at the colony level, this framework has the potential of taking into account every 229 aspect of the system. In the case of the [P SI + ] prion, we build a simple model with the 230 primary goal to qualitatively reproduce experimental observations from the literature. 231 We focus on curing experiments, and reproducing any curing experiment with our 232 framework requires introducing a very particular characteristic into the molecular 233 model. The kinetic scheme needs to be bi-stable, where the prion-free state and the 234 prion state are both simultaneously stable, and the transition between the two of them 235 is a bi-stable switch. This allows the possibility of curing, and it concomitantly explains 236 the phenomenon of sectoring in a deterministic way. This phenomenon is often 237 dismissed but is in fact instructive with regards to the molecular processes. By 238 investigating in more detail the case of GdnHCl curing, we have reason to question the 239 suggested effect of this agent. Indeed this experiment is reproduced by our model, but 240 the fragmentation of aggregates must not be completely inhibited contrary to the 241 commonly accepted effect of GdnHCl.

242
Overall, the framework of impulsive differential equations is versatile and could be 243 adapted to many different cases. Studying the [P SI + ] prion already revealed 244 instructive, even though work is still in progress. In the future, we aim to use this 245 framework to build and validate a complete model of aggregate replication including the 246 role of Hsp104 and possibly its co-chaperones, a size-distribution of aggregates, 247 stochasticity in the cell division events. This would be done through hypothesis testing 248 and close collaboration with biologists. Inferring parameters is a long-term goal, that 249 first requires understanding the very structure of the molecular processes. In particular, 250 it needs to be clear what is the mechanistic origin of the cooperativity in the replication 251 of aggregates. Another use of the model is to extend it to other yeast prions and 252 amyloid models, as yeast models are used to screen for anti-amyloid drugs, and a 253 specific modeling framework would help interpreting those experiments.

255
Using impulsive differential equations to model cell division 256 We first introduce impulsive differential equations (IDEs) and use them to model a population of dividing yeast cells. Given a system of ordinary differential equations f (Z, t), where Z ∈ R n represents the concentrations of the biochemical species we are tracking, and given Z 0 ∈ R n , Z 0 ≥ 0 an initial condition, we define the impulsive system The sequence (t k ) k∈N * represents the impulsion times, which we choose to be fixed. We 257 require 0 ≤ t 0 < · · · < t k . . . and lim k→∞ t k = ∞. Each impulsion corresponds to a cell 258 division, and the vector α k ∈ R n represents the partition rate at the moment of division 259 k. A solution to this system is a piece-wise continuous function Z on (t k , t k+1 ], for 260 k ∈ N * . The existence and the properties of such solutions are studied 261 theoretically [37,38]. IDEs are used in various modeling fields, including 262 epidemiology [39], population dynamics [40] and more recently T-cell differentiation [21]. 263 In order to model cell division with IDEs, we need to clarify how species are 264 distributed between mother and daughter cells. Consider a yeast cell about to bud and 265 produce a daughter. The volume of the mother cell is V 0 , which will split into two parts 266 right after division. The mother's volume becomes V M = πV 0 and the daughter gets where π is the volume asymmetry. For yeast cells it has been observed 268 that π = 0.6 [18]. Next consider a chemical component in the mother of mass M 0 and 269 concentration C 0 = M0 V0 . During division this mass is distributed between the mother After division, the mother concentration (C M ) and the daughter concentration (C D ) are given by We dot not consider degenerate cases where ε = −π or ε = 1 − π, corresponding to 276 the case where the transmission from mother to daughter is (respectively) full or null.

277
This hypothesis is biologically relevant in the sense that it seems highly unlikely that 278 exactly no misconformed protein be transmitted at the moment of division. 279 We assume that cells are exponentially growing, and that they divide when they reach the volume V 0 . In other words, the intracellular dynamics of prion aggregation have no impact on the division cycle. This is consistent with [17]. This assumption imposes a relation on the growth rate γ and the division time T . As before, if we denote the mother parameters with a subscript M and the daughter parameters with a subscript D, we have the following two relations As a consequence when π, T M and T D are fixed (as measured experimentally [9,17]), 280 γ M and γ D are determined by these relations. Figure 2A shows a visual representation 281 of the cell division model we are using. (The parameters we used in results are given in 282 Table 1.) 283 An impulsive differential equation model allows one to follow a cell lineage and its 284 intracellular components in time. The cell history is described by the sequence of 285 impulsions it undergoes, either mother or daughter impulsions. Analytically, we study 286 in detail the two extreme lineages, mother-only and daughter-only. Indeed those 287 lineages are subject to periodic and identical impulsions, which facilitates the 288 analysis [37]. If we use the IDEs to track all lineages starting from a single cell, we can 289 study the complete yeast colony. Having formalized the cell division dynamics model, 290 we next introduce the model for intracellular dynamics.

291
A bi-stable prion replication process 292 In this section, we introduce a two-dimensional bi-stable model of prion dynamics, at first without considering the effect of impulsions. In each yeast cell, we track the concentration of soluble Sup35 (V ), and of Sup35 in the prion conformation (S). The model is illustrated in Figure 2B, and described by the following system of ordinary differential equations The parameters λ and γ are respectively the Sup35 monomer production rate and the 293 cell growth rate. For simplicity we do not include the aggregate size dynamics but 294 model the recruitment of monomers into prion aggregates. It occurs with a cooperative 295 reaction, at maximal rate ρ and non-linear efficiency f (S). The function f is chosen to 296 be a Hill function of order n(> 1) and threshold K 297 f (S) = S n K n + S n .
This function is chosen for one distinct feature: it makes the system multi-stable as 298 soon as n > 1. The prion-free equilibrium (with S = 0) exists under any choice of 299 parameters and remains locally stable. Two other equilibria appear from a saddle-node 300 bifurcation, one of them is unstable and the other is locally stable. In conditions when 301 the three equilibria exist (one prion-free and two prion equilibria), the asymptotic 302 outcome of the model depends on the initial condition. Starting with a concentration of 303 aggregates too low will cause the solution to converge to the prion-free equilibrium, 304 whereas if the initial concentration of aggregates is sufficient their population will be 305 stably maintained. See S1 Appendix for the mathematical description and proof of this 306 property, as well as some numerical illustration of the bi-stability. As detailed in 307 Discussion, it is essential that we use a model with multi-stability because we are 308 investigating a curing experiment, where cells can lose the prion phenotype.

309
Complete bi-stable model with impulsions 310 When combining the bi-stable model and the impulsive differential equation framework, 311 we obtain a complete model of aggregate replication and transmission in growing yeast 312 colonies. The choice of parameters used by default is detailed in Table 1. The cell 313 division parameters (doubling time and mother-daughter volume ratio) are measured 314 experimentally and we choose values within the typical observed range [16,17]. The 315 mass transmission bias is a parameter that will require a thorough investigation in 316 further work but we choose a reasonable value of ε = 0.1 for the time being. The Sup35 317 monomer production rate is such that the prion-free steady-state concentration of Sup35 318 in cells is λ γ ≈ 2.5µM [19]. The prion replication parameters ρ, K and n are used as 319 adjustment parameters to investigate the behavior of the model.

320
Supporting information 321 S1 Appendix. Analytical study of the bi-stable impulsive system for yeast 322 prion propagation.