Mathematical Modeling of Prion Disease

Abstract

The prion hypothesis, once a heretical violation of the central dogma of molecular biology, has become an accepted mechanism used to explain a host of progressive neurodegenerative diseases in mammals and heritable phenotypes in yeast. From the beginning, mathematical models have been an essential tool in studying prion and other protein misfolding/aggregation processes. In this work, we review some of the major mathematical studies that have contributed to our understanding of prion disease and discuss trends in current and future studies.

1. Introduction

In the past century, the use of mathematical models to study biological phenomena has gone from an occasional dalliance of a theoretical mathematician to an established field of its own. Today, mathematics has impacted virtually every area in biology—from evolution (e.g., Fisher's Fundamental Theorem of Natural Selection) to biochemistry (e.g., Michaelis‐Menten Kinetics) [1]. But, the impact of biology on mathematics has been just as transformative and biology itself has served to motivate the development of novel mathematics [2].

In the latter part of the twentieth century, both biologists and mathematicians worked to identify and characterize mechanisms to explain a host of fatal neurodegenerative diseases in mammals ranging from scrapie—an infectious diseases observed in sheep—to fatal familial insomnia—a genetic disorder in humans. Initially, much of the focus of these studies centered on first the identification of the infectious agent of these diseases. The discovery of the prion—a proteinaceous infectious particle—originally represented a fundamental contradiction in the central dogma of molecular biology. But today there is increasing acceptance of protein‐only‐inheritance (see Figure 1) not only for mammalian diseases but also for heritable phenotypes in yeast. At present, mammalian prion diseases are untreatable and continued experimental, mathematical and interdisciplinary research offers the promise for identification of regulatory mechanisms and therapeutic targets.

Figure 1.

Prion and protein‐only inheritance. The central dogma of molecular biology stipulated that genotype (DNA) encodes phenotype (visible traits). However, prion proteins represent an important departure from this rule where inheritance may arise from proteins alone. Through adopting a stable misfolded conformation (square) a protein can go from harmless to capable of conferring a number of fatal, progressive neurodegenerative diseases.

Prion diseases offer a particularly intriguing biological phenomenon for mathematical analysis because such diseases cover many different systems and time scales. At the level of a population, such as a herd of sheep or population of deer, prion disease can be studied as a classical epidemic model where infections are spread among an initially uninfected (susceptible population). Prion disease can also be studied as a genetic disease whose phenotype is caused by a gain of function mutation in the gene coding for Prp. While the age of disease onset and death appear to be heritable, linking genotype to phenotype remains challenging [3]. Spontaneous prion disease is thought to be nucleation limited, with the formation of a stable minimal size aggregate (nucleus) of misfolded protein serving as the rate‐limiting step in the appearance of prion diseases. All prion diseases are characterized by aggregates of misfolded protein serving as templates to convert normally folded protein and amplifying through fragmentation. As such, many mathematical formulations have focused primarily on the dynamics of the aggregates themselves through modeling either discrete or continuous sizes using ordinary differential equations (ODEs) or partial differential equations (PDEs), respectively. Finally, in order to model the loss or reversal of the prion phenotype in certain experimental systems, prion dynamics are modeled as a stochastic process.

This chapter reviews the application of mathematical models to the study of prions. Our goal is to serve as a tool for both mathematicians and biologists interested in interdisciplinary research in prion disease. We first describe the time before the identification of the prion, when the work of mathematician Griffiths was central to proposing a protein‐only disease process. We next overview mathematical formulations focusing on the dynamics of prion disease through modeling the kinetics protein misfolding and aggregation as well as the coagulation and fragmentation dynamics of the misfolded aggregates themselves. We close by discussing recent advances and ongoing work in mathematical modeling of prions that are serving to further our understanding and motivate experimental studies and present some open questions.

2. Mysterious mammalian diseases, heritable yeast phenotypes and the mathematical origins of the prion hypothesis

No discussion of prion disease would be complete without discussion of the field prior to the establishment of the prion hypothesis, which stipulates that protein, rather than virus or bacteria, is the infectious agent of the prion disease. Here we give an overview of historical observations linking a variety of diseases in mammals, leading to the formulation of the prion hypothesis by Griffith [4], subsequent experimental validation by Prusiner [5] and discovery of prions in fungi. (For a more complete history of prion diseases refer to any of these reviews [6–10].)

Scrapie is likely the first prion disease to be observed with reports dating back to the 1500s [7, 11]. The first publication describing scrapie appeared in 1759 [12] and because scrapie was reported to be an infectious “distemper” from which sheep could never recover, shepherds of the time were advised to separate any animal observing symptoms from the rest of the flock. Publications at this time discussed and debated possible modes of transmission for this disease; ideas were wide ranging from inbreeding [13], humidity of the sheep pen [14] and even atmospheric events [15]. By the late nineteenth century, it was strongly believed that scrapie was a hereditary disease, but some reports noted spontaneous occurrences leading some to believe there was two forms of scrapie: hereditary and non‐hereditary [7].

In 1913, Sir Steward Stockman published “Scrapie: An Obscure Disease of Sheep” [16], which served as both a historical record of the disease as well as analysis of its symptoms and progression. In particular, he notes that scrapie has a long incubation time of 2–3 years. Research on the method of transmission of scrapie continued and by the early 1960s, it had been established that scrapie could spread through indirect contact between sheep (grazing in a field that had been occupied by an infected herd) [17], could transmit either as an infectious or heritable disease [18], could be transmitted through serum as when a vaccine for another disease (Louping‐ill) was prepared from sheep infected with scrapie [19] and could transmit between species (from sheep to goat [20] and sheep to mouse [21]). In combination, these observations suggested that scrapie did not behave as any previously observed disease‐causing agent.

Scrapie was not the only prion disease studied in the mid‐twentieth century. Around the same time that cross‐species infectivity of scrapie was demonstrated, researchers were studying kuru, a progressive neurodegenerative disease that appeared in Papua New Guinea. The disease was first reported in the scientific literature when anthropologists [22] and pediatricians [23] reported a deadly disease among the Fore people called kuru. The disease had an unusual distribution by sex and age; among children, both male and female could have the disease, but among adults incidence was nearly always limited to females [9]. It was also observed that the pathology of kuru was similar to a Cruzfeld‐Jacob disease, a very rare neurological disorder [24].

Researchers continued to conduct experiments to uncover the method by which scrapie and kuru were transmitted. In 1959, a critical connection was made between these seemingly separate disorders; Hadlow, a veterinarian, attended an exhibit at the Wellcome Medical Museum in London featuring images of neurological tissue from the brains of individuals who died from kuru. He noted the patterns and appearance of damage was extremely similar to what he had seen in scrapie. The similar pathology, combined with the apparent ability of kuru and scrapie diseases to be acquired or hereditary caused him to conjecture that a similar mechanism could be responsible for both diseases and advised researchers to see experiment with transmission of kuru from humans to other mammals (as had been done for scrapie) [25]. Indeed, soon after Hadlow's publication it was shown that, like scrapie, kuru could be transmitted to other mammals [26, 27].

While linking diseases such as scrapie, kuru and Cruzfeld‐Jacob was significant in formulating the prion hypothesis, it did not directly address the question of the infectious unit of the disease. In 1966, Alper and colleagues used radiation and filtration experiments on brains from mice scrapie and determined the infectious agent of scrapie appeared to be able to self‐replicate but without a nucleic acid code; they conclude by indicating the scrapie agent “is likely to be of an unusual nature” [28].

In 1967, Griffiths, a mathematician at Bedford College in London took the observations from Alper [28] and Pattison [29] and suggested the infectious agent of scrapie was “probably a protein without nucleic acid” [4]. While precise mathematical formulations were not given, Griffiths used the same type of reasoning that goes into the development of mathematical models to pose three possible mechanisms by which a host‐encoded protein could act as an infectious agent. Namely, he worked within the known rules of the underlying biological processes to pose hypotheses, which could then act to motivate further experimental design. It is precisely this form of interplay between the mathematical and biological sciences that serve to drive discovery.

It is worth noting that Griffith's proposed mechanisms for a protein infectious unit involved three distinct biological processes: gene regulation, protein aggregation and immune response. Because his second mechanism is closest to what we believe to be correct today, we postpone its discussion. First, he suggested a process by a gene encoding the prion protein was typically in the “OFF” state. If the prion protein was capable of acting as an inducer to this gene (i.e., turning it “ON”), then the introduction of prion protein would act infectiously by turning the gene “ON” and further production of the protein would maintain the gene in the “ON” state. As such a prion disease could occur spontaneously if the gene were perturbed to the “ON” state in an individual or be acquired through consumption of a protein. His third mechanism was one where the immune response (antibody) was itself equivalent to the foreign body (antigen) and thus prion disease could be the hosts immunity backfiring.

Remarkably, his second proposal quite closely depicts the dynamics of protein aggregation and fragmentation that today we believe was that of protein aggregation and fragmentation. He posed a simple model where proteins could exist as monomer, dimer, trimer and tetramer. Increase in size could occur through monomer addition and tetramers could split into two dimers. If the reaction to create a dimer from two monomers was itself required the catalytic influence of a dimer, the all‐monomer state would persist stably unless a dimer were introduced. Such a system he noted would be capable of self‐propagating as long as there were monomer (which could be produced by the host) and an initial infectious unit (a dimer, trimer, or tetramer).

Griffith's proposed “protein‐only” method of disease transmission spurred further experimental studies. Finally, in 1982, Prusiner demonstrated through several distinct lines of evidence (including sensitivity to proteases) that the infectious agent was a protein and coined the term “prion” to mean proteinaceous infectious particle [5]. Not long after, a team of researchers discovered the host gene coding for the prion protein, named PrP for prion protein, in mammals [30].

While mammalian disease was the driving force behind the investigations so far discussed, mammals are not the only organisms that today we know to exhibit protein‐only inheritance. In 1994, Wickner was investigating a heritable phenotype in yeast that did not appear to have a chromosomal determinant, but was associated with an altered form of a yeast protein Ure2p [31]; he proposed that this phenotype could be prion based. Thus, the prion hypothesis could plausibly explain a number of non‐Mendialian phenotypes discovered and studied by Cox [32]. The facility of yeast as an experimental system has spurred the identification of nearly a dozen prion proteins in yeast each of which is linked with a seemingly harmless phenotype [6, 33]. Thus, this opens the possibility that protein‐only inheritance may well have evolved as a regulatory mechanism.

While today there remain some scientists that reject the notion that a host‐encoded protein could be the infectious agent, increasingly sophisticated experimental studies continue to support the prion hypothesis. For example, in 2013 Zhang and colleagues demonstrated that prion diseases could be induced in mammals from recombinant prion protein produced in bacteria [34]. As such, the prion hypothesis has become the accepted view for both mammalian prion diseases and heritable yeast phenotypes.

Today we understand that proteins capable of propagating through a protein‐only mechanism do so by adopting an abnormal folded‐state (conformation) and forming aggregates each of which may act as a template to induce further misfolding among normally folded protein. (Note that we use the term “prion phenotype” to encompass both the concept of mammalian prion disease and harmless prion phenotypes in yeast.) Indeed, there are multiple possible prion phenotypes (in mammals these correspond to distinct incubation periods for disease symptoms) each of which corresponds to a distinct conformation typically called a prion strain. Finally, while all known mammalian prion phenotypes correspond to the same protein PrP, in fungi there are a number of prion proteins each linked to distinct phenotypes [6, 33]. However, as we will detail further, identification of this infectious agent is only the beginning in characterizing these processes.

3. Establishing a mathematical framework of prion aggregate dynamics

In this section, we discuss contributions of mathematical modeling in understanding the dynamics associated with prion disease (more generally phenotype). Because prion phenotypes can be either spontaneous or acquired, a distinction is often made between nucleation, the spontaneous appearance of an initially infectious unit and propagation of the infectious unit. However, since both phases involve aggregation of misfolded protein, similar mathematical formulations have been applied to both processes. Indeed both processes are also fundamental to other protein aggregation processes and disorders such as Alzheimer's and Parkinson's diseases. Because several reviews exist on mathematical models of aggregation in more general biological processes [35, 36], in this work we focus on mathematical methods of appearance and propagation as specifically applied to the in vivo dynamics of prion phenotypes. Although in spontaneous prion disease, nucleation occurs first, we will begin our discussion with propagation as this step has been better characterized.

Propagation. The first mathematical formulation of the autocatalytic propagation of prion aggregates was published by Eigen in 1996 [37] where, inspired by the dimerization process expressed in Griffiths’ third hypothesis [4] and observations from by Prusiner [38] and Lansbury [39–41], he developed systems of differential equations to analyze two theories on protein‐only amplification. Through his mathematical analysis, Eigen was able to demonstrate support for the idea that prion aggregates are themselves the infectious agent of prion disease but, as Eigen writes “aggregation of the prionic form is most probably a necessary, but not possibly sufficient, prerequisite of infection”.

In Eigen's first model, he explores the possibility suggested by Prusiner [38] that heterodimers act to template misfolding. He considers a system with two‐protein species: $A$, normal conformation and $B$, prion conformation; proteins of type $A$ are capable of forming heterodimers with proteins of type $B$ and through that interaction are irreversibly converted to type $B$. The resulting homodimer of $B$ would then resolve creating two proteins of state $B$, each of which may then act to template further conversion events. (Note that in this model the capacity of the system to convert protein from state $A$ to state $B$ depends linearly on the total concentration of $B.$)

The mathematical model resulting from these assumptions consists of two coupled differential equations. Eigen performed steady‐state analysis to determine the possible asymptotic concentrations of each protein species and how the local stability of each depended on the underlying kinetic values. He found two types of asymptotic behavior were possible and the one the system would converge to depended on the ratio of two kinetic parameters: the catalytic conversion rate and the death/decay rate of the prion conformation of the protein. If the death rate exceeded the conversion rate, the asymptotic concentration of prion proteins (type $B$) approaches 0 and nearly all the protein will be in then normal conformation (type A). When the conversion rate exceeds the death rate, the reverse happens, namely the amount of protein in the prion conformation (type $B$) will grow exponentially and most of the protein present is in the prion conformation. Since neither of these possibilities was consistent with the true behavior of prion disease, namely that the vast majority of individuals have primarily healthy protein and, even in the few individuals that do have prion diseases, still have detectable levels of normal protein. As such, Eigen concluded that a model where the conversion capacity was linear with the concentration of misfolded protein was not possible [37].

Eigen's second model, considered two mechanisms where the infectious agents were not individual misfolded protein monomers: a cooperative auto‐catalytic mechanism, which generalized his first model and aggregates of misfolded protein, in accordance with a proposed aggregation mechanism from Lansbury [39–41]. These assumptions result in their own—more complicated—sets of differential equations, but as for the previous model, steady‐state analysis revealed important properties of the asymptotic dynamics. Both models exhibited a “threshold” effect, that is, if the concentration of prion protein were low enough, the healthy state was maintained but the introduction of prion protein exceeding a threshold would cause the exponential growth of prion protein. While the results of Eigen's work did not definitively detail all necessary steps in the propagation of prion phenotypes, nor did he demonstrate global asymptotic stability of the prion phenotype, his work demonstrated that mathematical modeling—in particular systems of deterministic ODEs—could be used to theoretically interrogate biological hypotheses on prion dynamics. In particular, Eigen's analysis demonstrated that “aggregation is necessarily involved” [37] in prion propagation.

In 1998, Nowak and colleagues built upon Eigen's seminal work by incorporating additional experimental observations, in particular work demonstrating sensitivity of distinct Prp strains to protease cleavage. Their mathematical framework of prion infection dynamics was based on having prion aggregates act in two ways; first (as in Eigen's model) they would template additional misfolding, but now aggregates themselves could increase fragmentation [42]. Because this model forms the basis of most subsequent mathematical models on prion dynamics, we discuss its formulation in some detail. In this mathematical formulation, the state of the system at time $t$, is the concentration of proteins in the normal conformation, $x\left(t\right)$ and prion aggregates of every discrete size $i$, $y_i\left(t\right)$. They assume protein in the normal conformation is created at rate $\lambda$ and decays at rate $d$, aggregates of all sizes decay at rate $a$. Conversion occurs through contact between aggregates and normal conformers at a rate depending on the size of the aggregate, ${\beta }_i$. Finally, the total number of aggregates increases through fragmentation; in their most general formulation they specify the rate that aggregates of size $j$ fragment to create an aggregate of size $i$ as $b_{j,i}$ and that during fragmentation no mass is lost (i.e., if an aggregate of size $j$ is always fragmented into two aggregates of size $i$ and $(j-i)$). Translating these biochemical kinetic assumptions into a set of differential equations results in the following infinite system:

