Multi-Scale Mathematical Modeling of Prion Aggregate Dynamics and Phenotypes in Yeast Colonies

Abstract

Prion diseases are a multi-scale biological phenomenon that requires understanding intracellular processes as well as how cells interact with each other and their environment. In mammals, prion diseases are progressive, untreatable, and fatal. Yeast prion phenotypes are harmless and reversible, which suggests a deep understanding of the reversal of prion phenotypes in yeast may be informative to mammalian diseases. In yeast, the loss of some prion phenotypes appears to be stochastic and spatially dependent, suggesting a cell-based model of yeast prion dynamics would be a powerful tool for comparisons with experimental results and hypothesis generation. In this work, we consider the components necessary to develop such a model that depicts both the biochemical-, intracellular-, and colony-level scales in yeast prion phenotypes. We first review the literature of mathematical models of the intracellular processes of prion disease. We then review common approaches to cell-based modeling of multicellular systems and how they have led to biological insights in other systems. This chapter ends with a discussion of future studies aimed at motivating how these two types of models can be coupled to produce multi-scale models of prion phenotypes.

1. Introduction

Prion proteins are most commonly associated with fatal progressive neurodegenerative diseases in humans and other mammals [1, 2, 3, 4]. Intriguingly, all mammalian prion diseases are the result of a single host expressed protein, PrP. In mammalian prion disease, a misfolded form of PrP appears, and rather than being cleared by protein quality control mechanisms, this misfolded form persists and associates into aggregates. These aggregates then act as templates that convert normally folded protein to the misfolded form and may break into multiple aggregates further amplifying the conversion of normal protein [5, 6]. The initiating event of prion disease, the initial appearance of an infectious protein agent, can occur spontaneously, as in sporadic Creutzfeldt-Jakob disease [7]; because of a genetic mutation, as in fatal familial insomnia [8]; and through interaction with another infected host from the same (as in scrapie [9], chronic wasting disease [10], or Kuru [11]) or different species, as in when Bovine spongiform encephalopathy (“mad cow” disease) is acquired by humans [12] or scrapie is spread from sheep to goats [13]. Further, some prion diseases can occur through multiple modes, such as Creutzfeldt-Jakob which has been associated with each of these modes of transmission [14]. More generally, the amyloid structure of prion aggregates is common to non-prion proteins associated with neurodegenerative disorders such as Alzheimer’s, Parkinson’s, and Huntington’s disease, and evidence continues to suggest commonalities between such disorders and prions [15, 16], including their transmissibility [17]. At present, prion disease, and other amyloid disease such as Alzheimer’s [18], has no effective methods for treatment or early detection. However, one potential source to gain insight into these pathogenic phenotypes is to study instances in biology where amyloids confer beneficial consequences, such as growth advantages in yeast and long-term memory storage in mammals [19], and therefore the amyloid dynamics can be studied absent of the disease state [16, 20].

One promising candidate for probing the mechanisms behind prion-associated phenotypes is the yeast S. cerevisiae [21, 22, 23]. In yeast, prion proteins were first identified in the 1980s when interrogating the genetic basis for non-Mendelian phenotypes [21, 23]. Today, nearly a dozen harmless phenotypes in yeast are known to propagate through misfolded protein aggregates [6]. Unlike their mammalian counterparts, yeast prion phenotypes have been associated with many different proteins including Sup35 ($\left[{\mathit{PSI}}^{+}\right]$ phenotype), Rnq1 ($\left[\mathit{PIN}\right]$ phenotype), and Ure2 ($\left[\mathit{URE}3\right]$ phenotype). As for mammalian prions, each of these phenotypes has multiple variants (or strains) which themselves are associated with different properties such as infectivity and phenotype presentation. However, unlike mammalian prion disease, which can take months of years to present, yeast prion phenotypes manifest quickly and, most remarkably, can be made to appear or disappear through mild experimental manipulations which do not impact the host cell [6, 16, 24, 25]. The underlying aggregate dynamics are the same regardless of organism. As shown in Figure 1, the establishment of prion phenotypes in yeast (and of course mammals) is an inherent multi-scale process. First, there is the biochemistry of the interactions between normal and misfolded protein. Then, there is the level of an individual cell where these protein aggregates interact not only with the cellular environment including molecular chaperones and protein degradation factors. Then, these protein aggregates must be transmitted to other cells. Note that in yeast, this transmission occurs between mother and daughter cells during division [26, 27], while in mammals, this transmission typically occurs between cells in a tissue [28, 29]. Finally, prion phenotypes are determined by examining populations of cells. In yeast, prion phenotypes are assayed at the level of a yeast colony, and mammalian prion phenotypes are typically distinguished by their patterns of neurological decay [14].

Figure 1.

Multi-scale dynamics of yeast prion Phenotypes. (A) Prion aggregates are thought to be ordered linear structures (squares) capable of converting normally folded protein (circles) to the prion form and fragmenting into multiple aggregates. Prion aggregates are thought to have a minimal stable size $n_0$ below which the aggregate structure will not be stable. (B) Aggregate dynamics occur within an individual yeast cell (large circle) but are spread to daughter cells through transmission through the bud neck. (C) Phenotypes in yeast are observed at the level of a colony which consists of thousands of individual cells. $\left[{\mathit{PSI}}^{+}\right]$ phenotypes are associated with a color red ($\left[\mathit{psi}-\right]$ or no prion), pink ($\left[{\mathit{PSI}}^{+}\right]$), or sectored (mixture of phenotypes).

Although there is growing appreciation for the multi-scale nature of mammalian and yeast prion phenotypes, the vast majority of mathematical models developed for studying such phenotypes typically focus on only one scale. In particular, the majority of mathematical model have focused on the isolated prion aggregate dynamics. While these models have been particularly useful in comparisons to in vitro prion aggregation, they are limited in their ability to insight to the in vivo systems. In a population of living cells, different cell behaviors such as growth, diffusion, and division are known to impact the abundances and concentrations of reactants and could have a large impact on protein aggregation or more specifically propagation of prion aggregates. Moreover, in the case of yeast, since the proteins are only transmitted during cell division [6] and yeast have an asymmetric cell cycle [30], a colony of yeast with a prion phenotype will exhibit considerable heterogeneity in prion aggregates and potentially other relevant cellular constitutes [27].

One particular yeast prion phenotype that is intimately linked to population heterogeneity and yet to be mathematically characterized is that of colony sectoring for the $\left[{\mathit{PSI}}^{+}\right]$ yeast prion phenotype (Figure 1(C)). The $\left[{\mathit{PSI}}^{+}\right]$ phenotype is associated with the color of a colony: a red colony indicates a $\left[\mathit{psi}-\right]$ (prion-free) phenotype and $\left[{\mathit{PSI}}^{+}\right]$ colony is in a shade of white or pink. The actual color displayed by the colony is associated with a loss of function of the normal Sup35 protein as in when in a prion aggregate the Sup35 monomers have greatly reduced normal function [31]. A sectored colony represents a case when cells in the colony display both $\left[\mathit{psi}-\right]$ and $\left[{\mathit{PSI}}^{+}\right]$ phenotypes. Such a colony is thought to appear when the founding cell contains prion aggregates, but at some point through division of cells and growth of the colony, some cells lose the ability to propagate the prion phenotype. Because spontaneous appearance of Sup35 prion aggregates is a rare event, any subsequent daughter cells will also lack prion aggregates. These sub-colonies associated with the loss of the prion phenotypes are thought to correspond to the red sectors. Since the number and pattern of red sectors is fairly similar under the same condition, this suggests an intimate connection between the lineage of ancestry of the cells which were the first to lose the ability to propagate the prion phenotype. Indeed a previous study that did not consider the spatial distribution of individual cells suggested that loss of prion aggregates for a particular variant of $\left[{\mathit{PSI}}^{+}\right]$ occurred in a lineage-dependent fashion [27]. What is the mechanistic basis for this partial loss? Can these mechanisms be exploited in the context of mammalian prion diseases? Answering these questions requires developing mathematical frameworks that accurately represent all scales of the yeast prion system.

