The Role of Multiple Marks in Epigenetic Silencing and the Emergence of a Stable Bivalent Chromatin State

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.


Introduction
Understanding the design principles of epigenetic silencing phenomena poses challenges to experimentalists and theorists alike, primarily owing to the complexity of interactions between epigenetic marks discovered in recent years [1][2][3]. In the model organism Saccharomyces cerevesiae, since the discovery of telomeric position effect [4,5], intense experimental activity has continued to unravel the many facets of epigenetic silencing in telomeric, Hidden Mating (HML/HMR) loci and ribosomal-DNA (rDNA) regions. However, a complete picture of how epigenetic fate is established, maintained and inherited faithfully is still lacking.
Two universal characteristics of epigenetic phenomena are the switch-like behavior in the expression of genes proximal to silenced chromatin domains and the robust reestablishment of silenced domains through multiple cell cycles. Both of these phenomena impose constraints on any mathematical model of silencing. The nucleation and spreading of silencing in budding yeast, that strongly argues for the bistability of the mechanism, is welldocumented by now [6][7][8], though the mechanism of establishment, maintenance and inheritance is still a matter of active research [9][10][11][12][13][14][15]. At least two different mechanisms of silencing spreading have been proposed in the literature. We call the first mechanism the polymerization model of spreading [16][17][18][19], whereby interactions of proteins amongst themselves lead to a spread of silencing along the chromatin. We call the second mechanism the histone modification feedback mechanism, whereby the recruitment and spreading of silencers is controlled by histone modifications and silencing/transcriptional feedbacks; various anti-silencing marks etc. block entire regions of chromatin from being silenced. Are these two mechanisms redundant? We explore these mechanisms and argue that they play synergistic role in the robustness of the system; an important conceptual outcome of our approach.
Multiple chromatin modifications (marks) are thought to be responsible in carrying epigenetic information [2,3,20]. Antisilencing histone marks, primarily acetylation and methylation, are known to be key players in budding yeast silencing [21][22][23]. What role does multiple histone marks play in the dynamics and stability of the silencing phenomena? Why are multiple anti-silencing marks needed? Despite recent efforts at mathematical modeling of epigenetic chromatin modification in budding yeast [24][25][26] and in other systems [27][28][29], theoretical understanding on the role of multiple marks in epigenetic bistability remain incomplete. Moreover, currently there is a gap in the literature on mathematical models that explain epigenetic bistability in terms of local biochemical interactions. Building informative models is not necessarily impaired by our lack of knowledge of the biochemistry or reaction rates, and non-trivial qualitative predictions can be obtained from rather minimal models providing an integrative understanding of experiments in the field. We illustrate the role of multiple marks in epigenetic bistability and robustness of silenced/active regions using a minimal model.
Experiments largely rely on the behavior of knock-out mutants to unravel the components of the yeast silencing system. Less wellstudied are the inhibition/enhancement (of enzymatic activity) or over/under-expression (of proteins), though overexpression of Sir proteins and Dot1 have been examined [30][31][32][33][34][35][36]. We argue that compared to knock-out mutants, inhibition etc. is better suited at revealing the engineering design principles of the system. Using theory and simulations we infer the pivotal role of Dot1 in establishing stable silencing domains and resolve seemingly contradictory experimental observations on this matter. We also argue that as long as the set of known interactions included in our model are retained, our overall conclusions are robust to the addition of more complex interactions, either known or yet to be discovered.