$\frac{dx}{dt}=\lambda -dx(t)-{\displaystyle \sum _{i=1}^{\infty }{\beta }_{i}x(t)y_i(t),}$E1

$\frac{d{y}_i}{dt}={\beta }_{i-\text{1}}x(t)y_{i-\text{1}}(t)-{\beta }_{i}x(t)y_i(t)-a{y}_i(t)+{\displaystyle \sum _{j=i+\text{1}}^{\infty }(b_{j,i}+b_{j,i-j})y_j(t)}-{\displaystyle \sum _{j=\text{1}}^{i-\text{1}}{b}_{i,j}{y}_i(t)},$E2

for $i=\text{1,2,}\dots$, etc. While the model allows for quite general dynamics, under the simple assumptions that the conversion rate is independent of aggregate size, that fragmentation increases linearly with aggregate size and that fragmentation is equally likely between any two adjacent monomers in an aggregate, this infinite system of differential equations reduces to the following three‐dimensional system:

$\frac{dx}{dt}=\lambda -dx(t)-\beta x(t)Y(t)$E3

$\frac{dY}{dt}=bZ(t)-(a+b)Y(t)$E4

$\frac{dZ}{dt}=\beta x(t)Y(t)-aZ(t)$E5

where $Y(t)={\sum }_{i=1}^{\infty }{y}_i(t)$ represents the total number of aggregates and $Z(t)={\sum }_{i=1}^{\infty }i{y}_i(t)$ is the total amount of prion protein. We note that mathematically $Y\left(t\right)$ and $Z\left(t\right)$ correspond to the zeroth and first moments of the distribution of aggregate sizes and, as such, this demonstrates a moment closure of the aggregate size distribution. That is, the time‐evolution of the complete aggregate size distribution under these kinetic simplifications is determined by purely the zeroth and first moments. Nowak and colleagues remarked this reduced formulation was mathematically equivalent to prior viral models studied in mathematical epidemiology and derived an expression for the basic reproductive number of a prion aggregate. The basic reproductive number, or $R_{\text{0}}$ as is commonly denoted in the epidemic community, specifies the number of secondary infections (in this case infectious aggregate) created by an infection aggregate during its lifetime. In the case that $R_{\text{0}}>\text{1}$, we expect exponential growth of disease in a purely susceptible population and, as such, prion aggregate to persist stably. If $R_{\text{0}}<\text{1}$, we expect the infectious elements, in this case prion aggregates, to exponentially decay and ultimately be lost from the system. In this case the $R_0$ was shown to be a ratio of the underlying kinetic parameters: $R_{\text{0}}=\frac{\beta \lambda b}{da(a+b)}$. As such, the stability of prion phenotypes was now shown to be explicitly a function of biochemical properties offering the promise to interpret results in this new context.