Beyond mammalian prion disease, yeast prions offer a particularly unique opportunity to investigate and model the emergence and role of heterogeneity in cellular populations. In the past, our knowledge of individual cell behaviors was based on ensemble measurements that do not account for effects of cellular heterogeneity within a population. However, emerging technologies such as single-cell gene expression analysis, single-cell genome sequencing, and single-cell imaging technologies [32] provide an unprecedented opportunity to quantitatively measure heterogeneity at the single-cell level. Molecular and live imaging experiments investigating growth and development of multicellular systems provide rich data sets that can be used for the first time to understand the different signaling cascades and components that are necessary for cell behaviors to arise. This new level of insight makes it possible to reveal a more accurate picture of cellular behavior and highlights the importance of understanding cellular variation in a wide range of biological contexts. Developing models that describe heterogeneity throughout an entire population provides a quantitative way to understand the dynamical behavior of heterogeneous cell population characteristics/biological importance of heterogeneity.

In this chapter, we consider the different tools available for developing a mathematical framework for distinct scales of prion dynamics. We first review mathematical models for prion disease dynamics and discuss the few models which have considered multiple scales of prion aggregation. We then discuss cell-based models as a promising approach for coupling intracellular and intercellular scales. We discuss several common classes of cell-based models and emphasize the contributions they have made to the understanding of biological systems. We conclude with a discussion how one could build a multi-scale model for prion disease dynamics in yeast that couples current models for aggregation of proteins with a spatial model of yeast growth and proliferation on the scale of an entire population.

2. Intracellular dynamics

As we discussed in Section 1, prion dynamics is an inherently multi-scale process. However, the bulk of mathematical modeling on prion aggregation has considered only the microscopic scale, intracellular dynamics. Because the formation of the initial stable aggregate is a rare occurrence, most of these models have focused on the dynamics of the prion aggregates themselves in a well-mixed compartment leading to models that reflect what occurs microscopically within a single cell or are taken to be the average behavior of a population or tissue. In addition, these models assume either discrete or continuous aggregate sizes employing ordinary differential equations (ODEs) or partial differential equations (PDEs), respectively. This literature has emphasized establishing results on the global stability of prion aggregates as well as their asymptotic size distribution. More recently, mathematical approaches have been developed to consider the spread of prions between cells in a tissue or even organisms, but only under severely simplified molecular scales. At present, no single multi-scale framework has yet been formulated for prion dynamics leaving many unexplained phenotypes, such as yeast sectored colony phenotypes (see Figure 1(C)).

In this section we review the development of mathematical models in the study of prion dynamics. (Readers interested in a more detailed view on mathematical models of prion disease can consult [33].) Although many of these models share commonality with more general aggregation processes in neurodegenerative disease, we focus only on prion aggregation. (Readers interested in mathematical models of more general neurodegenerative phenomena can consult [34].)

2.1 Intracellular models of prion phenotypes

The prion hypothesis, the idea that proteins themselves could encode information about their structure without a DNA intermediary, was first suggested by experimental studies in 1966–1967 (see [35, 36, 37]). Remarkably, the first mathematical formulation of prion aggregation was published in 1967 when Griffiths [38] offered a simple mathematical autocatalytic process that was consistent with a protein-based infectious agent of scrapie. Since then, a large number of mathematical models have emerged which track the size density of prion aggregates under different biochemical processes but typically include conversion of normal protein to prion aggregates and fragmentation of aggregates. Most are based whole or in part on the notion of the nucleated polymerization model (NPM).

Nowak and colleagues developed the NPM in a series of papers [39, 40] based on a model previously introduced by Eigen [41]. 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)$. Protein in the normal conformation is assumed to be created at rate $\lambda$ and decays at rate $d$, while aggregates of all sizes decay at rate $a$ (see Figure 2) Conversion occurs through contact between aggregates and protein in the normal conformation 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, Nowak et al. 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., an aggregate of size $j$ is always fragmented into two aggregates of size $i$ and $\left(j-i\right)$. Translating these biochemical kinetic assumptions into a set of differential equations results in the following infinite system of ordinary differential equations:

Figure 2.

Nucleated polymerization model. The NPM (see [39, 40]) is the canonical model for prion aggregate dynamics. In this formulation normally folded protein (circles) can be created with rate $\lambda$, decayed with rate d, converted to the misfolded form, and added onto the end of an existing aggregate (linear complex of squares) with rate $\beta$. Aggregates have their own degradation rate, a, that is independent of aggregate size, and number of aggregates increases through fragmentation which occurs at a rate b times the number of fragmentation sites in an aggregate. (Note that in this diagram, we allow normal protein to be recovered if an aggregate is fragmented below the minimal nucleus size ${\mathit{n}}_{\mathit{0}}=\mathit{2}$ depicting the aggregate dynamics evolving according to Eqs. (7)–(9)).

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

$$\frac{{\mathit{dy}}_i}{\mathit{dt}}={\beta }_{i-1}x\left(t\right)y_{i-1}\left(t\right)-{\beta }_{i}x\left(t\right)y_i\left(t\right)-{\mathit{ay}}_i\left(t\right)+\sum _{j=i+1}^{\infty }\left(b_{j,i}+b_{j,i-j}\right)y_j\left(t\right)-\sum _{j=1}^{i-1}{b}_{i,j}{y}_i\left(t\right)$$E2

for $i=1,2,\dots$, since there could possibly be infinitely many aggregate sizes. They then considered two simplifying assumptions that allowed a moment closure. First, the conversion rate is assumed to be independent of aggregate size. This occurs when assuming that the templated conversion of normal protein to the misfolded state can only occur when interacting with either end of a linear aggregate. Second, the fragmentation of an aggregate is assumed to be linearly proportional to its size, and the fragmentation kernel is assumed to be uniform (i.e., fragmentation is equally likely to occur at any monomer junction in a linear aggregate). Under these two simplifications, the infinite system of differential equations has the following three-dimensional moment closure:

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

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

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

where $Y\left(t\right)={\sum }_{i=1}^{\infty }{y}_i\left(t\right)$ represents the total number of aggregates and $Z\left(t\right)={\sum }_{i=1}^{\infty }{\mathit{iy}}_i\left(t\right)$ represents 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, the time evolution under these kinetic simplifications is determined by the zeroth and first moments.

