Tunable, division-independent control of gene activation timing by a polycomb switch

SUMMARY During development, progenitors often differentiate many cell generations after receiving signals. These delays must be robust yet tunable for precise population size control. Polycomb repressive mechanisms, involving histone H3 lysine-27 trimethylation (H3K27me3), restrain the expression of lineage-specifying genes in progenitors and may delay their activation and ensuing differentiation. Here, we elucidate an epigenetic switch controlling the T cell commitment gene Bcl11b that holds its locus in a heritable inactive state for multiple cell generations before activation. Integrating experiments and modeling, we identify a mechanism where H3K27me3 levels at Bcl11b, regulated by methyltransferase and demethylase activities, set the time delay at which the locus switches from a compacted, silent state to an extended, active state. This activation delay robustly spans many cell generations, is tunable by chromatin modifiers and transcription factors, and is independent of cell division. With their regulatory flexibility, such timed epigenetic switches may broadly control timing in development.

(left). All histograms show that while the perturbations affect the all-or-none activation probability for the initially inactive alleles (top), the pertubations have no effect on the expression maintenance nor magnitude of the initially active alleles (bottom). (B) Mean percentage of cells remaining Bcl11b YFP-negative (n = independent experiments, *p < 0.05, **p < 0.01, two-sample t-test, two-tailed). (C) Relative mRNA levels of Eed were measured by qPCR. (D) Mean percentage of cells remaining YFP-negative after DN2 Bcl11b RFP+/YFPmonoallelic cells were transduced with retroviral constructs and recultured for 3 days (two-sample t-test, one-tailed, **p < 0.01, n = 3 independent experiments, error bars = 95% confidence interval). (E) Progenitors from Figure 2D (day 3) were stained with Annexin V to detect apoptotic cells. Segmented populations were analyzed using image processing tools in MATLAB, and specific features such as object area, perimeter, and CFP intensity were collected. Approximately 10% of the total individual cells were then manually labeled as live or dead and fed into classification tree machine learning algorithm to generate a classification model. The rest of segmented cells are classified subsequently via the trained model. (B) Mathematical model describing the population dynamics of Bcl11b RFP+/YFPcells transduced with c-Myc and empty vector (EV). The model includes a population of live cells (X) with a birth rate k b and a death rate k d that generates the dead cell population (Y). This population then has a permanent clearance rate of δ, indicating the process in which CFP degrades in these cells and the dead cells become undetected. (C) Experimentally determined decay of dead cells' detectability in time-lapsed movies. Individual dead cells were followed until their CFP level completely diminished and was undetectable by the segmentation algorithm, and the elapsed time was recorded. Thirty different individuals cells were recorded for each c-Myc and EV populations. Data was fit to an exponential decay function P(τ) = e -δτ with δ = 0.032 per frame for c-Myc and δ = 0.030 per frame for EV population. Parameter F dictates how sensitve nucleosome compaction affinity is to demethylation (i.e. when F is high, the compaction affinity is only moderately affected by demethylation; see Mathematical Appendix for more details). The activation energy barrier (Ea) is defined as the potential energy (V) height between the local maximum and local minimum of the potential energy landscape. Each potential curve was plotted with demethylation parameter set to 20 hrs -1 and methylation rate parameter as indicated by the curve's color see Mathematical Appendix). A nucleosome's methylation rate increases with the number of methylated nucleosomes in the system. (B) Average switching times as a function of methylation (β) and demethylation rate constant (α) ratio for the methylation read-write (black) and methylation-compaction (red) models. Sensitivity coefficient (ΔlogY/ΔlogX) for each plot was calculated by taking the slope of the linear fit y = ax+b for the methylation model data set and the last 5 data points for the compaction model.