Nowak and colleagues [42] were also the first to formalize what today is considered to be the standard prion aggregate kinetics, the nucleated polymerization model (NPM). In this model, the infectious units are aggregates above a critical size. Below this critical size, aggregates of the misfolded prion form of the protein are presumed to be highly unstable and are rapidly resolved into monomers (see Figure 2). (It is this nucleation process that forms the rate‐limiting step in the establishment of prion phenotypes and we discuss this extensively in the next section.) The dynamics of the NPM are similar to those presented in Nowak's first model; however, the minimal nucleus size modifies the resulting equations slightly. First, the quantities $Y\left(t\right)$ and $Z\left(t\right)$ now represent the aggregates above this critical minimal size, $n_{\text{0}}.$ That is,

$Y(t)={\displaystyle \sum _{i=n_0}^{\infty }{y}_i(t)}\text{and}Z(t)={\displaystyle \sum _{i=n_{\text{0}}}^{\infty }i{y}_i(t).}$E6

Figure 2.

Nucleated polymerization model of prion dynamics. This demonstrates the key steps in the nucleated polymerization model (NPM) of prion aggregate dynamics. (The description of the kinetic parameters is in the text.) This model is characterized by prion aggregates below a critical size $n_{\text{0}}=\text{2}$ (the nucleus size) resolving to protein monomers in the normal folded state.