Nowak and colleagues demonstrated this model (Eqs. (3)–(5)) was mathematically equivalent to previously developed viral models studied in mathematical epidemiology and derived a basic reproductive number for their model. The basic reproductive number, or $R_0$, is a common parameter in epidemiology and specifies the number of secondary infections (in this case infectious aggregate) created by an infectious aggregate in an initially susceptible population. In the case that $R_0>1$, we expect exponential growth of aggregates and therefore a stable nontrivial prion aggregate distribution. If $R_0<1$, we expect the prion aggregates to exponentially decay in number and the system to return to the aggregate-free state.

In their final formalization, Nowak and colleagues [39] considered the nucleated polymerization assumption of prion aggregates. Namely, it is believed that aggregates below a critical size are not thermodynamically stable. (This is thought to be part of the reason spontaneous appearance of prion disease is so rare because it is a nucleation-limited process [42, 43].) Since Nowak and colleagues did not consider spontaneous nucleation, the nucleus size enters only when aggregate fragmentation would create an aggregate of size below the minimum stable nucleation size $n_0$. In this case, the monomers in that aggregate return to the pool of normally folded protein. Under the previous simplifications on kinetic rates, this changes the resulting moment closure of the infinite system of ODEs as follows:

$$Y\left(t\right)=\sum _{i=n_0}^{\infty }{y}_i\left(t\right)\text{and}Z\left(t\right)=\sum _{i=n_0}^{\infty }{\mathit{iy}}_i\left(t\right)$$E6

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

$$\frac{\mathit{dY}}{\mathit{dt}}=\mathit{bZ}\left(t\right)-\left(a+b\left(2n_0-1\right)\right)Y\left(t\right)$$E8

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

Eqs. (7)–(9) (or appropriately modified Eqs. (1) and (2)) are referred to as the nucleated polymerization model (NPM). 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_0$.

The NPM, and its variants, has been extensively studied. In their paper, Masel, Jensen, and Nowak conducted an detailed analysis of the NPM [39] including extensive linking of experimental observations on the time to appearance of prion disease symptoms with the kinetic parameters of the NPM and determining a viable range of minimal nucleus sizes $n_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. [39] (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.

Mathematicians continued to formalize the NPM through the twenty-first century. 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 a 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]. Mathematical models of prion aggregate dynamics have been formulated under many more general kinetic assumptions such as nonlinear conversion rates, aggregate joining, and general fragmentation kernels (see [50, 51, 52, 53]). More recently, models have been developed which consider prion aggregates along with other intracellular species (such as molecular chaperones) which impact their biochemical dynamics [54, 55], and others have considered interactions between prions and other protein aggregation diseases such as Alzheimer’s [56].

2.2 Multicellular models

While some of the models discussed in the previous section were not explicit about their biological context, in nearly all cases the mathematical models of prion aggregation processes neglected the multicellular scale of the establishment of prion phenotypes. However, there is an emerging literature on multicellular models of prion disease dynamics. In all cases, the intracellular scale is either simplified or ignored.

2.2.1 Multicellular yeast colonies

Several models have emerged to consider the special case of heterogeneity of prion aggregates in yeast colonies. Such models have benefitted from the extensive work on yeast cell division which has facilitated the cellular scale of prion dynamics. We discuss two approaches taken to model this heterogeneity.

First, a number of models have emerged that allow for direct comparison between yeast “aggregate counting assays” [26, 57, 58, 59, 60]. In such assays, an original cell is sampled, and and it is exposed to GdnHCl which is thought to severely reduce the fragmentation rate without altering the ability of cells to divide. Because the number of aggregates remains (relatively) constant while the number of cells increases exponentially, eventually each cell in the population is assumed to have at most one aggregate. When these resulting cells are restored to normal growth conditions (i.e., aggregate fragmentation is normal) and allowed to find their own independent colony, the number of $\left[{\mathit{PSI}}^{+}\right]$ colonies observed is taken as reflective of the number of aggregates in the original cell. (For the purpose of narrative simplicity, we have greatly simplified both the experimental assay and its conclusions assay. Technically this assay counts the number of “propagons” in a cell, and here we assume that propagon is synonymous with prion aggregate although the true case is more complicated [6, 60].) Because the number of aggregates does not increase in most aggregate counting assays, this allows the intracellular processes to be greatly simplified. That is, the number of aggregates changes only through cell division when aggregates are distributed perhaps according to experimentally observed biases in cell division.

One recent study took a different approach by considering the related aggregate recovery assay using a generation and aggregate structured population model [61]. (Once again, the experimental process is too complicated to fully describe, so we offer a simplified description.) Aggregate recovery assays monitor the number of aggregates in a typical cell in time. This is a two-stage experiment where all cells are driven to have (ideally) at most one aggregate in a single cell. Then the cells are released, and aggregates will thus amplify within each cell, and cells will continue to grow and divide as normal. In this assay, cells from the recovering population are sampled in time, and the number of aggregates counted according to the assay in the previous paragraph. Both for mathematical simplicity and because this experiment only monitors the number of aggregates and not their size of concentration of protein in the normal conformation, Banks et al. [61] considered a simplified model of intracellular dynamics where the number of aggregates in a single cell increased according to either an exponential or logistic growth equation. They developed a structured population model that separately tracked the number of aggregates in cells having undergone 0, 1, 2, etc. cell divisions since the start of the experiment. Although this model greatly simplified the intracellular aggregate dynamics, researchers were able to recapitulate know biological relationships between prion strains.

At present, only two studies have considered both a detailed view of the intercellular aggregate dynamics and multicellular yeast population [27, 62]. Because these models were analytically intractable, both researchers developed a multi-scale stochastic simulation framework. Each cell in the population was modeled as a well-mixed compartment with discrete numbers of proteins in the normal conformation and in each possible aggregate size. The intracellular aggregate dynamics were forward simulated according to the Gillespie algorithm for sampling from the chemical master equation. During cell division aggregates were distributed between cells, and the lineage relationships between all cells were tracked. Tanaka et al. [62] used this stochastic multi-scale framework as a way to validate conclusions that drew from a system of ODEs they developed which considered population averages. Derdowski et al. [27] considered explicit comparisons with single-cell experiments in yeast. They used their model to conclude that in vivo fragmentation occurs in a rate-limiting fashion and that aggregate transmission was biased such that larger aggregates are retained by mother cells. However, neither study considered the spatial growth of cells in a colony and, as such, cannot be adapted to the question of yeast colony sectoring.

2.2.2 Mammalian tissue and organismal-level models

In mammals, mathematical models have been developed to study the spread of prion (and more generally protein aggregation) in the brain and between organisms. In the study of Alzheimer’s, researchers have leveraged knowledge about the connectivity between regions of the brain to develop a network diffusion model of disease progression [63, 64]. Remarkably, they discovered that characteristic patterns of atrophy in patients can be explained by eigenelements of their network model. However, as in the models of aggregate counting assays in yeast, they assumed a highly simplified model of protein dynamics. That is, they did not consider molecular scale processes of aggregation but only the diffusion of an agent through the brain and associated atrophy itself with higher concentrations of the agent. However, because the early phases of mammalian prion and neurodegenerative are poorly characterized in vivo, it is possible that their models are consistent with late stages of neurodegenerative disease.

Finally, because prion disease spreads between animals in the same population, epidemic models have been employed to study the population-level dynamics of these diseases. In particular, researchers have studied the spread of scrapie in flocks of sheep [65, 66, 67] and chronic wasting disease in elk [68, 69]. Although these models have been informative of the specific mechanisms which prion disease is likely to spread in these populations (mother to offspring, indirect environmental transmission), they again consider only the organismal scale and not the molecular mechanisms responsible for within the host aggregation dynamics.