Introduction
To understand the timed epigenetic switch controlling the Bcl11b activation, we used mathematical modeling to analyze a series of candidate biophysical mechanisms. This mathematical modeling analysis seeks to uncover the essential emergent properties of the switch, namely (1) its irreversible, all-or-none nature; (2) its long, stochastic time delay; (3) the heritability of its inactive and active states over DNA replication; and (4) its tunability with respect to changes in H3K27me3 levels and modifying-enzyme activity.
We consider two main candidate models. In the methylation read-write model (M), individual nucleosomes within in a one-dimensional lattice can be methylated or unmethylated. Gene expression is assumed to occur when the total fraction of methylated nucleosomes in this lattice falls below a threshold value. In the methylation compaction (MC) model, individual nucleosomes are also methylated and demethylated; in addition, these nucleosomes also interact to form a compacted assembly with rates dependent on their H3K27me3 state. Unlike the methylation read-write model, gene expression does not depend directly on H3K27me3 levels, but on the compaction state of the nucleosome assembly, which in turn depends on methylation states of individual nucleosomes. Both models explicitly model DNA replication as a process involving random segregation of modified nucleosomes into daughter strands. From our analysis, we find that the methylation-compaction mechanism, but not the methylation read-write mechanism, explains the emergent behaviors of the timed epigenetic switch controlling Bcl11b activation, and thus represents our favored model.

Model I: The Methylation Read-Write Mechanism (M)
Here, we adopt a standard framework for histone modification dynamics previously shown to generate multi-stability (Angel et al., 2011;Dodd et al., 2007). In this model, individual nucleosomes reside in a one-dimensional lattice, and exist in two states, a methylated state, corresponding to an H3K27 tri-methylated state, and demethylated state. We do not describe multiple demethylated states in our model (i.e. mono-methylation, di-methylation, and an un-methylated state), though our analysis, together with previous work (Dodd et al., 2007), indicates that our main conclusions should also hold in more complex models with additional states. As with previous models, the methylation rate of a given nucleosome depends on the number and distance of methylated nucleosomes in its vicinity, reflecting observations that PRC2 can bind and be activated by H3K27me3-marked nucleosomes to write H3K27me3 on neighboring nucleosomes. The positive feedback generated by this methylation read-write mechanism provides a basis for bi-stabilty in this model. Here, demethylation is taken to occur at a first order rate. We assume there is no spontaneous methylation in the absence of existing methylated nucleosomes; thus, once all nucleosomes are demethylated, the system irreversibly enters an activated state.
Methylation. We explicitly model mark binding and methyltransferase activities of the PRC2 complex, as well as the methylation state of each individual nucleosome. Take the gene locus to a linear array of N nucleosomes. Let ! = 1. . % denote the index for the ith nucleosome, and let & ! be its H3K27 methylation state. & ! = 0 denotes the de-methylated state while & ! = 1 denotes the methylated state.
Let )′ and ) denote the transitions between the methylation state and demethylation state, respectively. The model is set up as follows: For ! ∈ {1, . . , %}: The parameter L can be interpreted as the 'reach' of the PRC2 complex to neighboring nucleosomes. A large value of L indicates a long length scale for nucleosome interactions. This effect is set to have a gaussian shape so that nucleosome closest to the anchored PRC2 complex has the highest methylation rate. Similar distributions of activity have been reported for artificially tethered enzymes (Hass et al., 2015), as well as for histone modifications around transcription factor binding sites (Heinz et al., 2010). Moreover, we assume periodic boundary conditions for the onedimensional lattice, though similar results were observed with other non-repeating boundary conditions (not shown).
Cell division. To model the transmission of histone marks across cell divisions, we assume that methylated nucleosomes segregate randomly to the two daughter DNA strands upon replication; thus, each nucleosome position has a one-half probability of inheriting a nucleosome that is methylated. Experimental evidence suggests that approximately half of total global H3K27me3 partitioning of parental marks to the subsequent generations (Alabert et al., 2015).
From stochastic simulations of this model, we find that this methylation read-write mechanism can generate a time-delayed, stochastic switches from an inactive H3K27me3-high state to an active state without H3K27me3 (Fig. 4A); however, switching times are hypersensitive to mild changes to methylation and de-methylation rates (Fig. 4B), and therefore inconsistent with the graded changes in switching times observed upon inhibition of PRC2 methyltransferase or Kdm6a/b demethylases (Fig. 2). To understand the origins of this hypersensitivity, we re-formulate this model using a chemical kinetics framework amenable to analysis using transition state theory. To do so, we first consider the limit where < → ∞, such that each H3K27me3-bound PRC2 methylates all other unmethylated nucleosomes with the same reaction rate. In this limit, we can completely describe the state of the system by a single variable, the number of methylated nucleosomes N'. As the rates of adding or subtracting one methylated nucleosome from the system would reduce to become a function of N', independent of spatial arrangement. Consequently: where and % ) is the total number of nucleosomes. The master equation describing the time evolution of this system is given by: where & + is the probability of having N' methylated nucleosomes. When the total number of nucleosomes is large, we can approximate the number of methylated nucleosomes to be a continuous variable H. In this limit, we can rewrite the master equation as Fokker-Planck equation: where, we have ignored third and higher order terms, and where: Given the velocity and diffusion constants for this system as a function of methylated nucleosome number, the switching of the system is essentially given by the first-passage time of the system to reach the absorbing state x = 0. A closed-form solution of this first-passage time distribution for the given rate functions is hard to obtain; Nevertheless, we note that our system operates in the regime where the timescales of individual methylation and demethylation reactions are much shorter than switching times for this system. In this regime, switching times are well described by the Kramer's theory for escape of a Brownian particle over a potential well (Kramers, 1940), and would thus approximately scale exponentially with the height of a potential energy barrier. We can obtain the functional form of this potential barrier by relating it to the velocity function: From equations (3), (4), and (7), we can then integrate the system to explicitly derive the potential function: Where Q is an arbitrary number. A plot of O(H) is demonstrated in Fig. S3C. The energy landscape possesses a local minimum at a nonzero value of H, indicating the metastable state. The landscape has a local maximum near H = 0. State switching occurs when the system reaches the absorbing state x = 0. Thus, we define R . as the height in O(H) between the local maxima and the metastable state minima. When we plotted this potential energy for different values of β (Fig. S3C), we found that moderate changes in β led to significant changes in potential well height. As switching rates scale roughly exponentially with well height, we would expect this system would show extreme sensitivity in switching times with respect to changes in methylation rate changes.