Under the previous simplifications on kinetic rates, this changes the resulting moment closure of the infinite system of ODEs to the following three‐dimensional system of ODEs:

$\frac{dx}{dt}=\lambda -dx(t)-\beta x(t)Y(t)+b(n_0)(n_0-1)Y(t),$E7

$\frac{dY}{dt}=bZ(t)-(a+b(\text{2}{n}_{\text{0}}-\text{1}))Y(t),$E8

$\frac{dZ}{dt}=\beta x(t)Y(t)-aZ(t)-b(n_0)(n_0-1)Y(t).$E9

In this new formulation the basic reproductive number of a prion aggregate now also depends on the minimal nucleus size $n_{\text{0}}$. This form of the NPM has become the standard approach for modeling prion aggregate dynamics and inspired many future mathematical studies.

In 1999 Masel, Jensen and Nowak conducted an extensive analysis of the NPM [43]. In particular, they sought to link experimental observations on the time to appearance of prion disease symptoms with the kinetic parameters of the NPM. Among other contributions, Masel and colleagues determined a viable range of minimal nucleus sizes, $n_{\text{0}}$. Overall, there was remarkable consistency between parameters predicted from different experimental data sets analyzed providing support at the time for this mathematical formulation. In addition, Masel et al. [43] (and then Greer and colleagues with a generalization [44]) demonstrated that the dynamics of aggregates under the NPM are consistent with the long‐incubation time observed for prion phenotypes. If prion disease begins with the introduction of a small amount of prion protein (in the form of aggregates) those aggregates will first have to increase in size until there are enough fragmentation sites to permit aggregate amplification through fragmentation.