Cell-based models, a type of agent-based simulation framework, offer the ability to consider complex intracellular dynamics, distinguish between individual cells, and track the spatial distribution of cells within a colony. In the next section, we discuss common cell-based model formulations with a focus on the mathematical foundation for these models as well as their contributions to particular biological domains.

3. Cell-based modeling approaches

As described in Sections 1 and 2, prion diseases offer a particularly intriguing biological phenomenon for mathematical and computational analysis because such diseases cover many different systems and spatial scales. At the level of a population, prion diseases can be studied as a classical epidemic model where infections are spread among a susceptible population. Alternatively, prion diseases can be studied on a microscopic scale as a genetic disease whose phenotype is caused by prion aggregates that cause a gain of function mutation in certain genes. Current models of protein aggregation in yeast have successfully provided further insight into important mechanisms driving prion disease dynamics (i.e., conversion and fragmentation), but there is a need to develop models that consider the underlying microscopic processes of protein aggregation together with macroscopic properties of the environment in which they are taking place. This requires modeling frameworks that consider the impact of processes taking place on many different spatiotemporal scales. One such class of models that has been developed primarily for studying biological processes from many different scales is cell-based models.

In this section we review cell-based modeling approaches and their various applications. Cell-based models are mathematical models that represent individual cells as discrete entities and produce simulations to predict large-scale, collective behavior of a population or group of cells from the behavior and interactions of individual cells. The inputs to a cell-based model are experimentally observed cell behaviors including how individual cells respond to both intracellular and extracellular cues. Cell behaviors are encoded in a set of biologically relevant rules for cells to follow during simulations. The outputs of a cell-based model are measurable population-level and cell-level characteristics that follow nontrivially from cell–cell coordination and individual cell response to the changing microenvironment, including changes caused by the cells themselves.

Within this class of models a further distinction is made between lattice-based or cellular automaton (CA) models in which the particles live on the coordinates of a lattice and off-lattice methods that use real numbers to describe the coordinates of each particle. Biologically motivated CA models describe cell-cell and cell-environment interactions through phenomenological local rules, which allow for efficient simulation of many different biological systems ranging from bacteria, excitable media, and chemotactic aggregation to chicken embryonic tissues and tumors. (For reviews see [70, 71, 72, 73, 74, 75, 76, 77].) CA models are well suited for studying the dynamics of systems with a large number of cells, i.e., growing monolayers or tumors, because their computation time is efficient and they allow for systematic sensitivity analysis of different parameters [70, 73, 75, 76, 78, 79]. In addition, CA models provide an easy starting point from which to derive continuum equations for cell densities which are particularly helpful in studying the mechanisms driving pattern formation and the growth process of multicellular systems with billions of cells [73, 78, 80, 81]. However, CA models are rule-based, and model parameters do not represent a direct physical description of biophysical properties of cells [73, 78]. For this reason, we focus our discussion in this chapter on lattice-free modeling approaches.

In lattice-free modeling approaches, individual cells are represented using either a single particle or agent or a group of particles or agents that are free to move to any location in the computational domain [75, 76, 78, 82, 83, 84, 85, 86]. In simulations, each biological cell is modeled as a discrete entity and endowed with well-defined individual characteristics such as intracellular reaction kinetics or detailed biophysical properties. As stated by T. Newman, cell shape and cell response to local mechanical forces is one of the most challenging biological characteristics to build into a cell model [78]. For this reason, lattice-free modeling approaches were developed as a way to explicitly model cell shape as well as investigate how individual cells respond to local mechanical forces and interactions with neighboring cells.

In lattice-free modeling approaches, interactions between agents are described using forces or potential functions, and the position of an individual cell is determined by solving an equation of motion for each agent that belongs to that cell. Moreover, it is standard to make the following two simplifying approximations for the equation of motion of each agent: (1) motion may be described by considering each vertex to be embedded in a viscous medium that applies a drag force on it with mobility coefficient $\eta$, and (2) inertia is vanishing [83, 87, 88, 89, 90]. This leads to first-order dynamics with the evolution of the position ${\mathbf{x}}_i$ of agent $i$ to be determined by

$$\eta \frac{d{\mathbf{x}}_i}{\mathit{dt}}={\mathbf{F}}_i$$E10

where ${\mathbf{F}}_i\left(t\right)$ denotes the total force acting on vertex $i$ at time $t$. The mobility coefficient $\eta$ determines the timescale over which mechanical relaxation occurs and is often referred to as the damping coefficient.

Lattice-free modeling approaches vary in complexity from those that track the center of each cell as a single point [91, 92, 93] to those that represent individual cells as collections of membrane and cytoplasm nodes to allow for more biologically relevant, emerging cell shapes [89, 94, 95]. Although a main advantage of CA methods has been their computational efficiency, as computational power continues to increase, lattice-free modeling approaches offer equally powerful tools for realistically modeling biological cell behaviors such as shape change, polarized cell growth, division, differentiation, and intercellular signaling dynamics in great detail without adding computational complexity. Lattice-free models have been used to investigate the impact of individual cell behaviors and make predictions about mechanisms controlling many different biological and pathological processes such as determining shape of a tissue, size of a cell colony, or biological function of an organ [75, 76, 78, 82, 84, 86, 96, 97]. In this section, we review three main categories of lattice-free modeling approaches, namely, center-based models, vertex models, and subcellular element models, as well as provide examples of how lattice-free modeling approaches have been used to study various biological problems.

3.1 Center-based models

Center-based modeling (CBM) approaches, or cell center models, track the center of each cell as a single point and can be classified by two main components, a definition of cell connectivity (Figure 3) and a definition of how cells interact or a force law (Figure 4) [76, 78, 83, 91]. There are several commonly used methods for determining cell connectivity, including the overlapping spheres approach, triangulation, or Voronoi tessellation as described in Figure 3.

Figure 3.

Determining cell connectivity in center-based models. (A) In the overlapping spheres approach, each cell has an intrinsic radius, R_i, and two cells interact if their corresponding spheres overlap, i.e., if they are within a specified fixed distance from each other. (B) Triangulation involves computing a triangulation of the simulation domain using cell centers as nodes. Edges of the resulting mesh determine which cells are connected. (C) The Voronoi tessellation determines a region for each cell defined as a set of points in the plane that are closer to that cell center than any other.

Figure 4.

Center-based model (CBM). The CBM approach updates cell positions based on interactions with connected cells. We show an interaction between two cells where each cell is considered as a sphere with a given radius (although other types of connectivity are possible (see Figure 3)). The left hand side shows the position of node $r_i$ at time $t$ with corresponding radius $R_i$ and position of node $r_j$ at time $t$ with corresponding radius $R_j$. The position of each node is used to calculate the overlap ${\delta }_{\mathit{ij}}$ between cell $i$ and cell $j$ and the resulting linear force $F_{\mathit{ij}}$ to be applied to each node. The equation used to update the position of node $r_i$ is given as a function of the force $F_{\mathit{ij}}$, time step $\mathrm{\Delta }t$, and damping coefficient $\gamma$. The right hand side shows the new positions of node $r_i$ and $r_j$ at time $t+\mathrm{\Delta }t$ after each position has been updated. (See Section 3.1 for more information).