Model II: The Methylation Compaction Mechanism (MC)
Because the methylation read-write mechanism above does not account for the tunable characteristics of the Bcl11b activation timing switch, we considered a second model, where histone methylation facilitates interactions between nucleosomes to enable the stable maintenance of a repressed, compacted chromatin state at the Bcl11b locus. There are multiple mechanisms by which H3K27me3 could facilitate interactions between nucleosomes: H3K27me3 could recruit polycomb repressive complex 1 (PRC1), which could oligomerize through contacts on its Bmi1 or Phc subunits (Eskeland et al., 2010;Gray et al., 2016;Isono et al., 2013;Kahn et al., 2016), or undergo weak, multivalent interactions on its Cbx2 subunit that result in liquid-liquid phase separation (Howard, 2001;Larson et al., 2017). Alternatively, H3K27me3 could modulate affinities of weak multivalent interactions between nucleosomes (Gibson et al., 2019), and thereby modulate their ability to phase separate.
The methylation compaction model consists of two main modules: (1) a H3K27 methylation and demethylation mechanism, and (2) a dynamic chromatin decompaction mechanism linked to H3K27me3 modification state that ultimately underlies gene switching. In our description of compaction dynamics, we do not explicitly model the spatial extent of the compacted nucleosomal assembly; instead, we adopt a mean-field approach that is established in models of cytoskeletal polymer dynamics (Erickson and Pantaloni, 1981;Jackson and Berkowitz, 1980). With this approach, the numbers of un-methylated and methylated nucleosomes within a compacted assembly are given by C and C' respectively, along with those outside the assembly are given by D and D' respectively. As a result, the dynamical system is described by four states: 1) Compacted-Methylated 2) Compacted-Demethylated 3) Decompacted-Methylated and 4) Compacted-Demethylated: Here, C', C, D', and D denote the number of nucleosomes in these states, respectively. k1 to k8 denote the transition rates between them, which will be defined below. The gene is taken to be activated when all nucleosomes exist in a de-compacted state. The two mechanisms are intertwined so that methylation states affect compaction rates and vice versa. Detailed descriptions of the rates are given below: Methylation. In this model, un-methylated nucleosomes convert into a methylated state with a firstorder rate constant β. We assume this rate constant is the same regardless of whether nucleosomes are inside or outside the compacted assembly. Methylated nucleosomes convert into a demethylated state with a rate constant of : if the nucleosome is outside the assembly (D'), or a lower rate constant of S:, (S < 1) if the nucleosome is inside the assembly. This lower rate constant assumes that the demethylation reaction is less efficient on compacted nucleosomes, possibly due to competition for demethylase binding by compaction proteins, or due to the exclusion of demethylases through steric occlusion or phase separation. The rates describing these reactions on the four nucleosomal species are given by: Where: H3K27 methylation and demethylation rates are based on the catalytic activity of the Ezh2 subunit of the PRC2 complex and Kdm6a/b demethylases, respectively. Specifically, these rate constants were chosen to represent the conversion between H3K27me2 and H3K27me3. For simplicity, we do not model the H3K27me-binding dependent H3K27 methylation activity previously described (Margueron et al., 2009), though we show below that explicit modeling of this read-write effect would not significantly alter the conclusion of the model. Kdm6a/b demethylate H3K27me3 (Agger et al., 2007), and to our knowledge no cooperative activity of these complexes have not been reported.
Compaction. We adopt a mean-field description of the compacted nucleosomal assembly, following kinetic models of multi-stranded cytoskeletal polymer assembly (Howard, 2001). This description assumes that the nucleosome assembly is a roughly spherical structure held together by weak, multivalent interactions between individual nucleosomes, and can add or lose individual nucleosomes at its surface. Both methylated and demethylated nucleosomes can incorporate into the assembly; thus the assembly has a total size of: where U and U′ represent the number of methylated and demethylated nucleosomes in the assembly, respectively. Unlike other polymer models (MacPherson et al., 2018;Nuebler et al., 2018), we do not explicitly model physical connections between nucleosomes due to DNA; such connections would be expected to result in a spatial dependence of reaction rates within this chromatin domain; however, as the entire domain (100 nucleosomes) has a length scale greater than the persistence length of chromatin (~15-20 nucleosomes, from (Arbona et al., 2017)), and would thus enable free interactions between non-neighboring nucleosomes, we would expect the essential properties of our minimal model in a more realistic physical model that incorporates nucleosome connectedness.
The addition and removal of methylated and demethylated nucleosomes from the assembly is described by the following rate equations: Where: Here, methylated nucleosomes incorporate into the compacted assembly with a rate constant V; however, importantly, demethylated nucleosomes can also incorporate into the assembly with a reduced rate constant WV (where W < 1). The effect of methylation state on compaction rate is experimentally observed in instances such as recruitment of PRC1 complex by H3K27me3 marks (Kahn et al., 2016). The complex's subunits such as Ring1B and Phc-1 have been shown to be important in chromatin compaction and gene silencing (Eskeland et al., 2010;Francis et al., 2004;Isono et al., 2013). However, as PRC1 recruitment is not the only compaction mechanism in vivo, and because PRC1 can bind to nucleosomes independently of H3K27me3 (Francis et al., 2004), this model treats methylation as only a part, but not solely responsible for chromatin condensation. In choosing rate constants; we assume that compaction and decompaction is faster than histone methylation and demethylation rates, though timescales for both processes are assumed to be much faster than that for cell division. Fast compaction kinetics relative of modification is supported by in vitro studies of H3K27me3 methylation and demethylation kinetics, as well as in vitro DNA compaction by HP1α and chromatin condensation experiments (Kristensen et al., 2011;Ladoux et al., 2000;Larson et al., 2017;Sneeringer et al., 2010).
The reaction rates for nucleosome incorporation (loss) scales with assembly size as ~U ) % ( (~U ) % ' ( ), as these reactions only take place on the surface of the assembly. Assuming a compacted nucleosome complex is spherical, the compaction rate would thus be proportional to the surface area. Likewise, the decompaction rate is also proportional to the surface area but reversely proportional to the total number compacted nucleosome in the complex.
In this description, there is a critical threshold number of compacted nucleosomes, U 6 , below which the complex is thermodynamically unstable. The existence of a minimal nucleus size is a fundamental property of phase-separated assemblies held together by weak-multivalent, whereby addition of a new subunit to an already formed complex is thermodynamically more favorable than formation of the initial nucleus itself (Erickson and Pantaloni, 1981;Jackson and Berkowitz, 1980). Below this critical threshold U 6 , the compacted assembly disintegrates, and gene turns on.
Cell division. Heritability of histone marks and chromatin states are crucial in maintaining gene expression states across cellular generations. As with the methylation read-write model above, we assume that methylated nucleosomes partition randomly between two daughter strands upon replication; the total number of nucleated nucleosomes is then obtained by sampling from binomial distribution with one half probability and N equal to the total number of nucleosomes at the point of DNA replication. Furthermore, we assume that compacted nucleosomes persist within a compacted assembly reside upon passage of DNA polymerase. This model feature assumes that new nucleosomes rapidly incorporate into a compacted assembly after passage of DNA polymerase; however, in the subsequent version of this model below, we will relax this assumption to allow for disruption of compaction state by DNA polymerase passage (see below).
From Monte-Carlo simulations, we found that this dynamic methylation compaction model can recapitulate all the essential emergent properties of the Bcl11b activation switch. Specifically, this model shows the following dynamic properties: 1) Irreversible all-or-none switching to an H3K27me3-low, de-compacted state. From simulations, we found that the system adopts a stable compacted assembly of nucleosomes with higher H3K27me3 marking density, but switches abruptly to a de-compacted state with lower H3K72me3 levels. As there is no re-nucleation of the compacted assembly after its elimination, this de-compacted state represents an absorbing, permanently active expressing state. The abrupt decrease in the H3K27me3 levels arises because compacted nucleosomes demethylate at a lower rate; thus, upon total decompaction, the percent of methylated nucleosomes lowers to a new steady state level.
2) Noise induced gene activation. Transition to the completely decompacted state, or gene activated state, occurs via stochastic deviation of the system from its compaction meta-stable state. Activation is triggered when the system reaches below the threshold number of compacted nucleosomes.
3) Tunable activation rates. The model is able to generate a gene switch with slow, tunable activation rate. Delay in activation is in order of days and can be finely adjusted by modifying methylation and demethylation rates, and/or changing H3K27me3 levels at the gene locus (Fig.  4E), as experimentally observed (Fig. 2). This ability to tune activation rates by changing H3K27me3 densities distinguishes this methylation compaction model from the methylation read-write model above, and thus represents a more plausible model for describing the activation mechanism of this switch. Why is this model uniquely tunable? In this model, locus decompaction and gene activation are determined by a dynamic balance between rates of nucleosome entry or exit from a compacted assembly. The system still be sensitive to changes in these rates; however, as demethylated nucleosomes can still enter and exit a compacted assembly at a reduced rate, changes in the fraction of demethylated nucleosomes would cause a fine change in these entry or exit rates, and thus give rise to a plausible tuning parameter for controlling activation timing.