In early twenty‐first century, mathematicians continued formalizing the NPM. Prüss and colleagues [45] demonstrated that the prion phenotypes were globally asymptotically stable and not merely locally stable, through deriving a Lyapunov function. Engler et al. [46] analyzed the well‐posedness of the generalization of the NPM where aggregate sizes were continuous, instead of discrete. As such, rather than an infinite system of ordinary differential equations, the system consisted of a single ODE for protein in the normal configuration and a PDE specifying the distribution of aggregate sizes. While this formulation departs from the physically discrete nature of aggregates, in the limit of large aggregate sizes these formalisms are provably equivalent [47] and the use of PDEs permits a wider array of mathematical techniques. Most notably, the continuous relaxation on aggregate sizes has permitted determination of the explicit asymptotic density [44, 46]. (In comparison, the asymptotic density for the aggregate model with discrete aggregate sizes, while first approximated in 2003 by Pöschel et al. [48], was derived only recently by Davis and Sindi and required special functions [49].)

While today mathematical models of prion aggregate dynamics have been formulated under many more general kinetic assumptions (see [47, 50–52] for example) most of these models have been compared to only in vitro aggregation studies. For yeast, in vivo comparisons have been made for the Sup35/[PSI⁺] prion system [53–55], but linking experimental outcomes uniquely to specific kinetic parameters remains challenging.