In center-based modeling approaches, the cell-cell interaction force is usually written as a function of the location of each cell centers and almost always acts in the direction of the vector connecting the two interacting cells [76, 78, 83, 91]. Let $r_i$ be the position of the cell center of cell $i$, and let ${\mathbf{F}}_{\mathit{ij}}$ denote the force on cell $i$ due to cell $j$. If we assume that the force ${\mathbf{F}}_{\mathit{ij}}$ acts in the direction of the vector connecting the cells, then the force on cell $i$ due to cell $j$ is written as

$${\mathbf{F}}_{\mathit{ij}}=F_{\mathit{ij}}\frac{r_i-r_j}{\Vert r_i-r_j\Vert }$$E11

where $F_{\mathit{ij}}$ is the signed magnitude of ${\mathbf{F}}_{\mathit{ij}}$. The total force on cell $i$ is then

$${\mathbf{F}}_i=\sum _j{\mathbf{F}}_{\mathit{ij}}$$E12

where the sum is over all cells $j$ connected to cell $i$. As mentioned above, this force is usually taken to be balanced by a viscous drag as the cells move, so that the equation of motion for $r_i$, the position of the cell center of cell $i$, is:

$$\gamma \frac{{\mathit{dr}}_i}{\mathit{dt}}={\mathbf{F}}_i$$E13

where $\gamma$ is the viscosity coefficient which could, for example, represent adhesion between a cell and the underlying substrate. In simulations there are several choices for the numerical method used to update the position of each cell center, but the most simple is forward Euler discretization. Thus, the position of a cell center at time $t+\mathrm{\Delta }t$ is given by

$$r_i\left(t+\mathrm{\Delta }t\right)=r_i\left(t\right)+\frac{\mathrm{\Delta }t}{\gamma }{\mathbf{F}}_i.$$E14

The cell–cell interaction force, ${\mathbf{F}}_{\mathit{ij}}$, is usually a function of the overlap, ${\delta }_{\mathit{ij}}$, between two interacting cells which is given by

$${\delta }_{\mathit{ij}}=R_i+R_j-\Vert r_i-r_j\Vert$$E15

where $R_i$ and $R_j$ are the radii of cell $i$ and $j$, respectively, and $r_i$ and $r_j$ are the locations of the centers of cell $i$ and cell $j$, respectively. Different definitions of the contact area between cells are possible to further extend the model [76, 78, 83, 91]. There exists a wide range of force laws that have been used in the literature, ranging from simple linear laws to more complex nonlinear models that can incorporate nonhomogeneous properties such as cell-cell adhesion, cell-substrate adhesion, and cell bond breaking [76, 78, 83, 91]. For simplicity, we consider the linear law given by

$$F_{\mathit{ij}}=k_1{\delta }_{\mathit{ij}}$$E16

where $k_1$ is a stiffness parameter. Since the linear law implies there is a large attraction between distant cells, it is reasonable to define a cutoff range for interacting cells such that

