Paul Lemarre
email: plemarre@ucmerced.edu

Laurent Pujo-Menjouet

Suzanne S Sindi
email: ssindi@ucmerced.edu

Generalizing a Mathematical Model of Prion Aggregation Allows Strain Coexistence and Co-Stability by Including a Novel Misfolded Species

Keywords: Coexistence, Nucleated Polymerization, Prions, Protein aggregation, Protein misfolding, Strains, Subunits Mathematics Subject Classification (2000) 34D20, 34D23, 37N35, 92C40, 92C60

HAL is

Introduction

Prion diseases are a class of neurodegenerative disorders in mammals associated with a change in the folded shape (conformation) of the protein PrP. This misfolded conformation (PrP Sc ) is believed to induce transformation of proteins in the normal conformation (PrP C ), and proteins under this abnormal state may aggregate [START_REF] Bolton | Identification of a protein that pruifies with the scrapie prion[END_REF][START_REF] Collinge | Prion diseases of humans and animals: their causes and molecular basis[END_REF][START_REF] Aguzzi | Mammalian prion biology: one century of evolving concepts[END_REF][START_REF] Tuite | The prion hypothesis : from biological anomaly to regulatory mechanism[END_REF]. Prions have also been discovered in yeast [START_REF] Tuite | Propagation of yeast prions[END_REF], where they are mostly associated with harmless epigenetic phenomena [START_REF] Lindquist | Blessings in disguise: biological benefits of prion-like mechanisms[END_REF]. The dynamics of prion propagation have been studied extensively during the past decades, using experimental cases in mammals and yeast on one hand [START_REF] Sindi | Prion Dynamics and the Quest for the Genetic Determinant in Protein-Only Inheritance[END_REF][START_REF] Tuite | Propagation of yeast prions[END_REF][START_REF] Tuite | The prion hypothesis : from biological anomaly to regulatory mechanism[END_REF], and mathematical modeling on the other hand [START_REF] Sindi | Mathematical Modeling of Prion Disease[END_REF][START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF][START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF][START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF]. One feature of prion propagation that remains poorly understood theoretically is the ability of prions to develop a variety of infectious agents with different phenotypical properties, also called strains [START_REF] Tanaka | The physical basis of how prion conformations determine strain phenotypes[END_REF][START_REF] Derdowski | A size threshold limits prion transmission and establishes phenotypic diversity[END_REF][START_REF] Collinge | A General Model of Prion Strains and Their Pathogenicity[END_REF][START_REF] Colby | De novo generation of prion strains[END_REF][START_REF] Hunter | Learning ability of mice infected with different strains of scrapie[END_REF]. These differences are assumed to be associated with different conformations of the protein, which are then associated with strain-specific biochemical properties [START_REF] Tanaka | The physical basis of how prion conformations determine strain phenotypes[END_REF][START_REF] Weissman | Emerging Principles of Conformation-Based Prion Inheritance[END_REF][START_REF] Lindquist | Structural insights into a yeast prion illuminate nucleation and strain diversity[END_REF]. Mounting evidence suggests that strains can coexist and compete, but it is unclear what drives their interactions [START_REF] Sharma | De novo [PSI + ] prion formation involves multiple pathways to form infectious oligomers[END_REF][START_REF] Bradley | Interactions among prions and prion "strains" in yeast[END_REF][START_REF] Dur | Divergent prion strain evolution driven by PrP C expression level in transgenic mice. Nature Communications Li[END_REF]. We note that the phenomenon of strains applies not only to prion diseases, but also to other neurodegenerative disorders caused by prion-like mechanisms such as Alzheimer's or Parkinson's disease [START_REF] Watts | Serial propagation of distinct strains of Aβ prions from Alzheimer's disease patients[END_REF][START_REF] Cohen | Rapidly progressive Alzheimer's disease features distinct structures of amyloid-beta[END_REF].

Previous mathematical studies of spontaneous prion formation have led to widespread acceptance of the nucleated polymerization model [START_REF] Lansbury | The chemistry of scrapie infection: implications of the ice 9 metaphor[END_REF], which has been formulated under the assumption of either discrete or continuous aggregate sizes [START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF][START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF][START_REF] Doumic | Scaling limit of a discrete prion dynamics model[END_REF][START_REF] Davis | A study in nucleated polymerization models of protein aggregation[END_REF]. This model indeed reproduces qualitative and quantitative aspects of prion infection [START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF][START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF], all the while being analytically tractable [START_REF] Engler | Analysis of a model for the dynamics of prions II[END_REF]. However, as we show here, this model does not explain the interactions between prion strains observed in vivo and in vitro. In particular, various studies demonstrate that strains may coexist in both yeast [START_REF] Strbuncelj | The impact of genetic and epigenetic interactions on propagation of the PSI+ prion[END_REF]) and mammals [START_REF] Morales | Coexistence of multiple PrP Sc types in individuals with Creutzfeldt-Jakob disease[END_REF] but this is not supported by the nucleated polymerization model. While many mathematical studies have generalized the classical model to consider strain specific phenomena, such as [START_REF] Calvez | Prion dynamics with size dependency-strain phenomena[END_REF] who modeled aggregate size-dependent rates, there has been remarkably little mathematical work considering interactions between prion strains. A recent generalization of the nucleated polymerization model showed strain coexistence was possible in yeast when fragmentation was mediated by the molecular chaperone Hsp104 (Davis and Sindi 2016b). However, because chaperones are not known to be part of the in vivo dynamics of mammalian prions, we wanted to identify chaperone independent mechanisms that could explain coexistence of prion strains in both yeast and mammalian systems.