Nucleation. As mentioned in the previous section, the rate‐limiting step in prion phenotypes is thought to be the time to the appearance of stable nucleus, that is, an aggregate of misfolded protein that persists stably. (It is typically thought that this nucleus corresponds to a misfolded aggregate of a minimal stable size [42, 43].) The self‐assembly of particles into aggregates is fundamental to many physical, chemical and biological processes. Such a process is referred to in statistical physics as nucleation and mathematical models of nucleation have been studied for nearly a century [56]. In contrast to other biochemical models of protein aggregation, the spontaneous appearance of a prion nucleus is thought to be rare [57, 58]. As such, mathematical models of prion appearance are often framed as first‐passage processes; that is, these models focus on determining the amount of time until a critical event occurs, in this case the appearance of a prion nucleus.

One of the earliest models of self‐assembly of particles was proposed by in 1916 by Smoluchowski [59]. He considered the evolution of the density of clusters of discrete particle sizes under the assumption that clusters of any size could join together (coagulation). In 1935, Becker and Döring introduced kinetic equations for a similar process but where clusters could only change in size through monomer addition or removal [60]. More generally, models of particle self‐assembly are distinguished by their associated set of biochemical equations governing the evolution of cluster sizes. As such, the problem of prion nucleus appearance can be framed as: given a set of biochemical equations governing misfolded protein aggregate formation, determine the time it takes for a critical sized nucleus to form [61].

Broadly speaking, two mathematical formulations have been used to describe the time to nucleus formation: deterministic and stochastic. In a deterministic mathematical model, the predictions or model output is always the same for a given input. In such a formulation, the Law of Mass Action is used to convert the set of biochemical equations to a system of ordinary differential equations (ODEs) [62]. For the standard aggregation processes, like the Becker‐Döring process, systems of ODEs have been extensively studied [60, 63, 64]. In these ODEs, the mathematical model output is a continuously varying quantity approximating the concentration or number of aggregates of each possible size. The time to nucleation would then be specified as the time at which the value associated with the critical nucleus size exceeds a threshold value. When the number of total proteins present is large, a deterministic formulation describes the dynamics well; however, when the number of proteins is small, random effects begin to dominate and to capture these effects a stochastic formulation is required [65]. (We note that for in vitro experiments of prion aggregation, when the concentration of proteins far exceeds physiological settings, deterministic models have proven to be consistent with observed quantities [66, 67].)

Stochastic mathematical models allow for the possibility of the same input to produce different output. In this case, the state of the system is given not as a deterministic quantity, but a random variable that can take on different values [68]. For example, given a coin with two sides (heads and tails), the number of times a coin must be flipped until heads appears is a random variable; one might attain heads on the first try or require many trials before heads appears. Because the observed output can change, the quantity of interest is not the specific output but rather its properties. To continue our example, we might wish to know either what the mean (average) number of flips will be required from a fair coin to produce heads and possibly the variance in that quantity. Alternatively, we might wish to know the probability associated with observing any possible outcome (i.e., what is the probability we flip the coin $k$ times before observing a head); this corresponds to a probability density function. For our example of the coin, the number of flips required before heads appears is given as a geometric probability distribution. That is, the probability that $k$ coin flips are required before the first heads is observed is given by: ${(1-p)}^{k-1}p$ where $p$ is the probability of heads on any given trial. For all but simple systems, such as our coin flip example, it is not possible to obtain an explicit formula for our random variable in question. As such, an increasingly sophisticated set of mathematical and computational tools have been employed to aid in such processes.

We note that for nucleation problems, we are interested not in the state of our protein molecules at any particular time, but the first‐arrival time of the nucleus. That is, the time at which the first aggregate of minimal stable size appears. Below we refer to misfolded protein aggregates smaller than the critical nucleus as proto‐nuclei and any aggregate larger than the nucleus size as a propagon. (We note this is consistent with the definition of a propagon as being a prion aggregate capable of transmitting the prion phenotype upon transmission to an environment with normally folded protein [33, 69].)

