article: We introduce and analyze a minimal model of epigenetic silencing in budding yeast, built upon known biomolecular interactions in the system. Doing so, we identify the epigenetic marks essential for the bistability of epigenetic states. The model explicitly incorporates two key chromatin marks, namely H4K16 acetylation and H3K79 methylation, and explores whether the presence of multiple marks lead to a qualitatively different systems behavior. We find that having both modifications is important for the robustness of epigenetic silencing. Besides the silenced and transcriptionally active fate of chromatin, our model leads to a novel state with bivalent (i. e. , both active and silencing) marks under certain perturbations (knock-out mutations, inhibition or enhancement of enzymatic activity). The bivalent state appears under several perturbations and is shown to result in patchy silencing. We also show that the titration effect, owing to a limited supply of silencing proteins, can result in counter-intuitive responses. The design principles of the silencing system is systematically investigated and disparate experimental observations are assessed within a single theoretical framework. Specifically, we discuss the behavior of Sir protein recruitment, spreading and stability of silenced regions in commonly-studied mutants (e. g. , sas2, dot1) illuminating the controversial role of Dot1 in the systems biology of yeast silencing.
  To provide the relevant context of our approach to a wider audience, we briefly review the molecular biology of yeast silencing system and its central puzzles. In budding yeast, the Silent mating-type Information Regulator (SIR) proteins have been identified to play a pivotal role in heterochromatin formation [30], [37]–[39]. The Sir complex includes four proteins, of which Sir1 facilitates assembly of Sir2–4 at the Hidden Mating (HML/HMR) Cassettes, in concert with proteins factors like Orc, Rap1 an Abf1, all of which have specific binding sites at the silencer of HML/HMR. In telomeres, the players are slightly different; Sir1 is not needed for the nucleation of Sir complex [38], [39] and yKu DNA-end binding complexes is important. One of the proteins of the Sir complex, Sir2, is a NAD-dependent deacetylase [40]. The acetylation of Histone tails antagonizes silencing [21], [41]–[44]. The deacetylase activity of Sir2 is critical for the spread of silencing by presumably increasing Sir3-Sir4 affinity for nucleosomes [45]. Specifically, acetylation of H4K16 by Sas2 (Something About Silencing) activity undermined the spread of silencing [6], [7], [46], [47]. Recent studies have established that besides the interaction of Sir proteins with nucleosomes [48], interactions between Sir proteins themselves are crucial in the spread of silencing [18], [45], [49]. This spreading in telomeres appears to be unhindered by silencer boundary element and is arguably stochastic [33], . In contrast, in the HML/HMR loci boundary elements prevent the spread of Sir complex [7], [53]. Sir spreading happens in only a few of the natural telomeres [15], [33], [54]. Moreover, mechanisms like the clustering of telomeres anchored to the nuclear envelope is perhaps important in subtelomeric silencing [52]. Sir2 does not exclusively deacetylate H4K16 residue. Other residues like H3K9 and H3K14 are also deacetylated moderately [34], and H4K56 extensively by Sir2 [55]. However, the spreading of Sir2 compete primarily with H4K16Ac [46]. The compaction into silent heterochromatin requires further deacetylation of H4K56Ac [55]. The effect of spread of Sir proteins on silencing is subtle—spreading is necessary but not sufficient for heterochromatin formation and Sir binding is broader in scope [15]. Though the compaction of chromatin owing to stable Sir-complex-binding blocks association of RNA polymerase II thereby preventing transcriptional activity [7], Sir binding and spreading itself may be a very dynamic phenomenon. Not surprisingly, Sir protein binds at various loci which are not necessarily silenced, as deduced from Chromatin Immunoprecipitation (ChIP) studies [15], [56]. The association of Sir proteins does not immediately lead to gene silencing [57] arguing for added steps to the formation of heterochromatin [9], [10], [12], [58]. Nevertheless, Sir protein association is a prerequisite in establishing an inheritable pattern of silencing in the telomeric and HML/HMR loci. Active sites in budding yeast chromatin are not only hyper-acetylated at particular histone tail sites but also tend to be hyper-methylated at certain key residues [23], [32], [35], [51], [59]. One of the DOT (Disruptor Of Telomeric Silencing) proteins, Dot1, methylates H3K79 residue and competes in binding with Sir3 for the same basic patch on the histone core region (H3K79) [45], [59]–[61]. Dot1 can mono- di- or tri- methylate this residue distributively [62]. Studies have revealed that H4K16Ac displaces Sir3 binding, thereby aiding Dot1-mediated methylation of H3K79 [60]. What other histone marks (or other factors) regulate transcriptional activity is a complex question [63]. The enzymatic activity of Dot1 itself is modulated by various other factors. For example, Paf1 complex (RNA Polymerase Associating Factor) which is known to be an important factor in transcription elongation plays a crucial role in Dot1 methylase activity [64]–[66]. In fact, such positive feedback from transcriptional elongation into establishing transcriptionally-active marks is not uncommon. For example, H3K4 methylation [67]–[69], which is an active mark in yeast, requires Paf1 [66] for being methylated by Set1 [70]. The phenomenon of yeast silencing has several other components many of which are poorly understood. Though the structure and interaction domains of Sir proteins are well-studied [45], [48], [49], [71]–[73], protein factors, histone modifications etc. which interact with the system continue to be uncovered. For example, Rad6 dependent ubiquitylation of H2B-K123 influences methylation by Dot1 and Set1 [69], [74]. Ubiquitylation is a histone mark implicated in transcriptional initiation and elongation [75], and such a trans-histone regulatory pathway acts as a feedback into the establishment of transcriptional activity. Another such feedback in establishment of active marks is the interaction of Dot1 with a histone acetyltransferase [76]. Sumoylation is a silencing histone mark, and might play an important role in the silencing of subtelomeric regions [77], [78]. Variants of the histones [79], [80], particularly H2A. Z, plays a complex role in transcriptional activity, and has been shown to deter the ectopic spread of silencing [68]. The topic of histone variants invite a host of connections between chromatin assembly and transcriptional activity. These influences are perhaps peripheral to the central mechanism resolving epigenetic fates [81]–[83]. We argue that the epigenetic fates are reinforced both from transcriptional and silencing feedbacks—our focus is the core mechanism, summarized in Fig. 1. Recent experiments [15] have shown that binding of Sir3 protein is not indicative of heterochromatin and many genes in euchromatin, including highly transcribed genes, may show wide-spread Sir3 binding. Such a co-occurrence of silencing and active marks called bivalent chromatin states has been discussed in the context of Polycomb silencing [84] and that of HP1 binding in expressed exons [85]. We show how stable bivalent states emerge naturally on modelling the known interactions. In the next section we define the mathematical model based on this minimal picture.
  In this section we introduce the key components of our model in order to facilitate discussion, further details are relegated to Materials and Methods and Text S1. Continuing on the spirit of earlier work involving one of the authors [24], we define several local states that represent the various modifications of nucleosomes. Five possible nucleosome-states are considered at nucleosome: , , , and. The state is Sir-complex-bound (without making a distinction of the various Sir proteins). The state is acetylated histone (H4K16Ac). Multiple histone tails on the same nucleosome can get acetylated—we model the average level of acetylation. The is the methylated state (H3K79Me), and we treat the multiple levels of methylations on the average. The is unmodified state— it is neither acetylated, methylated or Sir-protein bound. The state is the transcriptionally active state and has multiple histone marks, particularly both methylation and acetylation marks. The model considers the following biochemical processes, see Fig. 2 for a pictorial representation of the model. Index are for nucleosomes within a local neighborhood of nucleosome: The mathematical details of the model are relegated to Methods and Materials and Text S1. Therein we establish that the key parameters are the rates of Sas2 activity, the cooperative Dot1 activity, the cooperative Sir2 activity and the cooperative Sir binding. The first result of our model is the phase space—a four dimensional parameter space with each point, being a distinct choice of parameters, produces solutions for the density of marks in the model. The distinct properties of such solutions, like stability and density of marks, define distinct phases. For example, regions of bistability are where stable active () and silent () states can coexist. Stability in this context is against the loss of epigenetic marks through various perturbations—histone turnover, DNA replication, stochasticity in enzymatic reactions, deacetylase/demethylase activity etc. and maintaining heritable distinct epigenetic fates for the chromatin loci. Bivalent region is where the stable solution exhibits simultaneously high local density of active and silenced marks. The ranges of values of the parameters we introduce above are unknown. Various loci in the wild type cell corresponds to different combinations of parameters and therefore maps to distinct points in phase space. These wild-type points should lie in regions where perturbations are unlikely to affect the epigenetic fates of the loci—a requirement of robust engineering design. We have studied aspects of this robustness elsewhere [25], [86], wherein we have discussed the requirements on engineering design for faithful reestablishment of epigenetic states from redistribution and dilution of modified histones during mitosis. A quantitative measure of robustness of a epigenetic states is the volume and shape of parameter space where the state is supported [87]. We first outline the phase space. The model has four types of uniform steady-state solutions (in mean-field analysis, see Materials and Methods and Text S1): the silenced state which has probability weight predominantly in the state, the active state with weight predominantly in the state, the intermediate state with weight distributed in, and, and the bivalent state with weight distributed predominantly over and. The model exhibits bistability between the silenced state and the intermediate state, and the silenced state and the active state. The bivalent state is monostable and is not to be confused with the bistable regions. A novel outcomes of our analysis is that only four distinct stable epigenetic fates emerge in all of phase space. Visualizing the four dimensional phase space is challenging—we present different two dimensional sections of the phase space for the purpose. Each section is defined by fixing the values of any two of the four parameters. The section shown in Fig. 3, for which and were fixed, is a representational one for the model' s phase space. Only stable solutions are shown. The region of stable silenced state is depicted by red diamonds, the stable active with green stars, stable bivalent with magenta circles, and stable intermediate with blue crosses in our phase plots. We adhere to this convention throughout the paper. The supply of Sir proteins is limited in order to prevent ectopic silencing in wild type. The reaction rates and, for basal and cooperative Sir binding respectively, are proportional to the concentration of ambient Sir proteins. In telomeres where no obvious boundary elements are present to limit the spread of silencing proteins, the concentration of ambient Sir proteins self-adjusts such that the spreading of silencing is stochastically stationary; the choice of parameters for which such stationarity is achieved is known as the zero velocity line [24]. Away from this line, but still within a bistable region, either the active state or the silenced state spreads. Denoting by the the total number of Sir proteins, and the density of Sir proteins, the equation that determines is (1) where is the volume of the cell nucleus. As silencing spreads, increases, and decreases. Such a self-adjusting is the titration effect on Sir spreading. Sir spreading needs to happen only in a few of all the telomeres [15] in order for this effect to be observed and be important for perturbations. The zero-velocity line determines how the local Sir recruitment rate evolves under perturbations, as ellaborated in the next sections. The epigenetically stable bivalent state is a novel outcome of our model, and we discuss the implications of its emergence under the inhibition of Sir2 deacetylase activity, i. e. , on reducing in our model. Because of the titration effect, whenever a single parameter is changed the parameter self-adjusts to maintain the system on the zero-velocity line. In effect, the silencing front progresses or retreats to be stochastically stationary again. In Fig. 4, we show the cross section of the phase space for fixed Dot1 activity () and Sas2 activity (). The zero velocity line is determined in stochastic simulations, as described in Materials and Methods, and is plotted over the mean-field phase diagram. As we decrease Sir2 deacetylase activity (), traces the zero-velocity line driving the system eventually out of the bistable region (wild-type) and into the region of monostable bivalent state, see Fig. 4. The implications of this observation are two fold— The first observation is a non-trivial prediction of the model. Though there has been reports of defective silencing boundary [46], [47], [88], [89] in telomeres for Sir2 and Sas2 perturbations, the nature of the defect for Sir2 inhibition has not been made precise. The second observation implies that effective Sir cooperativity () increases locally for decreasing deacetylase activity () of Sir2, and is relevant to the spatial feature of the bivalent state, which we investigate in lattice simulations. We study the spreading and steady-state occupancy of silencing marks with a nucleating center for Sir binding at one lattice end. We observe that the bivalent state is stochastically established Sir occupancy. The silencing boundary is ill-defined as shown in Fig. 5. A precise characterization of the typical size of these patches is the correlation length of Sir occupancy—a measure of how influential the state of a nucleosome is on the state of neighboring ones, see Fig. 6. We measure the correlation length of marks for the system self-tuning along the zero-velocity line shown in Fig. 4. The correlation length is high, as expected, in the bistable region where silencing domains are established. However, the bivalent state can maintain rather high correlation lengths resulting in local patches of silencing domains; see further discussion in the next section. The bivalent state is truly bivalent, in the sense that nearby nucleosomes carry opposite marks, only in selected regions of the parameter space. The precise role of Dot1 is the yeast silencing continues to be active debated [13], [36], [54], [69], [76], and a commonly used assay to report the telomeric position effect variegation (TPEV) has been brought to question for dot1 strain [90]. Experiments [23], [35], [54], [59], [91] have focused on both dot1 and overexpression of Dot1, but not inhibition. In this section, we summarize the bearings of our model on Dot1 perturbations. The dot1 strain is bistable in our model, but the bistability is for intermediate and silenced states, as seen in Fig. 7 to be the region of overlap of silenced-stable (red diamonds) and intermediate-stable (blue crosses). The bistability achieved in the absence of Dot1 activity may cast doubt on the necessity of Dot1 in our minimal model construction. We argue that Dot1 plays critical role in establishment of heritable silenced domains. We argue that measuring Sir occupancy, as opposed to transcriptional activity resolves the pivotal role played by Dot1. This is because transcriptional activity does not imply a unique state of local histone modifications. Though hyper-acetylation and hyper-methylation is associated with transcriptional activation in genome-wide studies [92], methylation is not required for moderate transcriptional activity [9], [69], [93]. A subtlety of our model, and a potential criticism, is our definition of an active state, viz. , states having both acetylation and methylation marks, as being too stringent. However, we can always consider domains to be moderately transcribed. We observe in Fig. 7 that there is a critical value of, which is given by (see Text S1), above which the system is active-silenced bistable and below which the the system is intermediate-silenced bistable. In the presence of non-zero basal rates of Dot1 activity () and Sir binding () the sharp transition between the two bistable regions becomes a crossover (see Text S1). Therefore, if we use transcriptional activity as a measurement, the model implies that increasing Dot1 level (β) from null simply makes transcription robust. Is Dot1 redundant in epigenetic bistability? We now argue that a clearer picture emerges on using Sir occupancy as a measurement. We drive the system in simulation through the silenced-active-bistable (wild-type) to the silenced-intermediate-bistable region by reducing, thereby recreating the effects of Dot1 inhibition. We also drive the system from the silenced-active-bistable (wild-type) bistable (wild type) to bivalent-monostable region thereby recreating the effects of Dot1 overexpression. Key observations from the model are summarized in the subsections below.
  We have introduced and analyzed a model which captures the distinct mechanisms of cooperativity at the molecular level in the budding yeast silencing system; (a) cooperativity in the histone modification states, (b) cooperative interactions of the Sir proteins in a chromatin-bound-complex. Both of these mechanisms have been proposed by different groups, sometimes as mutually-exclusive ones [1]–[3], [17], [18], [51]. We show that these two mechanisms complement each other in designing epigenetic stability. Our model is minimal—all the interactions included are essential for establishment and inheritance of stable epigenetic fates. The model is build entirely upon known biochemical interactions. Though we explore design principles of the system using inhibition/over-expression/knock-out of key proteins like Sir, Sas2 and Dot1, the model can be extended to include other silencing and active marks—and can inform, like it did for methylation, which of those marks are essential in resolving epigenetic fate, and which simply reinforces one fate over another. Our model makes the following predictions— Experiments focusing on direct measurement of Sir occupancy, preferably at single-cell resolution, as a function of tuning Dot1, Sas2 and Sir2' s enzymatic activity will put to test several observations of the model—the emergence of the bivalent state, and the rich behavior of silencing domains under Sir titration. We have argued that inhibition/over-expression experiments can reveal key information about the engineering design of the system not observable in knock-out experiments. We hope that our model will sharpen experimental questions on the systems biology of budding yeast silencing. For the sake of brevity we have presented results for single perturbations; the model is especially informative for multiple perturbations, for example, double mutants. We are currently investigating such perturbations experimentally. After our paper was submitted, an experimental study by Grunstein lab [94] has further revealed the role of H3K79 methylation in the silencing state of telomeres. Specifically, the transcription-mediated positive feedback of H3K79 methylation was discussed, and indications of our bivalent state where Sir occupancy and H3K79 methylation coexist were provided, though the response of this state to perturbations is still to be studied.
  The general model introduced in this paper is as follows. Define, , , and to be the probabilities of the nucleosome to be in one of the five mutually exclusive states, respectively. The following master equations define the time evolution of the probabilities (lattice model), where summation over nucleosome index is for a suitable neighborhood of interactions with nucleosome — (2) We first analyze the model in the mean-field limit. This implies that all spatial dependence—index and for the position along the chromosome—of all densities of marks are ignored and reactions occur in a well-mixed solution with infinite range of interaction. It is only in this limit that analytical solutions are obtainable. The mean-field analysis reveals all the phases of the model (phase space) and guides the results of lattice simulations. where the lattice is the linear chromosome allowing reactions to occur only in a local neighborhood. Coexistence of silenced/active domains and stochastic spreading of silencing front can be explored in lattice simulations [24]. In the Text S1 under Model I, we present the analytical mean-field solutions for the average densities of marks. All the phase diagrams presented are obtained from thse analytical solutions, which were numerically evaluated on a grid in the four dimensional model-parameter-space. The regions of bistability obtained in the mean-field picture typically shrink for the corresponding lattice model unless the range of cooperative interactions on the lattice is also infinite. Stochastic transition rates between silenced/active domains on the lattice reflect the timescales of stable establishment and heritability of such domains through cell-cycle perturbations. Therefore, we have computed the correlation length (typical domain sizes) of silencing domains in lattice simulations. The construction of the model is minimal in the sense that two of the three cooperative terms—cooperative Sir binding and cooperative deacetylation—are essential for achieving bistability with respect to active and silent states. Cooperative Dot1 binding is non-essential, however, in Text S1 we show that establishment of well-defined silenced and active regions is compromised for the model with cooperative Dot1 rate and tunable basal Dot1 rate. Intuitively, this is because bistability is achieved only for moderate to high basal Dot1 rate, which makes the system fragile to loss of silencing marks (from cell-cycle perturbations) and domains of silencing are disrupted easily by basal Dot1 activity. In the above sense, the model construction is minimal. Though we have presented a specific model for concreteness, many the qualitative features, which includes the number of states and bistabilities supported, remain unchanged under certain plausible additional interactions: a phenomena known as the equivalence of dynamical systems [95]. More precisely, the region of bistability and the cusp bifurcation [95] associated with it, remain even if the model is changed by introducing small perturbations. The transcritical bifurcations [95] in our model; namely, the bivalent state exchanging stability with the active or the silenced state, and the transition between the two types of bistable regions (see Fig. 3); become sharp crossovers in the presence of small positive basal rates and. We have verified that these basal rates need to be be small for robust bistability, see Text S1. Therefore, we expect the conclusions drawn from our simplified treatment to be applicable to a larger class of models. The simulations were performed using Gillespie algorithm on an one-dimensional lattice with lattice sites, with a fixed supply of Sir proteins and a fixed volume of cell. The specific choices for these parameters are quoted below. All rates are measured in units of the constant rate of loss of all marks—this loss rate models all cell-cycle perturbations as a continuous loss. One can, alternatively, consider discrete times where cell division halves all marks stochastically [25]. Our loss rate (set to one) is rather exigent given that all the other rates are in the range O (1) –O (10). The neighborhood of interaction is nucleosomes and is used to compute the local densities of marks. This interaction neighborhood reflects the polymer nature of chromatin where neighboring nucleosomes can come in physical contact often; sets the scale for all our correlation length computation. We only report the qualitative trends of correlation lengths which do not depend on the precise value of. Neighborhood weighting is as follows: For any nucleosome a window of neighboring nucleosomes indexed by is assigned exponentially decaying weights centered on; , where is the absolute value of the distance between and. Therefore, local densities of are where is the indicator variable for the mark. For the mark, , where, . On the lattice the first ten nucleosomes (“telomeric end”) are assumed to be in (silenced) and the last ten to be in (active). This choice simply sets the convention that in the bistable region, the silencing domain forms at the left end and active domain at the right end of the lattice. This boundary condition mimics the telomeric ends and the silencer regions of HML/HMR, which are ‘nucleation centers’ of Sir proteins with high basal rate of Sir binding (i. e. , high). The system is allowed to reach steady state (with stochastically fluctuating but non-propagating front) and the effective determines the zero-velocity line. Parameters: . In each iteration (time step) a reaction is executed, where the reaction is chosen with probabilities proportional to their reaction rates. We perform ten million time steps for each choice of the tuned parameter. Initial time allowed to approach to steady state is two hundred thousand steps. Samples are gathered every time steps, thereby providing 5000 steady state samples to ascertain the effective against the tuned parameter (gray line in Figs. 7 and Fig. 4). The correlation length of marks captures the typical sizes of silenced domains observed in the stochastic evolution of the system and can be computed using the Fourier Transform of the Green' s function for the marks, where is the Fourier space. The Green' s function is, where is the Fourier transform of the local density of marks where the local occupancy is a binary variable. We assume the standard form, where is the correlation length by definition and the' s are fitting constants. is computed by using a multi-taper (Slepian) estimate of the power spectrum using the Chronux toolbox [96]. The larger fluctuation for larger correlation length in the computation is owing to the long-lived nature of larger domains. Parameters: . Steady-state samples (5000) were generated from twenty million time steps sampled every, for each value of and at increments of. We utilize a nonlinear fitting tool in MATLAB ® to fit the estimated power spectrum and compute correlation length.
summary: Epigenetics is the study of heritable phenotypic variations that are not caused by changes in the genotype. Silent Information Regulator (SIR) silencing in budding yeast is an important model system for epigenetics. The standard model of silencing relies on feedback, mediated by chromatin modifications (for example, deacetylation of histone residues) which lead to enhanced recruitment of chromatin modifiers. However, the SIR mechanism is not completely understood and it is important to investigate whether as-yet-undiscovered components alter the systems design in a fundamental way. We address this question using minimal models constructed from experimentally known interactions. Rather than building a detailed network model with parameters to fit for quantitative predictions, we build an effective model and study its bifurcation diagram which leads to robust qualitative predictions on the nature of mutants. This minimal modeling delineates a phase space with qualitatively different epigenetic mechanisms and states; some of which arise from drug/genetic perturbations and exhibit large cell-to-cell variation in chromatin marks. Our methodology can be applied to the study of epigenetic chromatin silencing in other model systems, especially Polycomb silencing, and reveals engineering principles that may be of broad relevance.
keywords: systems biology
  model organisms
  theoretical biology
  biology
  computational biology
  biophysics
title: The Role of Multiple Marks in Epigenetic Silencing and the Emergence of a Stable Bivalent Chromatin State