4)
Division-independent timing control. When the cell cycle length is changed in this model, activation kinetics remain largely unaffected, implying that the methylation-compaction mechanism functions as a cell division-independent delay timer. These conclusions hold, as long as the dynamic methylation and compaction mechanisms operate on timescales much faster than the cell cycle length.
To gain insights into the origins of tunability for the methylation compaction model, we adopt an approach, where we reduce this problem to using the Fokker-Planck approach, as utilized to analyze the methylation read-write mechanism (Fig. S3D,E). The full system with both methylation and compaction reactions would correspond to diffusive motion of a particle in a three-dimensional state space describing both chemical and physical states of nucleosomes. However, to simplify this problem to gain intuition, we will first take the methylation and demethylation reactions to be fast compared to the compaction and de-compaction reactions, such that the system can be described a single parameter U ) , corresponding to the total number of compacted nucleosomes. At any given time, the number of methylated and demethylated nucleosomes in the compacted state is at quasisteady state, with values: and U = S: 5 + S: Similarly, assuming that the system is at quasi steady state, the number of methylated and demethylated nucleosomes in the uncompacted state is given by: and J = : 5 + : • J ) Let % ) = U ) + J ) . With this approximation, the averaged rate of adding or removing a nucleosome from the compacted assembly is then given by: Let the total number of compacted nucleosomes U ) be H. By writing down the master equation for this system, and by further applying the Fokker-Planck approximation, as performed in (5) and (6) As before, we define a potential energy for this system: The analytical solution for the potential energy O(H) for the methylation compaction model is:  Fig. S3D-E. We found that increasing methylation rate results in a much more attenuated increase in activation energy R . with the methylation compaction model. This confirms that the improved switching rate tunability in the MC model stems from the decreased sensitivity to changes in activation barrier height by methylation rate. This result intuitive explains why this system shows significantly more graded changes in switching times when methylation rates are changed.