In this work, we generalize the nucleated polymerization model by including a biologically motivated species we term subunits. More specifically, we add a conformational change step prior to polymerization, called templating and driven by subunits [START_REF] Igel-Egalon | Reversible unfolding of infectious prion assemblies reveals the existence of an oligomeric elementary brick[END_REF][START_REF] Igel-Egalon | Prion Strains and Transmmission Barrier Phenomena[END_REF]. Different studies have suggested similar approaches before (Alvarez-Martinez et al 2011; [START_REF] Hingant | A Micellar On-Pathway Intermediate Step Explains the Kinetics of Prion Amyloid Formation[END_REF], but not in the specific case of strain interaction and coexistence. We note that our model exhibits qualitatively distinct behavior from the classical nucleated polymerization model. Most importantly, it enables the coexistence of different strains (multiple strains present stably together), but it also allows for simultaneous co-stability of different equilibria (dependency on initial conditions). These two features will be our main focus throughout the study. Note that our model does not depict the spontaneous formation of prion aggregates because this process happens on a distinct time-scale. The appearance of the first stable prion nucleus is typically modeled as the first-arrival time of a stochastic process, as in D 'Orsogna et al (2012); [START_REF] Yvinec | First passage times in homogeneous nucleation and self-assembly[END_REF][START_REF] Yvinec | First passage times in homogeneous nucleation: dependence on the total number of particles[END_REF]; Davis and Sindi (2016a). In this work, our analytical and numerical results concern the propagation of prions from an initially established pool of aggregates.

In Section 2 we review the nucleated polymerization model and discuss its limitations with respect to multiple strain behavior. Section 3 introduces our model, with primary analytical results on its behavior in Section 3.1 and Section 3.2, a bifurcation analysis in Section 3.3 and a detailed analytical study of the two-strain case in Section 3.4. We end this section with numerical validations in Section 3.5. Finally, we discuss several implications of our model in the context of open questions in prion biology in Section 4.1, we identify the limitations and upcoming challenges in Section 4.2 and we present some justification for our choice over other possible modeling ideas in Section 4.3.

2 Background: The nucleated polymerization model can only predict single strain dominance

Nucleated polymerization model

We will first present the nucleated polymerization model describing prion dynamics when only a single strain is present. In this model, all protein monomers are either in the normally folded state or in the misfolded (prion) state. In the classical nucleation theory, the misfolded form of the protein is only stable when it belongs to an aggregate of size larger than a critical nucleus x 0 [START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF]. The dynamics of the aggregates and normally folded proteins evolve through several biochemical processes as depicted in Figure 1. Aggregates may increase in size when proteins in the normally folded state are converted to the misfolded form and incorporated into aggregates, this is called polymerization [START_REF] Bolton | Identification of a protein that pruifies with the scrapie prion[END_REF][START_REF] Collinge | Prion diseases of humans and animals: their causes and molecular basis[END_REF]. Aggregates may increase in number by fragmenting. In the original formulation, a prion aggregate is assumed to be a linear array of misfolded protein monomers. As such, polymerization can only occur as the ends of an aggregate interact with normal monomers at rate τ. Similarly, fragmentation is assumed to occur at rate β between any two adjacent misfolded monomers in an aggregate. In accordance with the nucleation assumption, if a fragmentation event creates an aggregate with size smaller than the nucleus size x 0 , it is immediately disassembled and the constituent proteins return to the normally folded state. Finally, normal monomers are produced with a source rate λ and degraded with rate γ and aggregates are degraded with rate µ.

In the mathematical literature, both continuous and discrete size aggregate distributions have been considered [START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF][START_REF] Davis | A study in nucleated polymerization models of protein aggregation[END_REF][START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF][START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF]. While the results and approaches taken are similar, in what follows we use the continuous (in aggregate size) formulation studied in [START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF] and [START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF] because our interest is mainly in qualitative study of the dynamics. In this case, the state space dynamics are governed by ordinary differential equation for V , the concentration of normally folded monomers and a partial differential equation, a transport equation, governing the aggregate size distribution u(x,t) as shown below:

         dV dt (t) = λ -γV (t) -τV (t) ∞ x=x 0 u(x,t)dx + 2 x 0 x=0 x ∞ y=x 0 β u(y,t)dydx, ∂ u ∂t (x,t) + τV (t) ∂ u ∂ x (x,t) = -µu(x,t) -β xu(x,t) + 2 ∞ y=x β u(y,t)dy. (1)
(2

)
This system is to be completed with boundary conditions u(x 0 ,t) = 0 and lim x→∞ u(x,t) = 0 for t ≥ 0, as well as initial conditions V (t = 0) = V 0 and u(x, 0) = u 0 (x) for x ∈ [x 0 , ∞[. In many cases, we are interested in the qualitative behavior of the system, such as whether misfolded protein will exist asymptotically, rather than the time-evolving aggregate density. To facilitate this analysis, we convert the equations ( 1)-(2) to a system of ordinary differential equations using moment closure. More specifically, for t ≥ 0 we define U(t) = ∞ x=x 0 u(x,t)dx to be the number of aggregates, the zeroth moment of the aggregate density, and P(t) = ∞ x=x 0 xu(x,t)dx to be the total mass of aggregated protein, the first moment of the aggregate density. We obtain the following system of ordinary differential equations involving the moments U(t), P(t) and the population of normal monomers V (t):

               dV dt =λ -γV -τVU + β x 2 0 U, dU dt =β P -µU -2β x 0 U, dP dt =τVU -µP -β x 2 0 U (3)
where t ≥ 0 and initial conditions V (0) ≥ 0,U(0) ≥ 0 and P(0) ≥ x 0 U(0).

As has been discussed in prior studies [START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF][START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF] this system asymptotically approaches one of two possible steady-states. Either the Disease Free Equilibrium (DFE) corresponding to the presence of only the normally folded protein, (V,U, P) = (λ /γ, 0, 0) or the Endemic Steady-State (ESS) where all three quantities are non-negative indicating the asymptotic stability of prion aggregates:

(V,U, P) = (µ + β x 0 ) 2 β τ , λ β τ -γ(µ + β x 0 ) 2 µτ(µ + 2β x 0 ) , λ β τ -γ(µ + β x 0 ) 2 β µτ . ( 4
)
The stability of the ESS depends on the basic reproductive number of the prion strain,

R 0 = λ β τ γ (µ + β x 0 ) 2 .
(5)

When R 0 > 1, as long as some aggregates are initially present (i.e., U(0) > 0), the system will approach the ESS. These global stability results are proved in [START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF]. This basic reproductive number plays the role commonly observed in epidemiological models [START_REF] Brauer | Mathematical Models in Population Biology and Epidemiology[END_REF].

Multi-strain nucleated polymerization model

As mentioned in the introduction, we are interested in the behavior of multiple prion strains and so we next consider a natural generalization of this model to multiple strains. Similarly to [START_REF] Tanaka | The physical basis of how prion conformations determine strain phenotypes[END_REF], we assume the same biochemical processes apply to all strains and that aggregates corresponding to distinct prion strains do not directly interact. Under these assumptions, strains each have their own aggregates and the aggregates dynamics are governed by strain-specific biochemical rates. That is, in a system with N different strains (N > 1) we have strain-specific parameters x i , τ i , β i and µ i for i ∈ {1, . . . N}. The only difference from the nucleated polymerization model is that strains compete for the same resource, the normally folded monomers. As such, we obtain a similar moment closure for each strain.

The system consists of a single equation for the population of normal monomers, V (t), and one for each moment of each strain, U i (t), P i (t)

                 dV dt =λ -γV + N ∑ i=1 (-τ i VU i + β i x 2 i U i ), dU i dt =β i P i -µ i U i -2β i x i U i , i = 1 . . . N, dP i dt =τ i VU i -µ i P i -β i x 2 i U i , i = 1 . . . N, (6)
where N ≥ 1, t ≥ 0 and initial conditions where for each strain U i (0) ≥ 0, P i (0) ≥ x i U i (0).

Let us now consider the different steady-states possible in our multi-strain nucleated polymerization model given by System (6). Note that different steady-states are possible depending on the parameters. First, the Disease-Free Equilibrium, when all prion strains become extinct, always exists and is given by (V,U 1 , P 1 , . . . ,U N , P N ) = (λ /γ, 0, 0, . . . , 0, 0).

Then for each strain i, i ∈ {1, . . . N}, there is a Strain-Specific Endemic Steady-State corresponding to the asymptotic presence of strain i only. It is given by

(V,U 1 , P 1 , . . . ,U i , P i , . . . ,U N , P N ) = (µ i + β i x i ) 2 β i τ i , 0, . . . , 0, λ β i τ i -γ(µ i + β i x i ) 2 µ i τ i (µ i + 2β i x i ) , λ β i τ i -γ(µ i + β i x i ) 2 µ i β i τ i , 0, . . . , 0 .
We note that the nonzero quantities in the steady-state correspond exactly to the nucleated polymerization endemic steady-state for the strain in isolation (Equation ( 4)). The strainspecific endemic steady-state is biologically feasible (i.e., all concentrations are positive)

Fig. 1: Nucleated polymerization model with a single prion strain, illustrated with a nucleus size of x 0 = 3. The normal form monomers (empty circles) appear with a production rate λ and are degraded with a rate γ. They are recruited by prion aggregates (crossed circles) with rate τ. Each link in an aggregate may break with probability β , and an aggregate below the nucleus size x 0 is disassembled into normal monomers. The aggregates are degraded with rate µ.

only when λ γ β i τ i > (µ i + β i x i ) 2 .
As we did for the single strain case, and has been done in prior studies [START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF][START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF], we define a strain-specific basic reproductive number

R i 0 = λ β i τ i γ(µ i + β i x i ) 2 , (7)
and note that biological feasibility of the strain-specific endemic steady-state corresponds to exactly R i 0 > 1. Next, we consider steady-states where multiple strains exist together. We define such steady-states to be Coexistence Steady-States. We note that for each strain at steady-state, the following strain-specific relation on the normally folded protein monomer density must be satisfied:

V = (µ i + β i x i ) 2 β i τ i = λ γ 1 R i 0 . (8)
As such, only strains with the same basic reproductive number can exist together at steady state. This means that as long as strains have different steady-state normal monomer concentrations, they cannot coexist at steady-state in our multi-strain model. Because the reproductive number of a strain R i 0 is a function of the strain-specific biological parameters, it is highly unlikely that these R i 0 coincide exactly. Moreover, even if by chance two strains had an identical R i 0 , coexistence would not be robust as the slightest perturbation of one parameter would completely remove the possibility of a coexistence equilibrium.

In addition, much like the original nucleated polymerization model, only one steadystate at a time will be asymptotically stable. The behavior of this multi-strain model ( 6) is driven by the following theorem.

Theorem 1 Assume all parameters λ , γ and τ i , β i , x i , µ i for i = 1 . . . N are positive. Then the system (6) admits a unique positive solution for each initial condition taken in

X = {(V,U 1 , P 1 , . . . ,U N , P N ) ∈ R 2N+1 : V,U 1 , P 1 -x 1 U 1 , . . . ,U N , P N -x N U N ≥ 0}.
Further suppose that the quantities R i 0 = λ β i τ i γ(µ i +β i x i ) 2 are distinct and define

R 0 = max i=1...N R i 0 .
If R 0 ≤ 1, the disease-free equilibrium (λ /γ, 0, 0, . . . , 0, 0) is globally asymptotically stable on X. If R 0 > 1, suppose Strain 1 verifies this maximum (renumbering the strains if necessary), there exists an endemic equilibrium for each strain with R i 0 > 1 (including Strain 1). The equilibrium involving only Strain 1 is given by

(µ 1 + β 1 x 1 ) 2 β 1 τ 1 , λ β 1 τ 1 -γ(µ 1 + β 1 x 1 ) 2 µ 1 τ 1 (µ 1 + 2β 1 x 1 ) , λ β 1 τ 1 -γ(µ 1 + β 1 x 1 ) 2 µ 1 β 1 τ 1 , 0, 0, . . . , 0, 0 ,
and is globally asymptotically stable on {(V,U 1 , P 1 , . . . ,U N , P N ) ∈ X : U 1 > 0 and P 1 > 0}.

Note that if U 1 (0) = P 1 (0) = 0, the outcome will be given by the same theorem without considering Strain 1.

Proof The existence and uniqueness of solutions is proved in the same fashion as in the single strain case, with direct adaptation of the proof provided in [START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF]. The global results rely on Lyapunov functions, and the proof is presented in Appendix A.

Biological interpretation and limitations

In the previous section, we demonstrated that the natural generalization of the nucleated polymerization model to multiple strains allows for at most one strain to persist, and the steady-state associated with this strain will be globally asymptotically stable. However, this is not consistent with biological studies suggesting a diversity of behaviors is possible with prion strains. Before we develop our template-assistance model, we first describe biological behavior observed both in yeast and mammalian prion systems that the previously introduced model does not support.

Coexistence of different strains cannot be explained. Numerous experiments demonstrate that distinct prion strains may be present stably at the same time in a host, even though they are associated with dramatically different disease properties [START_REF] Morales | Coexistence of multiple PrP Sc types in individuals with Creutzfeldt-Jakob disease[END_REF][START_REF] Dur | Divergent prion strain evolution driven by PrP C expression level in transgenic mice. Nature Communications Li[END_REF][START_REF] Langenfeld | Incongruity between Prion Conversion and Incubation Period following Coinfection[END_REF]. As emphasized in Section 2.2, coexistence equilibria for the multi-strain nucleated polymerization model only exist when strains have exactly the same values of R i 0 . We note that this is biologically unlikely because it describes a set of measure zero in parameter space. One could argue that, even if the coexistence steadystate is impossible, multiple strains could be present together for a long time before one takes over. This indeed happens when the R i 0 values are very close to each other. However, judging from phenotype differences between strains observed in vivo in mammals as well as yeast [START_REF] Dur | Divergent prion strain evolution driven by PrP C expression level in transgenic mice. Nature Communications Li[END_REF][START_REF] Strbuncelj | The impact of genetic and epigenetic interactions on propagation of the PSI+ prion[END_REF], it seems likely that these coexisting strains possess very different propagation properties. Under these conditions, the multi-strain nucleated polymerization model predicts the fast takeover by one strain. This is the first main flaw in the classical model, and its generalization with multiple strains.

Co-stability of different steady-states cannot be explained. In some experimental studies for mammalian prions, the observable outcome appears to vary based only on changes in the inoculum, i.e. the initial condition. For example see Le [START_REF] Dur | Divergent prion strain evolution driven by PrP C expression level in transgenic mice. Nature Communications Li[END_REF], where diluting the inoculum changes the strain which dominates, and repeating the same experiment (with possibly slight variations in the inoculum) yields different possible outcomes. In [START_REF] Langenfeld | Incongruity between Prion Conversion and Incubation Period following Coinfection[END_REF], superinfection leads to different outcomes depending on the delay between the first inoculation (with one strain) and the second (with another strain). Theorem 1 spotlights a weakness of the model regarding these observations, namely that its behavior is global. Regardless of the case, there is only one equilibrium at a time that is globally asymptotically stable. This means that with this model, the outcome does not depend at all on the initial conditions. In particular, in the endemic case, one strain always takes over the other ones.

Although the nucleated polymerization model is clearly insufficient to investigate multiple strain phenomena, it is still coherent with experimental data in the single strain case, and it is analytically well understood [START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF][START_REF] Engler | Analysis of a model for the dynamics of prions II[END_REF][START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF]. Based on these considerations, we will use the nucleated polymerization model as a "building block" towards a model to study multiple strain phenomena, as we explain in Section 3.

The template assistance model

Through mathematical investigation, we identified different generalizations of the nucleated polymerization model that could allow for multiple strain coexistence and co-stability (see Section 1 for precise definition). The simplest idea that appeared both encouraging and biologically consistent (see Section 4.3 for a parallel with ecological populations and a discussion that supports our modeling choice) was to add a step prior to polymerization, and to consider an additional biochemical species. This species represents subunits (potentially small oligomers) that are responsible for the templating activity (monomer conversion). The aggregation of subunits follows the nucleated polymerization dynamics and Figure 2 illustrates the principles of such a model.

We note that our concept of subunits has a biological basis. Indeed recent studies suggest that for mammalian prions, aggregates are in a kinetic equilibrium with subunits [START_REF] Igel-Egalon | Reversible unfolding of infectious prion assemblies reveals the existence of an oligomeric elementary brick[END_REF][START_REF] Igel-Egalon | Prion Strains and Transmmission Barrier Phenomena[END_REF]. It is still not clear how these subunits are involved in templating in vivo, but we here include them as the main templating agent. In yeast, there is also evidence of in vitro prion formation following different pathways [START_REF] Sharma | De novo [PSI + ] prion formation involves multiple pathways to form infectious oligomers[END_REF]. This modeling approach has been considered in previous studies [START_REF] Hingant | A Micellar On-Pathway Intermediate Step Explains the Kinetics of Prion Amyloid Formation[END_REF][START_REF] Yvinec | Probabilistic modeling in cellular and molecular biology[END_REF], although never with a continuous and deterministic formulation. Additionally, this modeling approach has never been used to investigate interactions between multiple prion strains. We will now formally introduce this generalization and investigate the behavior of our model with respect to multiple strain coexistence and interaction.

Model formulation for a single strain and primary results

We first introduce our template assistance model for a single strain. In our formulation we preserve the notation introduced for the nucleated polymerization model in Section 2, but we now introduce a new species, namely the subunit population S(t) at time t ≥ 0. As shown in Figure 2, the initial conformational change is reversible, and subunits return to the normal form at rate ω. They are degraded with a rate δ , most likely δ < γ representing the fact that the subunits are more resistant than normal monomers. As illustrated in Figure 2, the subunits are formed by direct interaction between subunits and monomers with speed ρV f (S). The function f models the efficiency of subunits to convert monomers. Our choice is

f (S) = S S K + S ,
which depicts a Michaelis-Menten process for the conformational change reaction (a mathematical justification for this choice will come later). Essentially, this reaction is not instantaneous, and the population of subunits S has to be above a certain threshold before reaching Fig. 2: Template assistance mechanism illustrated for a single prion strain. The normal monomers (empty circles) are converted by the subunits (squares) by direct interaction at rate ρ. This first conformational change is reversed spontaneously at rate ω, and the subunits are also degraded at rate δ . The aggregates are formed of subunits, and the dynamics driving the interaction between aggregates and subunits are those of the nucleated polymerization model (similarly to Figure 1, with subunits playing the role of monomers). The aggregates are degraded at rate µ.

full efficiency (K representing the density at which half the maximal efficiency is reached).

One could also say that the fraction of actively involved subunits is given by a Hill function (here of order 1 for simplicity), which is a commonly used class of functions to model enzyme-mediated biochemical kinetics (see for instance [START_REF] Segel | A Primer on Mathematical Models in Biology[END_REF], Chapter 8). A linear templating rate would prevent co-stability as shown in Section 3.2, and the non-linearity introduced appears a sufficient condition for co-stability although we can only show it numerically in Sections 3.3 and 3.5. Proving it analytically remains an open problem that we will investigate in future work. In addition, there is precedence for considering Michaelis-Menten dynamics in prion kinetics, Milto et al ( 2013) considered fibril elongation occurring as an enzyme-mediated reaction. This is not our choice here, but we suggest the templating by subunits is a cooperative reaction. The elongation rate could indeed be modified to obtain the same behavior as we obtain here, but it would require inconsistent hypotheses with the biological context we focus on, see Section 4.3 for a detailed discussion. As long as a quantitative study is not experimentally possible, to the best of our knowledge (see Section 4.2), we consider the simplest option for the elongation rate.

The full aggregate dynamics are now defined by two ordinary differential equations (for V (t) and S(t)) and a single partial differential equation for u(x,t), which will be identical to Equation (2) but with S(t) replacing V (t) in the transport term. Under these assumptions, we obtain a similar moment closure to the equations of the nucleated polymerization model

                       dV dt =λ -γV -ρV f (S) + ωS, dS dt =ρV f (S) -ωS -δ S -τSU + β x 2 0 U, dU dt = -µU + β P -2β x 0 U, dP dt =τSU -µP -β x 2 0 U. (9)
Similarly to the nucleated polymerization model, the System (9) has a unique positive solution. In particular, we prove the following Lemma.

Lemma 1 When the parameters λ , γ, β , τ, µ, x 0 , ρ, ω, δ , K are all positive, the system (9) admits an unique positive solution for each initial condition taken in

X = {(V, S,U, P) ∈ R 4 : V, S,U, P -x 0 U ≥ 0}.
Proof It is simple to verify that X is positively invariant. Furthermore, for t > 0, V (t) + S(t) + P(t) is bounded between 0 and λ ε + (V (0) + S(0) + P( 0))e -εt , with ε = min{γ, δ , µ}. The proof is concluded as in [START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF].

Steady-state analysis

We will next investigate the existence and, when possible, the linear stability of steady-states of our system. Because it will be useful in our analysis, we state here the Jacobian matrix of System (9) for clarity

J(V, S,U, P) =   -γ-ρ f (S) ω-ρV f (S) 0 0 ρ f (S) ρV f (S)-ω-δ -τU -τS+β x 2 0 0 0 0 -µ-2β x 0 β 0 τU τS-β x 2 0 -µ   .
Proposition 1 The disease-free equilibrium exists for any choice of positive parameters, and is given by

(V, S,U, P) = λ γ , 0, 0, 0 . It is locally stable if and only if λ γ f (0) < ω + δ ρ .
In particular when f (0) = 0, as with the choice f (S) = S 2 /(K + S) introduced in Section 3.1, the disease-free equilibrium is locally stable for any choice of positive parameters.

Proof A straightforward calculation yields the eigenvalues of the Jacobian matrix at this points {-γ, -µβ x 0 , -µβ x 0 , λ γ f (0)ωδ }. For positive parameters, the local stability depends only on the last one, and it is negative when the condition expressed above is verified.

The local stability of the disease-free equilibrium is the new feature that will enable costability of different equilibria. Any other function verifying f (0) = 0 would also yield the same property, but our choice was to introduce as few parameters as possible. On the contrary, we can show that a linear templating rate would not provide co-stability, as explained at the end of this section. A second type of steady-state is possible, when subunits are present but no aggregates.

Proposition 2 There exists at most two subunits-only equilibria, and they are given by

(V, S,U, P) = λ γ - δ γ S +/-, S +/-, 0, 0 with S +/-= 1 2 λ δ - ω + δ ρ γ δ ± 1 2 λ δ - ω + δ ρ γ δ 2 - ω + δ ρ γ δ K.
These equilibria are feasible (positive values) when the following condition is met

λ γ > ω + δ ρ + 2 δ γ ω + δ ρ γ δ K.
The equilibrium associated with S -is always locally unstable. The equilibrium associated with the higher value is locally stable if and only if

S + < (µ + β x 0 ) 2 β τ .
Proof If we impose U = 0 and P = 0, the remaining two equations on V and S lead to the equation

S 2 + ω+δ ρ γ δ -λ δ S + ω+δ ρ γ δ K = 0.
This gives directly the expression of the two possible values for S, and the condition for them to be real and positive. Their stability is once again analyzed through the Jacobian matrix. For both of the points S +/-, a pair of eigenvalues is given by -µβ x 0 + τβ S +/-, -µβ x 0τβ S +/-. These eigenvalues are both negative when S +/-< (µ+β x 0 ) 2 β τ

. The other pair of eigenvalues are the roots of a second degree polynomial (not shown for simplicity), and after simplification we show that they are negative when S +/-> ω+δ ρ γ δ K. Recalling the values of S + and S -, we see that this second condition is never met by S -and always met by S + . Hence the condensed result in the proposition.

This result shows that the two subunits-only equilibria emerge through a saddle-node bifurcation, and that only the one associated with a higher S value might be stable, with condition related to the last equilibrium. See Section 3.3 for a numerical illustration of this bifurcation.

Proposition 3 The endemic steady-state (V * , S * ,U * , P * ) is defined by the following relations

S * = (µ + β x 0 ) 2 β τ , V * = λ + ωS * γ + ρ f (S * ) , U * = 1 τS * -β x 2 0 (ρV * f (S * ) -ωS * -δ S * ), P * = µ + 2β x 0 β U * .
This equilibrium is feasible (positive values) when

λ > δ S * + ω + δ ρ γS * f (S * ) .
Proof The feasibility mainly relies on U * > 0, which with the value of V * reduces to the condition expressed above (noticing that τS *β x 2 0 = µ(µ + 2β x 0 )/β > 0). Recalling the conditions for the feasibility of the two subunits-only equilibria, we note that the endemic steady-state exists when the associated subunit density S * is comprised between the two subunits-only values S +/-. Indeed the condition expressed in Proposition 3 implies (but is not equivalent to) f (S * ) > ω+δ ρ γS * λ -δ S * , and the two subunit-only equilibrium are the solutions that verify equality in this relation. We investigate the stability of this equilibrium numerically through a bifurcation analysis in the next section.

Remark 1 Although we cannot show analytically that the condition f (0) = 0 yields costability of different equilibria, we can show that when f is linear there is no co-stability possible. Indeed, if f (S) = S, which is equivalent to K = 0 in all the developments above, there is only one subunits-only equilibrium associated with the subunit density Ŝ = λ δ -ω+δ ρ γ δ . The stability condition for the disease-free equilibrium is reduced to λ γ < ω+δ ρ or equivalently Ŝ < 0, which shows that these two equilibria cannot be both feasible and stable at the same time. The subunits-only equilibrium is in turn stable as long as Ŝ < (µ+β x 0 ) 2 β τ = S * (with the notation of Proposition 3. However, the feasibility condition of the endemic equilibrium is now reduced to V * > ω+δ ρ (V * defined in Proposition 3) or equivalently λ δ -ω+δ ρ γ δ = Ŝ > S * after simplification. This proves that the subunits-only equilibrium and the endemic equilibrium cannot be both feasible and stable at the same time. Essentially, when f is linear there are only three possible equilibria and they appear through a series of transcritical bifurcations, thus there can be no stability.

Numerical bifurcation analysis

As shown in Section 3.2, the subunits-only steady-states appear through a saddle node bifurcation (Proposition 2). The endemic steady-state appears when the associated subunit population S * crosses one of the two branches S + or S -(Proposition 3). If it appears through the higher branch, the high subunit only steady-state S + is stable at first, and then becomes unstable. If the endemic steady-state instead appears through the lower branch, none of the subunit only steady-states will be stable (see conditions in Section 3.2). These two scenarios are illustrated in Figure 3. As the monomer source rate λ increases, we can see the saddlenode bifurcation followed by the emergence of the endemic steady-state. This is a numerical illustration of the results proved in the previous section.

When the endemic steady-state exists, it can be locally stable or unstable. Numerical exploration shows that it undergoes a Hopf bifurcation. This is illustrated with a numerical two-parameter bifurcation analysis in Figure 4. We show dependence on two key parameters: the monomer synthesis rate λ and the polymerization rate τ. The bifurcation analysis is based on the properties of the Jacobian matrix at the endemic steady-state, and we characterize the presence of a Hopf bifurcation as when one eigenvalue crosses the imaginary axis with a non-zero imaginary part. Figure 4 illustrates regions of qualitatively distinct dynamics and the boundaries between these regions correspond to different bifurcation events. The endemic steady-state is only stable in Region 4, and the Hopf bifurcation occurs at the transition between Regions 4 and 5. In Region 5, only the disease-free equilibrium is locally stable. The solutions can either be attracted by this equlibrium or undergo stable oscillations. The exact nature of these oscillations is not known analytically, but the most important qualitative result here is that by crossing this Hopf bifurcation, and aggregate population can possibly be destabilized. The fact that oscillations have been observed experimentally (H. Rézaei, unpublished data) is another reason for interest in that behavior. Overall, the complexities of bifurcations in our model demonstrate the endemic steady-state is not always stable, and this is novel behavior to the nucleated polymerization model. 1 (Strain 1). The four different equilibria are depicted (disease-free equilibrium, endemic equilibrium, higher subunit only equilibrium (S + ), lower subunit only equilibrium (S -)). A thin line indicates an unstable equilibrium, whereas a bold line indicates a locally stable equilibrium Fig. 4: Two-parameter bifurcation diagram obtained for the monomer source rate λ and the polymerization rate τ, when the other parameters are described in Table 1 (Strain 1). The vertical line delimits the apparition of the subunit only equilibria, left of that line (Region 1) only the disease-free equilibrium exists and it is then globally stable (not proved). The branches delimit the existence of the endemic equilibrium. In Region 2, the subunit only equilibria exist but neither is stable. In Region 3, the subunit only equilibria exist and the higher one (S + ) is locally stable. In-between the branches, the endemic steady-state exists, but it is only locally stable in Region 4. In Region 5, the endemic steady-state is unstable. The boundary between Regions 4 and 5 corresponds to a Hopf bifurcation

Two-strain case

As we have emphasized throughout, our primary motivation in developing the template assistance model was to explore the system when multiple prion strains are present. We next generalize our model to two strains, as illustrated in Figure 5. The corresponding generalization to System (9) is the following:

                       dV dt =λ -γV -ρ 1 V f 1 (S 1 ) + ω 1 S 1 -ρ 2 V f 2 (S 2 ) + ω 2 S 2 , dS i dt =ρ i V f i (S i ) -ω i S i -δ i S i -τ i S i U i + β i x 2 i U i , i = 1, 2, dU i dt =β i P i -µ i U i -2β i x i U i , i = 1, 2, dP i dt =τ i S i U i -µ i P i -β i x 2 i U i , i = 1, 2. ( 10
)
The interaction functions, f i , are the same as in System ( 9), but now with a strain specific value of K i : f i (S) = S 2 /(K i + S).

As we will next detail, the behavior of this System (10) is dramatically different from that of the multi-strain nucleated polymerization model (System (6) presented in Section 2). Our new system allows for coexistence of prion strains; that is, each strain may exist in the same conditions as single strains did Section 3.2. Furthermore, two strains may coexist under different configurations: both as subunits only, or one endemic and the other as subunits only (we will refer to this case as semi-endemic) or both endemic. We note that subunit only coexistence is impossible for general values of the parameters, with the same argument that was used to rule out coexistence equilibria in the multi-strain nucleated polymerization model in Section 2.2. Finally, numerical results suggest that semi-endemic coexistence is very unlikely to occur.

The coexistence steady-state is given by

S * i = (µ i + β i x i ) 2 β i τ i for i =1,2, V * = λ + ω 1 S * 1 + ω 2 S * 2 γ + ρ 1 f 1 (S * 1 ) + ρ 2 f 2 (S * 2 )
,

U * i = 1 τ i S * i -β i x 2 i (ρ i V * f i (S * i ) -ω i S * i -δ i S * i ) for i = 1, 2, P * i = µ i + 2β i x i β i U * i for i =1,2. Notice τ i S * i = µ i (µ i + 2β i x i )/β i > 0.
Similarly to the single-strain case, this equilibrium exists when all values of the variables are positive yielding the conditions

f i (S * i ) > ω i + δ i ρ i S * i V * for i = 1, 2.
While we were unable to simplify these conditions, they are easily verified for a given combination of parameters. In particular, we note that for any choice of kinetic rates for prion strains (β i , τ i , ω i , ρ i ), the coexistence steady-state can always be made to exist by increasing the monomer source rate λ enough to satisfy the above conditions (because increasing λ increases V * without changing the values of S * 1 and S * 2 ). The second improvement of our System (10) over the multi-strain nucleated polymerization is that it allows for co-stability. As is the case with a single strain, the disease-free equilibrium is always locally stable (easily proved with the Jacobian matrix). Again as in the single-strain case, the co-stability of different types of equilibria (single-strain steady-states and coexistence steady-states) could not be proved analytically but can be observed through numerical exploration. As such, the outcome depends on the initial conditions as we will further investigate in the next section.

Numerical results

In this section we provide a detailed numerical study of the behavior of the template assistance model we developed in the previous section. We first take care to verify that the behavior in the single strain case remains similar to the original nucleated polymerization model. We next study the dependency on initial condition when multiple steady-states are locally stable by numerically investigating the basins of attraction. Notice that in the different time evolution figures we produce, the evolution of P, the first moment of the aggregates distribution i.e. the total mass or number of aggregated monomers, is not displayed because it is very similar to the evolution of U, the zero-th moment of the distribution i.e. the total number of aggregates. However, P is on a different scale of values than the other variables so for the sake clarity we do not show it.

The dynamics of the nucleated polymerization model are conserved. Because of its broad acceptance, it is crucial that the qualitative behavior of the nucleated polymerization model be maintained after modification. A generic study case is presented in Figure 6, comparing the dynamics of the classical model with our model. One can see the overall dynamics of aggregate formation are qualitatively and quantitatively very similar for long times. Even Fig. 5: Template assistance mechanism illustrated with two different strains (squares and triangles). The strains compete for normal monomers through the interaction with subunits (with rates ρ 1 and ρ 2 ). The dynamics of each strain taken individually are then described as in Figure 2. though it is not shown, the evolution of P(t) also coincides with the one observed in the nucleated polymerization model for long times. The behavior during early times is dramatically different, because a new species is introduced. Our focus here is not on transient dynamics because we lack the data to study them. This could be a potential way to discriminate between different models. This shows that our model brings new possibilities without eliminating previously supported behavior [START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF][START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF][START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF], especially on aggregate formation, size and numbers.

The template assistance model exhibits co-stability. Even though the nucleated polymerization dynamics can be reproduced, our model offers a new variety of behaviors even for the single strain case. Because the disease-free equilibrium remains locally stable, our system exhibits dependency on the initial conditions. Figure 7 represents the basins of attraction of the endemic steady-state and the disease-free steady-state, for the case studied in Figure 6. This helps to visualize the complexity and the diversity of behaviors allowed by our model. A biological interpretation for the shape of the basins plotted in Figure 7 would be that an outbreak of prion propagation is possible only by having a specific initial mix of subunits and aggregates. In particular, starting off with too many aggregates (U(0) high) would prevent the subunits pool from successfully building up. We want to stress out that the important Table 1 (Strain 1). The initial condition is given by V (0) = λ γ = 200,U(0) = 0.1, P(0) = (2x 0 + µ β )U(0), S(0) = 0.2 aspect of Figure 7 is the shape of the basins, the precise values for the borders are irrelevant due to the arbitrary choice in parameters. Furthermore, this shape is subject to change when modifying the parameters, although so far we have no way of predicting these changes analytically. A numerical investigation could be undertaken, but in order to make biologically relevant predictions we require experimental observations or experimentally derived biochemical parameter values (see Section 4.2).

The two-strain case exhibits coexistence and co-stability. To support the analytical results on strain coexistence (Section 3.4) and give numerical evidence of co-stability we show some simulations with two strains. Figure 8 illustrates the different possible outcomes, depending on the initial conditions. The set of parameters used (see Table 1) allows for four Fig. 7: Representation of the basins of attraction of the endemic steady-state (dark color) and the disease-free steady-state (lighter color), according to the initial condition. The horizontal axis represents the initial number of subunits, and the vertical axis represents the initial number of aggregates. The parameters used are those of Strain 1 in Table 1, the initial value for V is set to λ /γ (disease-free value), and the initial value for P is set to (2x 0 + µ/β )U(0) (steady-state mean size) different equilibria to be co-stable (disease-free equilibrium, two single strain equilibria and coexistence equilibrium). Depending on the initial mass of each strain, the outcome can be the takeover of one strain or the stable coexistence of both, or the extinction of both.

The shapes of the different basins of attraction are complex and non-intuitive. Indeed, in some cases increasing the initial amount of one prion strain can lead to its extinction (see the right-hand side of Figure 8b). Although such behavior seems highly complex, we note that coexistence and co-stability are phenomena supported by experimental studies on mammalian prions [START_REF] Langenfeld | Incongruity between Prion Conversion and Incubation Period following Coinfection[END_REF].

Discussion

In this work, we have developed a novel model of prion aggregate dynamics which, in contrast to the traditional model in the field, is capable of supporting recent biological observations in both yeast and mammalian prion systems (Le [START_REF] Dur | Divergent prion strain evolution driven by PrP C expression level in transgenic mice. Nature Communications Li[END_REF][START_REF] Strbuncelj | The impact of genetic and epigenetic interactions on propagation of the PSI+ prion[END_REF]. While a lack of detailed experimental findings have made precise quantitative comparisons impossible, the qualitative analysis of our model provides insight into a number of prion phenomena as we discuss below.

(a) Linear scale (b) Logarithmic scale Fig. 8: Representation of the basins of attraction of the different equilibria in the two-strain case, in linear scale (8a) and logarithmic scale (8b). The parameters used are described in Table 1. Each strain is initialized with its steady state proportions (see Section 2.1) and diluted with a specific rate (between 1 and 1.10 -5 ), and V is initialized to the disease-free value λ /γ. The axes represent the initial mass of each strain, i.e. the initial value of P + S for each strain (Strain 1 on the horizontal axis, Strain 2 on the vertical axis). The colors represent the outcome of the simulation (with colorbar above), one strain only or the two strains stably coexisting. Notice that the disease-free equilibrium is observed only in the bottom left corner

Qualitative insight into open problems in prion biology

Our model suggests a mechanism for prion template-mediated conformational change. The essential idea is to involve a small oligomeric subspecies as the main conformers, as supported by recent experimental results [START_REF] Igel-Egalon | Reversible unfolding of infectious prion assemblies reveals the existence of an oligomeric elementary brick[END_REF][START_REF] Igel-Egalon | Prion Strains and Transmmission Barrier Phenomena[END_REF]. This change in the nucleated polymerization model induces qualitative behavior that could account for a number of so far unexplained features of prion propagation and transmission. For most of these aspects, the idea is that separation between the templating efficiency and the polymerization efficiency, associated with a non-linearity in the templating activity, can cause intricate outcomes. Some strains will be associated with high numbers of subunits but low numbers of aggregates, or vice-versa. Most importantly, strain behavior would be modified depending on the other strains present and the initial configuration of aggregates versus subunits. Overall, judging from analytical and numerical investigation, we hypothesize that our model provides new and useful insight into prion dynamics. In particular, while our model was designed specifically to investigate mammal prion phenomena, similarities with yeast prions are indisputable. The general idea of a more complex polymerization pathway could apply profitably to the case of yeast prions. In particular, different oligomer pathways have been uncovered for yeast prions [START_REF] Sharma | De novo [PSI + ] prion formation involves multiple pathways to form infectious oligomers[END_REF], which suggests that our model could be modified to study yeast prion strains. The aspects of prion dynamics that could be investigated with our model include, but are not limited to, the following elements.

Strain coexistence is made possible by the consideration of a new prion species, as we emphasize in Figures 8 and9 and in Section 3.5. Indeed, under favorable environment configuration and initial conditions, two strains can coexist stably. This a major improvement compared to the nucleated polymerization model, which does not allow for coexisting strains in a robust parameter range, the competition for normal monomers being too strict (see Section 2.2). As was observed for the Creutzfeldt-Jacob Disease, and other prion associated disorders, multiple strains or substrains (Li et al 2010) may be present simultaneously in a host [START_REF] Langenfeld | Incongruity between Prion Conversion and Incubation Period following Coinfection[END_REF]. It remains to be investigated quantitatively if the coexistence exhibited by our model can indeed correlate with the levels of prion coexistence in vivo.

Coexistence of strains was also observed in yeast prions [START_REF] Strbuncelj | The impact of genetic and epigenetic interactions on propagation of the PSI+ prion[END_REF]), and our model shows that including an additional step in the polymerization pathway would give a reasonable explanation for this behavior. Other ideas have been suggested for the case of yeast prion strains, in particular the influence of molecular chaperones (Davis and Sindi 2016b).

Inoculum influence is another new feature of this model that could explain various experimental results of sub-propagation. The basins of attraction shown in Figures 7 and8 can only give a glimpse of how complex our model's behavior really is. The intricateness of the different basins suggests that a slight variation of inoculum (i.e., introduced aggregate composition) and initial state could yield dramatically different outcomes, even more so when multiple strains are present. Changing slightly the initial proportions between the two strains can yield any of the three possible outcomes (one or the other strain taking over, or both strains coexisting). These outcomes in our model are supported by different experimental studies [START_REF] Langenfeld | Incongruity between Prion Conversion and Incubation Period following Coinfection[END_REF][START_REF] Dur | Divergent prion strain evolution driven by PrP C expression level in transgenic mice. Nature Communications Li[END_REF] as we develop in Section 2.2. Overall it seems biologically relevant that the inoculum heavily influences the outcome, which is not the case with the classical nucleated polymerization model (see Section 2.3).

Monomer expression level is a factor that influences strain selection and sub-propagation (Le Dur et al 2017). Our model does exhibit qualitatively different behaviors for different levels of monomer synthesis, as suggested by the bifurcation analysis in the single strain case, see Figure 4. We also observed qualitatively different behavior in the two strain behavior when modifying the monomer source rate λ . In addition, for yeast a stable prion strain could be made unstable by increasing the rate of monomer synthesis. Although this has been attributed to balance between monomers and chaperones [START_REF] Derdowski | A size threshold limits prion transmission and establishes phenotypic diversity[END_REF]Davis and Sindi 2016b), our template assistance model provides another possible hypothesis. We plan to investigate the dependency on the monomer source rate in future studies.

Strain interaction could be investigated with our model. Prior experimental studies reported that coinfection or superinfection by different prion strains can yield various and remarkable results. For instance Shikiya et al ( 2010) and [START_REF] Langenfeld | Incongruity between Prion Conversion and Incubation Period following Coinfection[END_REF] sequentially inoculate different strains in a single host. Depending on the delay between two inoculations and the concentration ratio between the two inocula, the second strain can be totally suppressed or both strains can be made to coexist. Based on numerical studies, we suggest that such behavior could be explained with our model. For instance, disturbing a single strain equilibrium by adding a small inoculum of another strain can yield various results, blocking the new strain or eliminating the existing strain. This behavior is illustrated in Figure 8 (in particular the bottom-right corner). Similar behavior has been observed in yeast prions [START_REF] Bradley | Interactions among prions and prion "strains" in yeast[END_REF], and this suggests once again that the scope of our model could extend to the case of yeast prions.

Decorrelation between aggregation and infectivity rise is a feature reported numerous times in experimental studies [START_REF] Rubenstein | Scrapie-Infected Spleens: Analysis of Infectivity, Scrapie-Associated Fibrils, and Protease-Resistant Proteins[END_REF][START_REF] Mays | Prion Infectivity Plateaus and Conversion to Symptomatic Disease Originate from Falling Precursor Levels and Increased Levels of Oligomeric PrP Sc Species[END_REF]. It has been observed that the infectivity level of a prion-infected host's brain can become significantly high before any detection of aggregated mass. The behavior of the template assistance model could account for such a decorrelation. Indeed with this model, the formation of aggregates is always preceded by a transitional phase with a very high number of subunits and almost no aggregates, as shown in Figure 6. This transitional phase would not be detected by classical PrP Sc monitoring methods, because of the low stability of these subunits. During this phase, the infectivity could already be high because of the high number of subunits. This transient phase could also account for the non-null delay observed even in seeding experiments, as suggested before [START_REF] Alvarez-Martinez | Dynamics of polymerization shed light on the mechanisms that lead to multiple amyloid structures of the prion protein[END_REF][START_REF] Hingant | A Micellar On-Pathway Intermediate Step Explains the Kinetics of Prion Amyloid Formation[END_REF]. This raises some critical questions for future study as we discuss below. Finally, in yeast prions, the notion of propagon could be related to this idea [START_REF] Derdowski | A size threshold limits prion transmission and establishes phenotypic diversity[END_REF]. The minimal infectious agent, the propagon unit, would be a proper mix of aggregates and subunits.

Challenges and limitations

In order to compare the numerical results we observe with in vivo behavior, we must be able to infer realistic values for the model parameters. Few studies have investigated parameter inference before [START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF][START_REF] Derdowski | A size threshold limits prion transmission and establishes phenotypic diversity[END_REF], and then only in the case of a single strain. As such, more detailed experimental studies involving single and multiple strains will provide valuable data for refining parameter values. However, before we can hope to fit parameters, it must be clear how to link biological observables (i.e., prion phenotypes) with quantities in our mathematical model. What biologists observe and quantify are global properties related to the disease (incubation time, deposition pattern, molecular density distribution, migration patterns, see [START_REF] Morales | Coexistence of multiple PrP Sc types in individuals with Creutzfeldt-Jakob disease[END_REF] for a detailed description), but it is not clear how these properties are linked to the biochemical kinetics of aggregate formation. Accumulating studies suggest that phenotypes (disease properties) are not necessarily linked to prion strains (conformational states of the protein) in a one-to-one relation. Some phenotypes could potentially be composed of a mixture of different strains, and in some cases, strains are propagated "silently" (not influencing the phenotype) (Le [START_REF] Dur | Divergent prion strain evolution driven by PrP C expression level in transgenic mice. Nature Communications Li[END_REF]. This makes it even more complicated to determine each strain's characteristics.

Finally, by introducing a novel species (prion subunits) we have continued to complicate the question of what is the infectious agent of prion phenotypes. It is not clear what causes disease and what causes the onset of symptoms. Some insight from our primary study would be that infectivity might be optimal when the inoculum consists of a specific mix of aggregates and subunits. Disturb these proportions and the propagation fails. Our work also suggests that these proportions are strain-specific, and dependent on the environment or on the presence of other strains.

Ecological parallel and other modeling ideas

The coexistence of strains allowed by our model is similar to the phenomenon of "predation-mediated coexistence" in ecology. First, in the basic model (Section 2) the aggregates, predators, compete for the same resource, monomers. As such, they cannot coexist which is coherent with the competitive exclusion principle in ecology. In our model (Section 3), the aggregates are still predators but they do not compete for the same subunits, the equivalent of preys. These subunits do compete for a shared resource, but the predation limits their proliferation and thus their resource uptake. This allows different strains of subunits to coexist (as long as they are in presence of aggregates). This type of coexistence was exhibited in some ecological systems, for instance the Daphnia planktons [START_REF] Declerck | Impact of fish predation on coexisting Daphnia taxa: a partial test of the temporal hybrid superiority hypothesis[END_REF][START_REF] Gliwicz | Predation-Mediated Coexistence of Large-and Small-Bodied Daphnia at Different Food Levels[END_REF].

Judging from this rationale, other ideas could potentially yield the same outcomes, namely coexistence and co-stability. One idea would be to bound the monomer uptake of aggregates by including a non-linear aggregation rate in the equations of the nucleated polymerization model (see Section 2.1). Indeed if the polymerization speed saturated with high numbers of aggregates, we could potentially obtain coexistence. It could for example model the fact that the aggregates are in a constrained volume and their probability to encounter free monomers goes to 0 as their numbers grow. In the case of mammal brains or yeast cells, the aggregate densities at play are very low and the available volume is not limiting. This hypothesis is not relevant in our biological context, and this is why we do not develop it here.

An other idea is related to the biology of yeast prions. In yeast, the propagation of prions, specifically the fragmentation process, has been shown to be dependent on different molecular chaperones, including Hsp104 [START_REF] Jones | Chaperoning prions : the cellular machinery for propagating an infectious protein ?[END_REF][START_REF] Kryndushki | Yeast [PSI + ] Prion Aggregates Are Formed by Small Sup35 Polymers Fragmented by Hsp104[END_REF][START_REF] Winkler | Hsp70 targets Hsp100 chaperones to substrates for protein disaggregation and prion fragmentation[END_REF]. While chaperones have not been identified as part of in vivo mammalian prion propagation, it is possible that they may be involved [START_REF] Romanova | Hsp 104 and Prion Propagation[END_REF]. Davis and Sindi (2016b) modeled a chaperone-limited fragmentation process, and numerically observed the possibility of coexistence of different strains. We intend to explore such possibilities in future studies, in particular for the case of yeast prions.

Conclusion

We first presented the natural generalization of a classical model of prion dynamics in the context of multiple strains. We have shown that it does not reasonably explain the different experimental behaviors and phenomena associated with prion strains. We have then suggested a more complex mechanism of aggregate formation. Our model involves a novel misfolded species termed subunits, and allows for a new variety of behaviors. In particular, we show analytically that coexistence of different strains and co-stability of different equilibria are possible under our model. This suggests interpretations for experimental results where prion strains can interact, coexist or out-compete each other, mainly using the dependency on initial conditions and the presence of bifurcations.

In this work, most of our results are qualitative, but we plan to dedicate future work to quantitative investigation of this model in comparison with experimental results. To do so, it may be necessary to consider different formulations of the model, including discrete aggregate sizes but also stochastic simulations. More experimental results will be needed in order to conduct such a study. In addition, during our mathematical investigation, other ideas came up that could be potential ways of modeling multiple strain interaction. We chose the simplest idea that was consistent in the biological context of focus, namely mammalian prions in vivo. We were driven towards simplicity by the lack of data, but we are confident that collaboration with biologists will allow for quantitative studies and adjustments to the model.

A Proof of Theorem 1

For the existence and uniqueness of solutions we refer to [START_REF] Prüss | Analysis of a model for the dynamics of prions[END_REF], from which the results can easily be extended to the general N strain case.

A.1 Global stability of the DFE

When R 0 = max i=1...N R i 0 ≤ 1, we define the function L given by L(V,U 1 , P 1 . . . ,U N , P N ) = 1 2 V - λ γ 2 + N ∑ i=1 b i µ i β i U i + P i , with b i = 2 (µ i + β i x i ) 2 β i τ i - λ γ - β 2 i x 2 i β i τ i .
This function is a Lyapunov function for the DFE (λ /γ, 0, 0, . . . , 0, 0). Indeed it is positive, because

R i 0 ≤ 1 implies (µ i +β i x i ) 2 β i τ i ≥ λ γ and thus b i ≥ (µ i +β i x i ) 2 -β 2 i x 2 i β i τ i > 0 for i = 1 . . . N.
It evaluates to 0 at the DFE, and its derivative along trajectories is

L = V - λ γ λ -γV - N ∑ i=1 (τ i VU i + β i x 2 i U i ) - N ∑ i=1 b i U i µ 2 i β i + 2µ i x i -τ i V + β i x 2 i , = - 1 γ V - λ γ 2 - N ∑ i=1 τ i U i V 2 -V λ γ + β i x 2 i τ i + b i + x 2 i λ β i γτ i + b i (µ i + β i x i ) 2 τ i β i . The choice of b i yields L = - 1 γ V - λ γ 2 - N ∑ i=1 τ i U i V 2 -2V (µ i + β i x i ) 2 β i τ i + (µ i + β i x i ) 2 τ i β i 2 + µ i (µ i + 2β i x i ) β i τ i (µ i + β i x i ) 2 β i τ i - λ γ , L = - 1 γ V - λ γ 2 - N ∑ i=1 τ i U i V - (µ i + β i x i ) 2 β i τ i 2 + µ i (µ i + 2β i x i ) β i τ i (µ i + β i x i ) 2 β i τ i - λ γ .
It is simple to verify that, under the hypothesis that R i 0 ≤ 1 for each strain, the derivative of L is negative along the trajectories, and evaluates to 0 at the DFE only.

A.2 Global stability of the endemic steady-state

As in the Theorem 1, let us assume all strains have distinct R i 0 values, and that Strain 1 has the maximal value R 0 = R 1 0 = max i=1...N R i 0 . It is convenient for this part to consider a modification of the model, namely choosing to study W i = P ix i U i instead of P i for each strain. The equations of the model now read

dV dt =λ -γV (t) + N ∑ i=1 (-τ i V (t)U i (t) + β i x 2 i U i (t)), dU i dt =β i W i (t) -(µ i + β i x i )U i (t), i = 1 . . . N, dW i dt =τ i V (t)U i (t) -(µ i + β i x i )W i (t), i = 1 . . . N.
In this system, the disease steady-state (involving Strain 1) is given by V * = (µ 1 +β 1 x 1 ) 2 β 1 τ 1 , U * 1 = λ β 1 τ 1 -γ(µ 1 +β 1 x 1 ) 2 µ 1 τ 1 (µ 1 +β 1 x 1 ) ,W * 1 = µ 1 +β 1 x 1 β 1 U * 1 , ∀ j ∈ {2 . . . N},U * j = P * j = 0. When R 0 = R 1 0 = max i=1...N R i 0 > 1, we define the function L given by

L(V,U 1 ,W 1 , . . . ,U N ,W N ) = 1 + β 1 x 2 1 τ 1 V * -β 1 x 2 1 (V -V * -V * log(V /V * )) + µ 1 + β 1 x 1 β 1 (U 1 -U * 1 -U * 1 log(U 1 /U * 1 )) + (W 1 -W * 1 -W * 1 log(W 1 /W * 1 )) + 1 + β 1 x 2 1 τ 1 V * -β 1 x 2 1 N ∑ j=2
(1 -1 η j )( µ j + β j x j β j U j +W j ), with η j = β 2 j x 2 j β j τ j V * + (µ j +β j x j ) 2 β j τ j -V * (µ j +β j x j ) 2 β j τ j -β 2 j x 2 j β j τ j (µ j +β j x j ) 2 β j x j -

β 2 j x 2 j β j τ j -V * 2 for j = 2 . . . N.
This function is a Lyapunov function for the disease-steady state mentioned above. First it is positive, recalling that the classical function φ : ξ → ξξ *ξ * log(ξ /ξ * ) for ξ * > 0 is positive on R * + , with a global minimum of 0 in ξ = ξ * . Notice also that τV * > β 1 x 2 1 . The terms in the sum are also all positive, because for j ∈ {2 . . . N}, η j > 1. Indeed, the numerator of η j can be expressed as (µ j + β j x j ) 2 β j x j -β 2 j x 2 j β j τ j -V * 2 denominator of η j +V * (µ j + β j x j ) 2 β j τ j -V * + β 2 j x 2 j β j τ j (µ j + β j x j ) 2 β j τ j -β 2 j x 2 j β j τ j .

The fact that Strain 1 has the highest R i 0 value means it is associated with the lowest (µ i +β i x i ) 2 β i τ i

. With the formulation above, this proves that for all j = 2 . . . N, η j > 1 and thus 1 -1 η j > 0. It is easy to verify that L evaluates to 0 at the point of the disease steady-state.

To compute the derivative of L along the trajectories, recall that d dt φ (ξ ) = ξ -ξ * ξ dξ dt . Straightforward calculations lead to, after simplifications, L(V,U 1 ,W 1 , . . . ,U N ,W N ) = γ + (γ + τ 1 U 1 )

β 1 x 2 1 τ 1 V * -β 1 x 2 1 (V -V * ) 2 V + τ 1 V * U * 1 3 - V * V - W 1 W * 1 U * 1 U 1 - V V * U 1 U * 1 W * 1 W 1 -1 + β 1 x 2 1 τ 1 V * -β 1 x 2 1 N ∑ j=2 τ j U j V 1 η j V 2 +V η j (µ j + β j x j ) 2
β j τ j -V * -β 2 j x 2 j β j τ j -(µ j + β j x j ) 2 β j τ j + η j β 2 j x 2 j β j τ j V * .

Details of the calculations are not shown but can be made available. The first term in this expression is obviously negative along trajectories. The second term can be studied through the function defined on

R * + × R * + (x, y) → 3 -x -y - 1 xy .
This function is negative if x + y > 3 or if 1 xy > 3. The complementary of this region is closed and bounded, and the function being continuous on this region, it admits a maximum value. Canceling the gradient leads to the only extremum x = y = 1 where the function evaluates to 0. If we now replace x by V V * and y by W 1

W * 1 U * 1 U 1 ,
we have shown that the second term in L is negative along trajectories, with a maximum value of 0 being reached at V = V * ,U 1 = U * 1 and W 1 = W * 1 . The last term remains to be studied, especially the second order polynomials in V . For j ∈ {2 . . . N}, the polynomial in the j-th term of the sum can be expressed (using the choice for η j ):

V 2 + 2 β 2 j x 2 j β j τ j V * (µ j +β j x j ) 2 β j τ j -

β 2 j x 2 j β j τ j -V * V + β 2 j x 2 j β j τ j V * β 2 j x 2 j β j τ j V * +
(µ j +β j x j ) 2 β j τ j -V * (µ j +β j x j ) 2 β j τ j -β 2 j x 2 j β j τ j (µ j +β j x j ) 2 β j x j -

β 2 j x 2 j β j τ j -V * 2 V.
Its discriminant is given by ∆ =4

β 2 j x 2 j β j τ j V * (µ j +β j x j ) 2 β j τ j -

β 2 j x 2 j β j τ j -V * 2 × β 2 j x 2 j β j τ j V * - β 2 j x 2 j β j τ j V * + (µ j + β j x j ) 2
β j τ j -V * (µ j + β j x j ) 2 β j τ j -β 2 j x 2 j β j τ j , = -4

β 2 j x 2 j β j τ j V * (µ j +β j x j ) 2 β j τ j - β 2 j x 2 j β j τ j -V * 2 (µ j + β j x j ) 2
β j τ j -V * (µ j + β j x j ) 2 β j τ j -β 2 j x 2 j β j τ j .

For the same reasons as before, this quantity is negative. The polynomial thus has no real roots, and is always positive. This proves that each term in the sum in L is positive. Overall, we have shown that every term in derivative is negative along trajectories and they cancel out at the point of the disease steady-state. This concludes the proof of the global stability of this equilibrium.

Fig. 3 :

 3 Fig. 3: Bifurcation diagram obtained for the steady-state subunit populations S when the monomer source, λ , varies for two values the polymerization rate τ. The other parameters are set as inTable 1 (Strain 1). The four different equilibria are depicted (disease-free equilibrium, endemic equilibrium, higher subunit only equilibrium (S + ), lower subunit only equilibrium (S -)). A thin line indicates an unstable equilibrium, whereas a bold line indicates a locally stable equilibrium

  Fig. 6: Numerical simulation using the nucleated polymerization model (a) and the template assistance model (b). The parameters used are described in Table 1 (Strain 1). The initial condition is given by V (0) = λ γ = 200,U(0) = 0.1, P(0) = (2x 0 + µ β )U(0), S(0) = 0.2

  Fig. 9: Illustrations of typical cases from the diagram shown in Figure 8, representing the time evolution of the two strains (blue line is Strain 1, orange line is Strain 2). The parameters used are shown in Table1. In each case, V is initialized to the disease-free steady-state λ /γ, and the strains are initialized with the specified conditions. The top panel (a) shows takeover of Strain 1, the middle panel (b) shows takeover of Strain 2, whereas the bottom panel (c) shows asymptotic coexistence of both strains.

Table 1 :

 1 Parameter definitions and values used for numerical simulations (unless specified otherwise). The values are chosen arbitrarily, with magnitudes consistent with values from[START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF], where time is expressed in days. The aggregate-associated quantities are in fibril number per volume unit

	Parameter	Definition

Acknowledgements We thank Human Rézaei (INRA, France) for interesting discussions and careful review, and all the reviewers for relevant comments.

Funding: This research was sponsored in part by the Joint DMS/NIGMS Initiative to Support Research at the Interface of the Biological and Mathematical Sciences (R01-GM126548) and the NSF-INSPIRE program (BCS-1344279). This work was supported by the Association France-Alzheimer, SM 2014. Conflict of Interest: The authors declare that they have no conflict of interest.