$$F_{\mathit{ij}}=\left\{\begin{array}{lll}{k}_1{\delta }_{\mathit{ij}}& \text{for}& {\delta }_{\mathit{ij}}\ge {\delta }_{\mathit{min}}\\ 0& \text{for}& {\delta }_{\mathit{ij}}<{\delta }_{\mathit{min}}.\end{array}\right.$$E17

Some other choices for force laws include the Johnson-Kendall-Roberts (JKR) model [83, 91], the Hertz model [98, 99], and harmonic-like interaction [100, 101, 102]. The JKR model describes the interaction between two isotropic homogeneous spheres that are strongly adhesive, i.e., as soon as the two cells come into contact range, they immediately form a contact area of finite size by active deformation of the cell membrane. However, a unique property of the JKR model is that when two cells are pulled apart, they continue to interact for distances greater than the standard contact range but less than a specified distance (i.e., ${\delta }_{\mathit{min}}>{\delta }_{\mathit{ij}}>{\delta }_{\mathit{ij}}^c$). This phenomenon is referred to as hysteresis behavior. Similar to the JKR model, the Hertz model approximates each cell as a homogeneous sphere, but in the Hertz model, hysteresis behavior is not included. In the Hertz model, the potential interaction is represented as the sum of a repulsive and an attractive force (i.e., $F_{\mathit{ij}}=F_{\mathit{ij}}^a\mathit{ttr}+F_{\mathit{ij}}^r\mathit{ep}$). On the other hand, harmonic-like force laws provide a simple model for the short-range cell-cell interaction due to cell-cell adhesion and elastic deformation that approximate two adhesively interacting cells by cuboidal objects with a non-deformable core, linked by linear springs. Each of these definitions for cell-cell interaction force offers their own set of advantages and disadvantages.

3.2 Application of center-based models in biology

Center-based models have been used to analyze multicellular processes in tumors [91, 99, 100], intestinal crypts and epithelial tissues [85, 103], cell migration in extracellular matrix [85, 100], and tissue regeneration, growth, and organization [91, 99, 100, 103, 104]. In addition, several center-based models have been developed for studying macroscopic properties of yeast colonies [92, 93, 105]. (See Section 4 for more details.) Although the mechanical information that can be extracted from these models is rather approximate, various biophysical aspects in these problems have explained important mechanisms for proper cell organization and speed of growth in tissues [78, 97, 106].

3.3 Vertex models

Two-dimensional vertex models are a well-known class of lattice-free, cell-based models that provide a more detailed description of cell shape [84, 87, 107, 108]. In vertex models, each cell membrane is represented as a polygon composed of vertices and edges that are shared between adjacent cells (Figure 5A). A set of rules or an equation of motion defines how each vertex moves, and collective movement of vertices leads to changes in cell shape over time (Figure 5A). Most often, the equation of motion is a force or potential function based on the current configuration of vertices, i.e., location and connection between pairs of vertices, as well as other geometrical features such as area or perimeter of cells [84, 87, 107, 108]. Many vertex models also incorporate rules to govern changes in connections between vertices to allow for rearrangements in cell neighbor relationships [107]. The precise equations of motion used differ between models and may be adapted to suit each particular biological problem. Below we review some of the common forms of the equation of motion and typical rules for rearrangements of vertex connectivity.

Figure 5.

Vertex model. Vertex models depict each cell membrane as a polygon of vertices and edges that are shared between adjacent cells. At each time step, an equation of motion determines how each vertex moves, and collective movement of vertices over time results in changes in cell shape. (A) Example forces acting on vertex $i$ due to interactions with shared vertices in cell $\alpha$ according to Eq. (18). (B–D) Cells in the model form and break bonds due to neighbor rearrangements. Here we describe three basic cell neighbor rearrangements accounted for by most vertex models [87, 107]. (B) T1 transitions occur when two vertices sharing a short edge (defined by a minimum distance much smaller than the average length of an edge) merge into a single vertex which is then reassigned to two new vertices, thus changing the local network topology. (C) T2 transitions occur when a given cell shrinks to an approximately zero area and is removed to represent apoptosis. (D) T3 transitions occur when the intersection of a vertex with an edge is avoided by replacing the approaching vertex with two new vertices. For more detail see Section 3.3.

In most cases, the equation of motion for vertices is deterministic since a reasonable assumption for tightly packed cell aggregates is that stochastic motion of cells is mitigated by strong interactions between cells [84, 87, 107, 108]. One difference among vertex models in the literature lies in the definition of the force. The choice of form for the function ${\mathbf{F}}_i$ reflects which forces are thought to dominate cell mechanics for the system being studied. Some commonly modeled forces include tension or elastic forces, due to the combined action of a cell’s actomyosin cortex [107] and adherens junctions [107], or pressure, due to hydrostatic pressure [87]. Note that the forces acting on each vertex may either be given explicitly, or else an energy function may be specified, whose gradient is assumed to exert a force on each vertex.

For a concrete example, we will consider an equation of motion built from an energy function given by Farhadifar et al. which has recently seen extensive use in modeling wing disk epithelia [87]. This energy function encodes constraints associated with the limited ability of cells to undergo elastic deformations, volume changes, and other movements due to adhesion to other cells. The energy function is defined by

$$E=\sum _{\alpha =1}^N\left(\frac{K_{\alpha }}{2}{\left(A_{\alpha }-A_{\alpha }^{\left(0\right)}\right)}^2+\frac{{\mathrm{\Gamma }}_{\alpha }}{2}{P}_{\alpha }^2\right)+\sum _{i,j}{\mathrm{\Lambda }}_{\mathit{ij}}{l}_{\mathit{ij}}$$E18

where the area elastic modulus $K_{\alpha }$, bond tension parameter ${\mathrm{\Lambda }}_{\mathit{ij}}$, perimeter coefficient ${\mathrm{\Gamma }}_{\alpha }$, and preferred cell area $A_{\alpha }^{\left(0\right)}$ are all model parameters whose values must be chosen in a biologically relevant manner.

In vertex models that use this energy function, the choice of preferred cell area $A_{\alpha }^{\left(0\right)}$ depends on how cell growth is incorporated into the model. For a given energy function $U$, the force on each vertex is determined by the negative derivative of the energy with respect to the coordinates of that vertex $F_i=-{\nabla }_{i}U$, where ${\nabla }_i$ denotes the gradient operator evaluated at ${\mathbf{x}}_i$. Computing the gradient of Eq. (18) and exploiting the fact that the movement of vertex $i$ affects only the energy of the cells associated with it, the force $F_i$ may be written explicitly as

$$F_i=\sum _{\alpha \in S_i}\left(-K_{\alpha }\left(A_{\alpha }-A_{\alpha }^{\left(0\right)}\right){\nabla }_i{A}_{\alpha }+{\mathrm{\Gamma }}_{\alpha }{P}_{\alpha }{\nabla }_i{P}_{\alpha }\right)-\sum _j{\mathrm{\Lambda }}_{\mathit{ij}}{\nabla }_i{\mathrm{ł}}_{\mathit{ij}}$$E19

where now the first sum runs over cells associated with vertex $i$ and the second sum runs over vertices sharing an edge with it.

Many of the earliest vertex models simulated cross-sections of epithelial tissues and so did not account for cell neighbor rearrangements [87]. However, later models were formulated to simulate top-down dynamics of epithelial cells in a plane which required modeling cell neighbor rearrangements [107]. Namely, to accurately describe planar epithelial cell dynamics, cells in the model must be allowed to form and break bonds by changing connectivity among vertices as described in Figure 5.

3.4 Application of vertex models in biology

Vertex models were initially developed to study the packing of bubbles in foams [109]. Similar to foams, cells in epithelial tissues are tightly packed and mechanically coupled to their neighbors by adhesion molecules along their common interfaces [84, 87, 107, 108]. In addition, cells exert forces onto each other and their environment. Since vertex models are well suited for modeling tightly packed cell ensembles with negligible intercellular space, the vertex modeling framework was later adapted to study two-dimensional packing and rearrangement of apical cell surfaces in planar epithelia [84, 87, 107, 108]. In these studies, epithelial cells are described by a planar two-dimensional network of vertices that defines the apical cell surfaces as polygons with straight interfaces between neighboring cells. The underlying assumption of two-dimensional apical vertex models is that the main forces acting to deform the cells are generated along the apical cell surfaces or can be effectively absorbed in the apical representation. The work of Odell et al. [110] first demonstrated the use of vertex models to study how spatial patterning in either active forces or passive mechanical properties may lead to tissue deformation. Subsequent studies by several groups have examined a variety of alternative patterns of force and material properties that can give rise to similar tissue deformations [84, 87, 88, 107, 108, 111].

An important behavior of biological cells that does not affect the dynamics of foams but plays a large role in the structure of biological tissues is the ability of cells to grow, divide, and die. Model components representing cell division and death require the modification of tissue connectivity which can add significant computational complexity, but these cell behaviors are essential mechanisms of many of the biological process being studied. Thus, one biological question that has been addressed using vertex models extended to include cell death and division is the control of packing geometries in an epithelial sheet [84, 87, 88, 107, 108, 111]. A highly cited example is the work of Farhadifar et al. [87] who performed a systematic analysis of the equilibrium cell packing geometries and their dependence on cell mechanical and proliferative parameters using the Drosophila wing disk as a model system. By comparing simulations with experimental results, the authors arrived at a set of parameter values for which their model accounts for vertex movements based on biologically relevant cell area variations, cell rearrangements due to laser ablation, and epithelial packing geometries seen in vivo. This work demonstrates how vertex models may be parametrized and tested using experimental data.

In many biological systems, patterns of mechanical stress and resulting macroscopic tissue shape changes may happen concurrently and can feed-back into each other. Vertex models can be easily modified to incorporate both mechanical and chemical feedbacks. For example, Odell et al. [110] included a simplified feedback in their earliest model by assuming that contraction of individual cells was activated by high levels of stretching due to tissue deformation. Although vertex models have been successful in testing many of the mechanisms that govern the macroscopic behavior of tightly packed epithelial tissues, this class of model generally ignores contributions from important components such as cell-matrix interactions [112], actomyosin contractility and other biophysical properties of the cell membrane [113], and active remodeling of cytoskeletal components. Thus, there have been several notable extensions of the vertex model to address these important factors [110, 114, 115]. Although vertex models provide a promising modeling framework for further investigation of different biological systems, in the next section, we review a different modeling framework that was developed to take into account more detailed representations of cell deformation and provide an extension to the types and complexity of mechanical stress patterning and feedback considered.

3.5 Subcellular element model

The subcellular element (SCE) model provides a framework for simulating lattice-free multicellular structures in which the shape of each cell dynamically emerges from interactions with the local environment due to model assumptions about the underlying mechanical properties of the system [82, 86, 89, 90]. In the model, each cell is composed of a large, and possibly varying, number of small nodes called subcellular elements (Figure 6). Each subcellular element of a cell is modeled as a single point at its center of mass, which changes position over time subject to three processes: (i) weak random fluctuations, (ii) elastic interaction with elements of the same cell, and (iii) elastic interaction with elements of other cells.

Figure 6.

Subcellular element model (SCE). The SCE model uses sets of elements/nodes to represent individual cells. At each time step, an equation of motion is used to update the location of each node based on resulting forces from interactions between nodes due to cell growth, deformation, etc. (A) Illustration of a SCE model where membrane and cytoplasm of individual cells are represented by one set of homogeneous elements/nodes. Elements/nodes of the same cell are held together via intracellular interactions (solid lines) modeled using short-range potential functions ($\nabla {\mathit{E}}_{\mathit{\text{intra}}}$). Elements/nodes of two neighboring cells interact via intercellular interactions (dashed line) modeled using potential functions that weakly bind adjacent elements/nodes ($\nabla {\mathit{E}}_{\mathit{\text{inter}}}$). The force applied to node ${\mathit{\alpha }}_{\mathit{i}}$ is calculated based on interactions with all other nodes ${\mathit{\beta }}_{\mathit{i}}\in$cell_i and adjacent nodes ${\mathit{\alpha }}_{\mathit{j}}$ in neighboring cell, cell_j. (B) SCE model where the membrane (filled nodes) and cytoplasm (clear nodes) of one cell are represented separately using two different sets of elements/nodes. The force applied to cytoplasm node ${\mathit{\alpha }}_{\mathit{i}}$ is calculated based on intracellular interactions which include interactions with all other cytoplasm nodes in the same cell (${\mathit{E}}^{\mathit{II}}$) and interactions with all other membrane nodes of the same cells (${\mathit{E}}^{\mathit{MI}}$) via short-range potential functions. The force applied to membrane node ${\mathit{\beta }}_{\mathit{i}}$ is calculated based on both intracellular interactions and intercellular interactions. Intracellular interactions for membrane node ${\mathit{\beta }}_{\mathit{j}}$ include interactions with all cytoplasm nodes of the same cell (${\mathit{E}}^{\mathit{MI}}$) and interactions with adjacent membrane nodes of the same cell (${\mathit{E}}^{\mathit{MMS}}$), and intercellular interactions for membrane node ${\mathit{\beta }}_{\mathit{j}}$ include interactions with adjacent membrane nodes from neighboring cells. For more details about these interactions, refer to Section 3.5.

In the SCE modeling approach, the membrane and cytoplasm of each cell can be represented together as one set of homogeneous elements/nodes or separately using two different sets of elements/nodes (Figure 6). Collective interactions between pairs of nodes from the same cell represent bulk cytoplasmic contents, and collective interactions between pairs of adjacent nodes from neighboring cells represent volume exclusion of cells as well as adhesive properties. In models that consider the membrane and cytoplasm separately, collective interactions between pairs of cytoplasm nodes and membrane nodes represent the cytoplasmic pressure of the cell applied to the cell membrane. Biomechanical and adhesive properties of cells are modeled through viscoelastic interactions between elements represented by phenomenological potential functions such as linear springs or Morse potential functions, which are used to simulate close-range repulsion (modeling volume exclusion of neighboring segments of cytoskeleton) and medium-range attraction between elements of the same or different cells (modeling the adhesive forces between segments of cytoskeleton) [116].

The potential functions described above are used in model equations to calculate the displacement of each element/node at each time step based on their interactions with neighboring elements/nodes resulting in the deformation of cells within the tissue (Figure 6). As mentioned previously, nodes are assumed to be in an overdamped regime so that inertia forces acting on the nodes are neglected [83, 87, 88, 89, 90]. This leads to the following two equations of motion describing the movement of element/node ${\alpha }_i$ in cell $i$:

$$F_{{\alpha }_i}=-\left(\sum _{{\alpha }_i\ne {\beta }_i}\nabla E_{\text{intra}}-\sum _{i\ne j}\sum _{{\alpha }_i}\nabla E_{\text{inter}}\right)$$E20

where $\eta$ is the damping coefficient, $x_i$ is the position of node $i$, $N$ is the total number of nodes, and $\nabla E$ is the potential function used to describe the interaction of node $i$ with node $j$. This equation is discretized in time using the forward Euler method, and the position of node $x_i$ is incremented at discrete times as follows:

$$x_{{\alpha }_i}\left(t+\mathrm{\Delta }t\right)=x_{{\alpha }_i}\left(t\right)-F_{{\alpha }_i}\frac{\mathrm{\Delta }t}{\eta }$$E21

where $\mathrm{\Delta }t$ is the time step size.

One of the important features of the SCE modeling approach is the ability to change parameters of the potential functions that are used to describe interactions between elements to calibrate model representations of biomechanical properties of a particular type of a cell directly using experimental data. More specifically, the SCE model can be used to perform in silico bulk rheology experiments on a single cell in order to scale the parameters such that the passive biomechanical properties of each cell are independent of the number of elements used to represent each cell [117]. As a result, SCE simulation output captures the underlying biomechanical properties of the real biological system being studied and remains relevant regardless of the choice of the number of elements used in the model.

As indicated in Fletcher et al. [84], computational experiments follow a creep-stress protocol in which a constant extensile force is applied to the end of an SCE cell whose opposite end is fixed. Before the extensile force is released, the strain is measured as the extension of the cell in the direction of the force relative to its initial linear size. In silico estimates of the viscoelastic properties of cells modeled using the SCE approach have been shown in many biological applications to agree with in vitro rheology measurements [117, 118]. This indicates that the simple phenomenological dynamics of the SCE modeling approach are enough to capture low to intermediate responses of cytoskeletal networks over short timescales ($\sim$10s) [118]. Over longer timescales ($\sim$100 s), cells respond actively to external stresses by undergoing cytoskeletal remodeling, and this phenomenon can be incorporated into the SCE modeling approach by inserting and removing subcellular elements of a cell in regions under high or low stress [86].

The generalized Morse potential functions implemented in SCE modeling approach are commonly used in physics and chemistry to model intermolecular interactions [119] and in biology to represent volume exclusion of neighboring regions of the cytoskeleton [94, 95, 120, 121, 122, 123, 124, 125]. While it is difficult to associate specific potential functions directly with specific cytoskeletal components of cells, computational studies of bulk properties at the tissue level have suggested that the precise functional form of the potential used in modeling has a small impact on overall system dynamics [97, 117].

3.6 Application of SCE models in biology

The subcellular element (SCE) modeling approach has been successful in modeling mechanical properties of individual cells as well as their components and determining individual cell impact on the emerging properties of growing multicellular tissue as well as describing cellular interactions with mediums such as the extracellular matrix and fluids [76, 84, 86, 89, 90, 95, 96, 111, 117, 120, 122, 124, 125, 126, 127]. The general modeling framework was initially developed by Newman et al. [89] for simulating the detailed dynamics of cell shapes as an emergent response to mechanical stimuli. Recent applications of the SCE modeling approach show that it is flexible enough to model additional diverse biological processes such as intracellular signaling [122], cell differentiation [94], and motion of cells in fluid [124].

To date, the SCE has not been widely used to study biological processes outside the area of epithelial morphogenesis. Christley et al. [122] developed a model of epidermal growth on a basal membrane that incorporates cell growth through the addition on new elements and division by redistributing a cell’s group of elements between two new daughter cells. This mode was a novel implementation because it was coupled to a subcellular gene network representing intercellular Notch signaling. The SCE modeling framework has also been coupled to a fluid flow model to simulate the attachment of platelets to blood vessel walls [124]. Using the SCE framework for modeling individual platelets in simulations provides a detailed mechanical model of individual platelet behavior that allowed the authors to examine the relationship between platelet stiffness and movement in the fluid. This facilitated the generation of new hypothesis about mechanisms platelets use to adhere to injured sites on the blood vessel walls. In the context of developmental biology, the primary application of the SCE to date has been a computational study of primitive streak formation by Sandersius et al. [86]. In addition, the SCE model has been applied to model development in both animal [95] and plant [94] systems.

4. Conclusion

In this chapter, we described the medical and scientific importance of studying misfolded protein diseases as well as provided a broad description of several classes of mathematical models that can be used to further investigate the underlying mechanisms governing protein aggregation and propagation in a multicellular system. In Section 2, we gave an overview of different modeling frameworks that have been developed and validated for studying protein aggregation in both yeast and mammalian systems. In Section 3, we described in detail various cell-based modeling frameworks that have been successfully used to gain insight about the impact of individual cell behaviors on macroscale properties of tumors, developmental tissues, and yeast colonies. However, our main goal in this chapter was to familiarize the reader with each class of model in order to facilitate further discussion about how these two types of models could be combined to develop a more complete representation of prion disease dynamics within an actively growing and dividing yeast colony.

While current models of protein aggregation in yeast have successfully provided further insight into important mechanisms driving prion disease dynamics (i.e., conversion and fragmentation), they have been unable to recreate a number of important physiological characteristics including variable phenotype induction rates that result in sectored phenotypes among a single yeast colony (see Section 1 for details). One reason for this may be that a large number of models representing intracellular dynamics of protein misfolding diseases were developed for studying protein aggregation dynamics in isolation and these models disregard the contribution of individual cell behaviors within the growing yeast colony as a possible mechanism governing prion disease dynamics. A major open question in prion biology is to understand how prion aggregates spread between cells within a whole colony or tissue. Experimental observations such as sectoring provide compelling data that transmission mechanisms other than what is addressed by current aggregation-only-models must play a role in the presence and persistence of prion disease phenotypes in yeast colonies, i.e., processes such as conversion, fragmentation, nucleation, and even enzyme-mediated fragmentation alone cannot entirely explain the spread of prion disease throughout a yeast colony. Thus, in order to test hypotheses about the impact of individual cell behaviors on the spread of misfolded proteins, it is necessary to develop a novel modeling framework.

Developing a modeling framework for investigating prion disease dynamics within an entire yeast colony is challenging because it requires capturing the physical processes on an individual cell level that determine colony growth (i.e., budding and variable cell cycle length) as well as capturing the interplay of individual cell processes with protein aggregation dynamics (i.e., asymmetric protein distribution at the time of division, persistence of diseased phenotype/aggregate that was given to daughter while it grows to begin a new cell cycle, lineage-dependent protein propagation). Several models have already been developed to investigate physical mechanisms controlling patterns of yeast colony growth such as cell division polarity, mother-daughter size asymmetry, and cell-cell adhesion via budding of the new daughter cell [92, 93, 105]. Jönsson and Levchenko developed an off-lattice, center-based model in which cells are modeled as elastic spheres of variable size [92]. Their work showed that cell growth inhibition by neighboring cells and polar division growth patterns were the most significant factors driving appreciable differences in the shape and size of yeast colony development. In addition, Wang et al. built an off-lattice, center-based model that incorporates key biological processes in yeast cell colonies including budding, mating, mating type switch, changes in cell cycle length and cell size due to aging, and cell death [93]. The main results of their work include proposed mechanisms for how budding patterns in yeast cells affect colony growth. Both of these studies serve as excellent starting points for developing a modeling framework that includes the physical processes and individual cell behaviors of yeast cells that govern colony growth.

One way forward would be to extend current cell-based models for yeast colony growth by combining them with a detailed model of protein aggregation dynamics as well as implementing specific model components that capture the interplay of individual cell behaviors with protein aggregation dynamics. Specifically, the new modeling framework would need to include asymmetric protein distribution at the time of division and varied cell cycle lengths to represent how growth of a new daughter cell impacts persistence of aggregates that were inherited at time of division and track protein concentration through different cell lineages to investigate the impact of lineage on colony phenotypes in yeast. In the model, prion dynamics would need to be simulated within each individual cell using models of intracellular dynamics described in Section 2. In addition, using a center-based model would allow for the spatial arrangement of cells to account for the effect of biophysical properties such as increased adhesion between a mother and daughter cell during budding.

Preliminary simulations (Figure 7) show that individual cell behaviors indeed impact yeast colony structure as well as protein aggregation dynamics in a growing yeast colony. For our preliminary framework, we extended current cell-based models for yeast colony growth by combining them with a simplified model of protein aggregation dynamics as well as implementing asymmetric protein distribution at the time of division and strong adhesion between mother and daughter cells during budding. These are two cell behaviors we hypothesize to play an important role in protein aggregation dynamics in yeast colonies. Our simulations begin with a single founding cell that divides to give birth to the first daughter. All cells continue to grow and divide throughout the simulation and stop once the colony reaches about 1500 cells total. Successive daughter cells of the founder cell are enumerated in the order of their division. For example, the first daughter is the cell born after the first division of the founder cell, the second daughter is the cell born after the second division of the founder cell, and so on.

Figure 7.

Impact of individual cell behaviors on yeast colony structure and protein aggregation. In our preliminary simulations, we extended current cell-based modeling approaches for yeast colony growth by combining them with a simplified model of protein aggregation dynamics along with asymmetric protein distribution at cell division and strong adhesion between mother and daughter cells during budding. Simulations begin with a single founder cell and stop after 15 hours when the colony contains about 1500 cells total. (A) Example of sub-colony structure for each of nine daughter cells born from the original founder (black cell). Each sub-colony is depicted with a different color indicated below. (B) We show in white the fifth-generation sub-colony corresponding to the fifth daughter of the founder cell and all of its offspring. (C) Example of the interplay between protein aggregation dynamics and individual cell behaviors in a simulation leading to a sectored colony. Protein aggregation dynamics for individual cells were modeled using a logistic growth function with a strong Allee effect which results in stable extinction of protein aggregates in the sub-colony of an early daughter cell. The number of aggregates per cell is represented by different colors given in the key below.

Figure 7A shows sub-colony structure for each of nine daughters cells born from the original founder cell of the colony represented in black. The first daughter sub-colony (grayish black) occupies a larger area than any of the other sub-colonies. In addition, the average location of each successive generation is located closer to the boundary (Figure 7B). Protein aggregation dynamics were modeled using a logistic growth function with a strong Allee effect. Figure 7C shows the number of aggregates per cell in a simulation that resulted in a sectored colony. The interplay between protein aggregation within each cell and individual cell behaviors results in complete loss of aggregates in an early daughter cell, and subsequently, the sub-colony generated by that daughter remains free of aggregates throughout the simulation. These results suggest that the unified model may have the potential to predict mechanisms underlying experimentally observed phenomena such as sectored prion phenotypes in yeast colonies in addition to serving as a tool for future hypothesis generation and testing.

Although our preliminary studies appear promising, there is still a large gap to fill in terms of developing models that consider underlying microscopic processes of protein aggregation together with macroscopic properties of the environment in which they are taking place. Indeed, as mentioned in Section 2 there have been recent efforts to model the impact of tissue structures on the spread of protein misfolding diseases in mammalian systems [63, 64]. However, most of the large-scale models for protein aggregation lack details of components at the microscopic scale that would allow for the study of the interplay of different temporal and spatial scales. We hope that future modeling studies will begin to incorporate both scales in order to test hypotheses about mechanisms that can explain unresolved experimental data and yield new strategies for treating protein misfolding diseases.