We will first frame this problem as a continuous‐time stochastic process and then discuss how statistical properties of the first‐arrival time may be computed. (For a detailed discussion of stochastic processes and first arrival times in biological systems, refer to [61, 68].) For simplicity, let us assume that our system consists of a total number of $m$ molecules of our protein in question and that this number remains constant (i.e., no synthesis or degradation). In this case, if we observe the system at any particular time the state of the system consists of the number of protein aggregates of each possible size. If $n_i(t)$ is the total number of aggregates of each size $i$ at time $t$ ($i=\text{1}$ corresponds to monomer), then for all time we have:

${\displaystyle \sum _{i=1}^{m}i{n}_i(t)=m}.$E10

We use $\Omega$ to denote set of all possible molecular configurations and observe that the size of this space increases exponentially in $m$ but, because there is no synthesis or degradation, $\Omega$ is finite. We distinguish between two sets of configurations in $\Omega$; $A$: those configurations with only proto‐nuclei and $A^C=\Omega$_: those where the system has at least one propagon. While not required, it is often assumed that the system begins in the all‐monomer state (i.e., $n_{\text{1}}(\text{0})=m$, $n_i(i)=\text{0}$ for all $i>1$) [61, 70]. At any given time, aggregates of any size may increase or decrease in time through dynamics such as monomer addition, fragmentation and/or coagulation as allowed by the biochemical assumptions (for example, the Smoluchowski [59] or Becker‐Döring [60] assumptions). Refer to Figure 3 for a visualization of the stochastic model. The first‐arrival time is the time that the stochastic process reaches any configuration in $A^C$.

Figure 3.

Stochastic model for prion nucleus appearance. The rate‐limiting step in the appearance of prion phenotypes is thought to be the waiting‐time until the appearance of a nucleus, an aggregate of misfolded prion protein that exceeds a critical size. In a stochastic formulation, the number of prion subunits of each size is tracked in time. As detailed in the text, we are considering a reduced system where the total number of protein molecules ($m$) remains fixed in time.

There are three methods for computing first‐arrival times and their associated moments (mean, variance, etc.) in this stochastic formulation. First, these quantities may be defined directly by analyzing solutions of the chemical master equation (CME). The CME is the first‐order linear differential equation that describes the time‐evolution of the probability of the system to occupy any particular configuration [68, 71–74]. Because the size of the state‐space is exponential in the number of monomers, the CME is computationally intractable for all but very simple systems. Second, computational simulations representing individual realizations of the nucleation process are generated in silico; the mean and moments of the first‐passage times are then calculated from these empirical results [75, 76]. While Monte Carlo approaches are typically easy to code and highly parallelizable, they suffer from slow convergence and in the case of rare events, like nucleation, individual realizations may take arbitrarily long to terminate [77]. Third, heuristics may be used to simplify the dynamics in particular regimes. For example, in a series of studies Chou and colleagues [65, 78–80] approximated the mean first‐arrival time to a critical nucleus for the Becker‐Döring model (only monomer growth or detachment) in for two parameter regimes (strong growth and weak growth) by computing the arrival time for the dominant pathway from the all‐monomer state to the appearance of the first aggregate of a minimal stable size. While dominant path approaches are readily apparent for some models of aggregation, like Becker‐Döring, they are difficult to determine for more general sets of reactions. Further, because these results rely on particular parameter combinations, the approximations pose challenges to parameter inference—where we want to determine the kinetic parameters that best match available data.

We note that beyond the mathematical challenges in modeling nucleation there remain many practical challenges. Experimentally, it is typically not possible to separate the spontaneous appearance of a propagon, an initial infectious aggregate, from the prion phenotype itself. Finally, critical events in the underlying biochemical kinetics of nucleus formation in prion disease are unknown. In our formulation above, we described protein subunits, but the protein itself is only capable of aggregating when in a particular conformational state. As such, an accurate predictive model of spontaneous nucleation must also include a model of protein misfolding. While it is clear that particular prion variants (distinct conformations) are favored under particular experimental conditions [81], the connection between nucleus formation and conformation has yet to be fully explored. As such, nucleation remains a challenge on experimental, mathematical and computational fronts.

4. Present state and challenges in prion disease modeling

As described above, a combination of mathematical and experimental studies over the past few decades have led to the formulation of a protein‐only form of inheritance associated with prion phenotypes. We first summarize our present knowledge and then remark on present day studies and challenges that remain in modeling prion disease. Today we believe that prion phenotypes are established through two distinct phases, nucleation and amplification. Once an initial nucleus—prion aggregate above a critical size—is introduced to a host, four steps are required for successful in vivo propagation of prion phenotypes (see Figure 4). First, normal folded protein must be continuously created. Second, aggregates of the misfolded form of the protein act as templates by converting normally folded protein to the same misfolded state. Third, the total number of templates increases through fragmentation where a single aggregate is split into two (or more) smaller aggregates. Finally, misfolded protein must spread through other cells. For yeast, this transfer of misfolded protein occurs through cell division where for mammals this likely involves extracellular diffusion [33].