Primer on molecular biology of budding 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 matingtype Information Regulator (SIR) proteins have been identified to play a pivotal role in heterochromatin formation [30,[37][38][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][42][43][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,44,[50][51][52]. 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 hyperacetylated at particular histone tail sites but also tend to be hypermethylated 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][60][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][65][66]. In fact, such positive feedback from transcriptional elongation into establishing transcriptionally-active marks is not uncommon. For example, H3K4 methylation [67][68][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][72][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][82][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

Author 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.
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.

The model and its phase space
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 i:S i , A i , M i , E i and U i . The S state is Sir-complex-bound (without making a distinction of the various Sir proteins). The A state is acetylated histone (H4K16Ac). Multiple histone tails on the same nucleosome can get acetylated-we model the average level of acetylation. The M is the methylated state (H3K79Me), and we treat the multiple levels of methylations on the average. The U is unmodified state-it is neither acetylated, methylated or Sirprotein bound. The E 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 j are for nucleosomes within a local neighborhood of nucleosome i: N Basal Sir complex binding at rate r 0 : U i ? r 0 S i : N Basal Dot1-mediated methylation at rate b 0 : N Basal Sas2-mediated acetylation at rate a: N Cooperative deacetylation by Sir2 at rate C: N Cooperative Sir binding rate r: ''Silencing polymeriza- N Global and local rates of loss of marks g that converts all other states to U.
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 a, the cooperative Dot1 activity b, the cooperative Sir2 activity C and the cooperative Sir binding r. 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 (E) and silent (S) 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 S, the active state with weight predominantly in the state E, the intermediate state with weight distributed in U, A and M, and the bivalent state with weight distributed predominantly over S and E. 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 C and r 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 Sir titration effect results in a local rate r of Sir recruitment The supply of Sir proteins is limited in order to prevent ectopic silencing in wild type. The reaction rates r 0 and r, 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 selfadjusts 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 S the the total number of Sir proteins, and s the density of Sir proteins, the equation that determines r is where V is the volume of the cell nucleus. As silencing spreads, S bound increases, and r decreases. Such a self-adjusting r 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.

Sir2 inhibition leads to bivalent state
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 C in our model. Because of the titration effect, whenever a single parameter is changed the parameter r 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 (b) and Sas2 activity (a). 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 (C), r 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-N Inhibition of Sir2 activity eventually drives the wild-type bistable system, where silencing and active states are stable, into a monostable region where bivalent states are stable.
N The inhibition of Sir2 activity leads to an excess of ambient Sir proteins.
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 (r) increases locally for decreasing deacetylase activity (C) 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 occupancya 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 S marks for the system self-tuning r 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.

Role of Dot1 in the systems design
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 dot1D strain [90]. Experiments [23,35,54,59,91] have focused on both dot1D and overexpression of Dot1, but not inhibition. In this section, we summarize the bearings of our model on Dot1 perturbations.
The dot1D 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 A domains to be moderately transcribed. We observe in Fig. 7 that there is a critical value of b, which is given by b crit~( 1za) 2 a(2za) (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 nonzero basal rates of Dot1 activity (b 0 ) and Sir binding (r 0 ) 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 (b) 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 b, thereby recreating the effects of Dot1 inhibition. We also drive the system from the silenced-active-bistable (wild-type) bistable (wild type) to bivalentmonostable region thereby recreating the effects of Dot1 overexpression. Key observations from the model are summarized in the subsections below.
Inhibition of Dot1 leads to defective region boundaries. Dot1 inhibition (bvb crit ) leads to lower local density of Sir occupancy in telomeres but higher fraction of Sir proteins to be chromatin-bound, leading to Sir depletion and illdefined silenced domains, see Fig. 8 and Fig. 9. In our model, Dot1 inhibition leads to the system crossing the silenced-activebistable (wild type) to the silenced-intermediate-bistable region. In model simulations, we show that the wild type system has higher density of Sir occupancy-Dot1 inhibition leads to lowering of local Sir occupancy in the silenced-intermediate-bistable region, see Fig. 8. However, the fraction chromatin-bound Sir exhibits the opposite trend with increased fraction for Dot1 inhibition. More chromatin-bound Sir leads to depletion of ambient Sir proteins.
The consequence of Sir depletion in the dot1D strain by the telomeres can lead to reduced silencing in hidden mating loci [69,91]. Strikingly, the lower density of Sir occupancy in dot1D has recently been observed in experiments [54] in agreement with our result.
In the dot1D strain, the silencing regions are variegated and the domains are ill-defined, see Fig. 9. Methylation-and multiple histone marks in general-play a key role in establishing and maintaining heritable silenced and active domains distinguished by sharp boundaries. The effect of tuning Dot1 activity from null is presented in Fig. 7, corresponding to Fig. 4. Recall that the self adjustment of r(s ambient ) owing to titration effect maintains the steady-state system on the zero-velocity line. A sample of the relevant cuts through the phase diagram is also shown in Fig. 7.
In order to characterize both the bistable regions further, we measure the correlation length of the S marks as the system travels along the zero-velocity line; see Fig. 10 and refer to Fig. 7. For clarity of the discussion below, we quote numerical values of b specific to the section of phase space presented. Ovbiously, these values are in general not meaningful. With that caveat-the range b[½0,8 in Fig. 10 and Fig. 7 is relevant for the present discussion, where the crossover between the two bistable regions is Figure 6. Effect of Sir2 perturbations: Correlation length of silenced mark (S) state of the system along the zero-velocity-line in Fig. 4. Sir inhibition pushes the system to the bivalent-stable phase. The transition from active-silenced-bistable state (wild-type) to the bivalent-stable phase is approximately at C&4. The correlation length is high, as expected, in the bistable region where silencing domains are established, however, the system continues to enjoy relatively strong correlation lengths in the bivalent-stable phase. The scale on the y-axis depends on the (unknown) biophysical parameters of the wild-type system, we are only reporting the trend. doi:10.1371/journal.pcbi.1003121.g006 approximately at b&1:5. In the wild type (1:5 * v b * v 9) silenced state is stable and we expect strong correlations in Sir occupancy, which should survive Dot1 inhibition (lowering b) even after active state loses stability for b * v 1:5). Accordingly, we observe in Fig. 10 that the correlation length drops but does retain moderate values for strong Dot1 inhibition. The drop is explained by the defective establishment of the silencing domain in telomeres, see Fig. 9.
The survival of correlations indicates that Dot1 perturbation does not necessarily eliminate silencing in shorter loci like HML/ HMR with boudary elements that may nucleate silencing, provided their length is comparable to the correlation length of Sir occupancy [36,54,59,90]. We cannot determine this correlation length without further knowledge of the wild-type parameters and other details of the chromatin polymer. Nevertheless, we make a strong case for reporting the local level of Sir occupancy in the loci, as opposed to reporting transcription to study the role of Dot1.
Over-expression of Dot1 leads to patchy silencing and bivalent state. Over-expression of Dot1 eventually pushes the system to the bivalent state, which is stable, and therefore, potentially heritable. We observe that the state enjoys large spatial correlations, implying that Dot1 overexpression also leads to patchy silencing, but unlike Dot1 inhibition, large patches of silencing domains established stochastically may be very long-lived in the cell. Such domains may be easier to establish and maintain in smaller chromatin regions, like HML/HMR, providing a possible explanation as to why Dot1 overexpression may have a weak effect at such loci.
In the case of over-expression of Dot1 the system is driven out of the bistable region and is eventually pushed to a bivalent- monostable regions. In the vein of the previous discussion, refer to Fig. 10 and Fig. 7 once again. The range b * > 9 in Fig. 10 and Fig. 7 is currently relevant. In the wild type (1:5 * v b * v 9), we have strong correlations as expected. However, b * > 9 the system loses bistability and is bivalent-monostable. The correlation length in this region is found to be comparable as in the bistable region (wild-type) well beyond the crossover. The strong persistence of correlations in spite of loss of stable silencing domain is a nontrivial prediction of our model. It suggests that for Dot1 overexpression silencing marks are established in patches that are long-lived and possibly heritable. A mechanistic explanation is: In the bistable region (wild-type), silencing domain is established and ambient Sir protein concentration is low. In the bivalentmonostable region no such stable domain exists. However, increased effective rate of cooperativity in Sir recruitment (r) increases for increasing Dot1 activity (b) along the zero-velocity owing to depletion of chromatin-bound Sir proteins, see Eq. 1. Therefore, Sir proteins that are chromatin bound stochastically form long-lived patches.
Dot1 and Sas2 perturbations lead to distinct phases. It is now germane to contrast the effects of Dot1 perturbation and Sas2 perturbations in the light of robustness caused by multiple histone modifications. Both of the acetylation and methylation marks we consider are active marks. However, Sas2 inhibition pushes the Figure 8. Effect of Dot1 inhibition on local Sir density and net chromatin-bound Sir: The density of Sir occupancy, determined from simulations, in telomeres (blue solid line) and fraction of Sir chromatin-bound (green dashed line) as a function of increasing Dot1 activity (b) starting from dot1D upto wild-type values. The system is bistable for all values b explored here, however, the nature of the bistability transitions from intermediate-silenced-bistable to active-silenced-bistable at b&1:5. Note that for decreasing b, the local density of Sir protein decreases though the net fraction of chromatin bound Sir protein, counter-intuitively, increases. This is because compact silenced region are lost owing to Dot1 inhibition. doi:10.1371/journal.pcbi.1003121.g008 wild-type system to regions of phase space where bistability loses robustness because the region of bistability narrows (the zerovelocity line runs towards a cusp), see Ref. [24,25] and Text S1. Dot1 inhibition pushes the system to regions where it continues to enjoy bistability, and the stable silenced state persists. However, a heritable pattern of silenced and active domains is compromised in Dot1 perturbations, as discussed above. Intuitively, the polymerization model of Sir recruitment which ensures continuous Sir occupancy also necessitates a histone modification feedback mechanism to counteract ectopic spread. The latter is, owing to the engineering limitations of local interactions and cell-cycle perturbations, easier to achieve using multiple histone marks [25].

Discussion
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 mutuallyexclusive ones [1][2][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-N Dot1 over-expression and Sir2 inhibition can push the system to a novel bivalent state exhibiting patchy silencing. We argue that the contradictory role of Dot1 in HML/HMR loci and telomeres can be reconciled by studying the correlations of Sir occupancy in the bivalent state. Sir occupancy on the chromatin, instead of transcriptional activity, is a better reporter of the effects of Dot1 and Sir2 perturbations.
N We compute Sir occupancy in simulations and compare its qualitative behavior to existing experimental work [36,54,59,90]. Specifically, we observe that Dot1 inhibition can lead to higher fraction of chromatin-bound Sir and lower density of Sir in a loci owing to poor establishment of silencing domain, thereby causing Sir depletion. On the other hand, Dot1 overexpression can result in compact patches of Sir binding established stochastically. These patches are long-lived and may be heritable.
N In comparing the behavior of Dot1 and Sas2 inhibition, we predict that while Dot1 inhibition does not eliminate stability of silenced domains, Sas2 does. We argue that methylation and acetylation play distinct roles in meeting the design requirements of bistability of states and establishment of inheritable domains of those states.
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.

Materials and Methods
The general model introduced in this paper is as follows. Define P i (S), P i (A), P i (M), P i (E) and P j (U) to be the probabilities of the nucleosome i to be in one of the five mutually exclusive states, S,A,M,E,U respectively. The following master equations define the time evolution of the probabilities (lattice model), where summation over nucleosome index j is for a suitable neighborhood of interactions with nucleosome i- We first analyze the model in the mean-field limit. This implies that all spatial dependence-index i and j 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 meanfield 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 b~0 and tunable basal Dot1 rate b 0 . Intuitively, this is because bistability is achieved only for moderate to high basal Dot1 rate b 0 , 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 r 0 and b 0 . 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 L lattice sites, with a fixed supply of Sir proteins S tot and a fixed volume of cell V . 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 N 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; N 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 N. Neighborhood weighting is as follows: For any nucleosome i a window of N neighboring nucleosomes indexed by j is assigned exponentially decaying weights w ij centered on i; , where Dx ij D is the absolute value of the distance between i and j. Therefore, local densities of On the lattice the first ten nucleosomes (''telomeric end'') are assumed to be in S (silenced) and the last ten to be in E (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 r 0 ).
The system is allowed to reach steady state (with stochastically fluctuating but non-propagating front) and the effective r(s ambient ) determines the zero-velocity line. Parameters: V~50, L~200,S tot~1 00,N~10. 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 10|L time steps, thereby providing 5000 steady state samples to ascertain the effective r(s ambient ) against the tuned parameter (gray line in Figs. 7 and Fig. 4).
The correlation length of S 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 functionG G(k) for the S marks, where k is the Fourier space. The Green's function is G(k)~r r(k)r r({k), wherer r(k) is the Fourier transform of the local density r(s) of S marks where the local occupancy s is a binary variable. We assume the standard form G G(k)~C 1 k 2 zj {2 zC 2 , where j is the correlation length by definition and the C's are fitting constants.G G(k) 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: V~150,L~400, S tot~3 00,N~10,r 0~0 :01,b 0~0 :01.
Steady-state samples (5000) were generated from twenty million time steps sampled every L|10, for each value of b[½0,20 and C[½0,10 at increments of 0:5. We utilize a nonlinear fitting tool in MATLAB H to fit the estimated power spectrum and compute correlation length j.

Supporting Information
Text S1 We have relegated the mathematical details of the model in Text S1. Analytical solution of the model in the mean-field limit is presented therein, and the nature of the different phases is summarized. We also compare related models with fewer or more biochemical interactions (parameters) in order to establish that the interactions we have considered are minimal for the engineering design criteria of epigenetic silencing. We also elaborate on Sir perturbations and Sas2 perturbations and discuss the nature of bifurcations/crossovers observed in the phase space. (PDF)