Figure 4.

In vivo prion propagation. The infectious agent of prion disease is aggregates of misfolded proteins (squares). Four steps are essential to stable propagation of the prion form of the protein: (1) new normally folded protein must be created, (2) prion aggregates act as templates to convert normally folded protein to the misfolded state, (3) prion aggregates are fragmented into smaller aggregates each of which may be capable of acting as a template, and (4) prion infectious units must spread through other cells. In the case of yeast prions, this transmission occurs through the normal process of cell division while for mammalian disease transmission corresponds to diffusion and transport through tissue.

However, beyond this basic understanding remain many challenges, both mathematical and biological. We briefly outline some open questions in prion biology we believe are amenable to interdisciplinary approaches.

Consideration of the cellular environment. Prion phenotypes are established through protein misfolding, but protein misfolding itself is not rare. Eukaryotic cells have developed a complex network of molecular chaperones and protein degradation factors that act continuously to identify and clear misfolded proteins [33]. As such, understanding the in vivo propagation of prion phenotypes requires considering the environment in which they appear. In the case of yeast prions, the molecular chaperone Hsp104 has been shown to be essential for the propagation of [PSI⁺] prion phenotype in yeast. While comparatively few mathematical studies have considered the role of Hsp104 as enzyme catalyzing fragmentation [55, 82], the results from these studies have resolved previously unsupported results on shifts in aggregate size distributions and (as we will discuss further below) rates of loss (curing) of specific [PSI⁺] strain phenotypes.

Spread of prion aggregates. A major open question in prion biology amenable to mathematical analysis is how prion aggregates spread between cells either through cell division (yeast prion phenotypes) or within a mammalian tissue. Although this question was first explored when Nowak and colleagues [42] presented simulations of their NPM where prion aggregates could move between distinct cells in a population and when Payne and Krakauer [83] considered prions spreading in tissue as a traveling wave, our understanding of the spread of prion aggregates is still incomplete. Because the long incubation period and complicated physical domains involved in mammalian prion disease will pose significant experimental and mathematical challenges, yeast prion phenotypes may provide a useful tool in this question. While stochastic models have been developed for yeast, which link cellular levels of prion aggregates with computational simulations [55], formulations amenable to analytical treatment need to be developed to allow for a more systematic characterization of the spread of prion aggregates.

Reversing prion phenotypes. Mammalian prion diseases remain untreatable and ultimately fatal and as such the identification of clinical treatments and methods of early detection remain important scientific and technical challenges. Of particular challenge may be that drugs which act to promote aggregate fragmentation, with the goal of causing all aggregate to drop below the critical nucleus size, may in fact promote prion amplification at low doses when they would merely act to accelerate the exponential growth of prion aggregates [84].

A promising avenue toward finding approaches to manage prion diseases in mammals, might be to more clearly understand the biochemical processes responsible for a number of reversible prion phenotypes in yeast [33]. As was demonstrated by Derdowski et al. [55], a combination of enzyme‐limited fragmentation and aggregate‐size transmission bias appeared to be responsible for the observed natural rates of curing for the [PSI⁺] weak phenotype. In addition, for yeast treatment with GdnHCl has been shown to significantly slow aggregate fragmentation leading to a natural reversal of prion phenotypes by dilution of the aggregate during cell division [69, 85–88].

Further, biological and mathematical researchers should consider mechanisms of curing prion disease through studying mutations in prion proteins, which are known to slow or halt the disease progression. Such mutations are known to exist in mammals [89] and yeast [90].

Evolvability of prions. While prions were originally implicated in mammalian disease, the fact that they persist as a number of harmless heritable phenotypes in yeast raises intriguing questions about how and why prions may have evolved [57, 91, 92]. Because the yeast phenotype [PSI⁺] is associated with a decreased efficiency in stop‐codon recognition, it is thought to serve as an evolutionary capacitor by promoting the generation of novel transcripts [57]. More recently, a combination of mathematical and experimental studies have demonstrated that smaller [PSI⁺] aggregates still retain function associated with the normal Sup35 conformation offering the possibility this system evolved to tune stop‐codon recognition [54]. While we are still far from understanding the forces behind prion phenotype evolution, the evidence continues to mount for possible beneficial examples of prion‐like mechanisms [93].

5. Conclusion

Many questions remain about prion phenotypes and it is essential once again for scientists with different backgrounds to utilize their disciplinary expertise and methods to address these questions. As we have discussed, two critical points in the history of prion disease came from researchers that were not primarily biologists, namely the mathematician Griffith [4] and the veterinarian Hadlow [25]. If the past is any predictor, future studies in prion phenotypes will continue to benefit from an interdisciplinary approach.