Mechanistic stochastic model of histone modification pattern formation
© Anink-Groenen et al.; licensee BioMed Central Ltd. 2014
Received: 10 July 2014
Accepted: 2 October 2014
Published: 27 October 2014
The activity of a single gene is influenced by the composition of the chromatin in which it is embedded. Nucleosome turnover, conformational dynamics, and covalent histone modifications each induce changes in the structure of chromatin and its affinity for regulatory proteins. The dynamics of histone modifications and the persistence of modification patterns for long periods are still largely unknown.
In this study, we present a stochastic mathematical model that describes the molecular mechanisms of histone modification pattern formation along a single gene, with non-phenomenological, physical parameters. We find that diffusion and recruitment properties of histone modifying enzymes together with chromatin connectivity allow for a rich repertoire of stochastic histone modification dynamics and pattern formation. We demonstrate that histone modification patterns at a single gene can be established or removed within a few minutes through diffusion and weak recruitment mechanisms of histone modification spreading. Moreover, we show that strong synergism between diffusion and weak recruitment mechanisms leads to nearly irreversible transitions in histone modification patterns providing stable patterns. In the absence of chromatin connectivity spontaneous and dynamic histone modification boundaries can be formed that are highly unstable, and spontaneous fluctuations cause them to diffuse randomly. Chromatin connectivity destabilizes this synergistic system and introduces bistability, illustrating state switching between opposing modification states of the model gene. The observed bistable long-range and localized pattern formation are critical effectors of gene expression regulation.
This study illustrates how the cooperative interactions between regulatory proteins and the chromatin state generate complex stochastic dynamics of gene expression regulation.
KeywordsChromatin structure Histone modification patterns Epigenetics Stochastic mathematical model Bistable dynamics Boundary formation Cooperative interactions
Eukaryotic gene activity is strongly related to the epigenetic state [1–4]. Covalent histone modifications, chromatin folding into higher-order structures, and binding of large regulatory protein assemblies to the chromatin play a critical role in generating the epigenetic state. Changes in chromatin structure are guided by post-translational covalent histone modifications. Histones carry a multitude of posttranslational chemical modifications (for example, acetylation, methylation, ubiquitination, and phosphorylation) with dedicated patterns along the genome that correlate with defined gene expression patterns [2, 4, 5]. During differentiation, the large-scale chromatin structure and its histone modification composition defines whether the chromatin holds a permissive or restrictive chromatin state, which determines gene sensitivity to transcription factors; thereby establishing cell-type specific expression [2, 6].
Besides cell differentiation, covalent-histone modifications are proposed to mark transcription initiation, for example, involving ATP-dependent nucleosome remodeling, establishment of the pre-initiation complex, and allowing RNA polymerase II to start and proceed transcriptional elongation [7–14].
A prominent feature of core histones is the large amount and diversity of covalently modified residues they can possess [2, 15]. These marks guide the recruitment and binding of defined histone modifying enzymes and other regulatory proteins . The histone modification composition along chromatin is often bordered by insulators, which likely demarcate the end of a local histone modification pattern [17, 18]. All of the above mentioned processes play a key role in regulation of transcription at the level of single genes [19–22].
A large body of evidence indicates that transcription regulation is dynamic and that it is inherently stochastic at various levels to give rise to large cell-to-cell heterogeneity  in transcriptional activity and in the resulting number of transcripts per cell . The importance of the local epigenetic neighborhood of genes for stochastic transcription is becoming more apparent [21, 22, 25, 26]. The stochastic aspects of the establishment, maintenance, and decay of histone modification patterns are still largely unresolved. Some experiments indicate that the dynamics of these patterns may explain the large cell-to-cell variability of eukaryotic gene expression [22, 25, 27–29].
Insight into the mechanisms underlying the dynamics of histone modification pattern formation is largely lacking, mainly because many factors act in concert and simultaneous measurement or control of these factors is still experimentally challenging. In such a case, mathematical models can be effective to provide additional mechanistic insight into the molecular behavior of the system. Such models help to stimulate intuitive understanding and suggest new experiments to test the predicted behavior. Several computational efforts have enhanced our understanding of histone modification spreading and bistability [30–34]. For instance, the initial theoretical model of Dodd et al.  showed that the cooperativity of histones is essential to obtain bistability in distinct epigenetic states. However, the current body of computational models of histone modification patterning and spreading describes histone modification kinetics and cooperative processes using phenomenological process descriptions [30–36]. Although this is advantageous for reaching qualitative understanding of the dynamic modes of histone modification mechanisms, it does not allow us to predict the system consequences from basic biochemical and biophysical parameters. As a result, the exact time-scales of histone modification pattern formation and how protein-chromatin interactions and the movement of enzymes along the chromatin affect such patterns often remain unclear. Such knowledge is becoming more and more relevant as the field is moving in the direction of single-cell studies that address the influence of histone modification dynamics on gene regulation.
In this work, we study how the biochemical mechanisms of histone regulatory proteins give rise to spreading of histone modifications along the body of a gene . We focus on the recruitment of regulatory proteins to chromatin [38, 39], the diffusion of regulatory proteins along the chromatin [40–42], and gene-connectivity through chromatin-chromatin interactions. To achieve this, we present a novel stochastic model with physical model parameters and mechanisms for chromatin binding, 1D diffusion, and recruitment of regulatory proteins, and gene connectivity based on biochemical mechanisms and literature-based kinetic parameters. This model is implemented in a dedicated software package that can be run as a plugin of the stochastic simulation software StochPy . We study how this basic biochemical characterization of covalent histone modification turnover generates the establishment, maintenance, and decay of histone modification patterns along the body of a single gene. We find that histone modification patterns can be established within minutes by combined recruitment and diffusion mechanisms. A strong synergism between diffusion and weak recruitment leads to nearly irreversible transitions in histone modification patterns providing stable epigenetic patterns even in the absence of chromatin connectivity. In the absence of chromatin connectivity, spontaneous and dynamic histone modification boundaries can be formed but results in a highly unstable, fluctuation-driven state. Chromatin connectivity introduces bistability in the epigenetic state involving spontaneous, erratic switching between the two opposite (that is, a fully active and inactive) histone modification states of the model gene, and the formation of stable epigenetic histone modification patterns along the entire body of the gene. We demonstrate that such behavior can be found within biologically plausible parameter regimes. This study suggests an important role for epigenetics in gene activity regulation. The stochastic switching of the epigenetic state that we observe can contribute to cell-to-cell heterogeneity of transcription while the observed formation of a stable epigenetic pattern could play a role in robust gene silencing.
The histone modification spreading model
In the simulations, the modification process is initiated from a specific site in the nucleosomal array. This targeted nucleosome either corresponds to a nucleosome with a specific histone modification composition, or to a nucleosome-free region (for example, a 3’ or 5’-NFR) with an accessible protein-DNA binding domain. Such 5’-NFRs often demarcate the transcription start site. They are typically enriched in transcription factor binding sequences  and known to be involved in the induction of transcription [22, 45]. Here, we assume this initiation site or DNA binding site at the position of one designated nucleosome in the center of the array and we refer to this nucleosome as the initiation site (Figure 1A).First, we monitor only the methyltransferase explicitly and we assume the acetyltransferase as background activity. The modification mechanism of a single nucleosome is given in Figure 1B. The choice for the methyltransferase here is arbitrary, as the reverse situation - that is, an explicit acetyltransferase and background methyltransferase - would give identical results with the selected parameters. Therefore, conclusions that are drawn from our model simulations can be applied to either activation or silencing of the modeled gene. Later we relax this assumption and track the fate of both the acetyl- and the methyltransferase, and add ‘jumping’ of the transferase by means of chromatin connectivity.
We distinguish several mechanisms for transferase movement along the nucleosomal array: diffusion (Figure 1C), modification-induced recruitment (Figure 1D), and a combination of both (Figure 1E). Transferase diffusion is based on the findings that proteins can display 1D diffusion along DNA and chromatin [40, 46–49], which is implemented in the model as a 1D random walk. A diffusion-event is reflected in the simulations as a discrete step from one nucleosome to a neighboring nucleosome (Figure 1C). Note that when a transferase encounters another transferase on the neighboring nucleosome it cannot push or diffuse over this obstacle. Before a next diffusion step, the transferase can modify an encountered unmodified nucleosome. The next mechanism, modification-induced recruitment (Figure 1D), corresponds to a mechanism observed during the formation of heterochromatin induced by Heterochromatin Protein 1 (HP1) (reviewed by ). HP1 binds with its chromo domain to chromatin at histone H3 trimethylated at lysine 9 (H3K9me3) . With its chromo shadow domain HP1 can homodimerize with HP1 isoforms or heterodimerize with histone methyltransferase enzymes such as SUV39h1 (the mammalian ortholog of Drosophila suppressor of variegation 3-9 (Su(var)3-9)), thereby inducing spreading of H3K9Me3 . Here, we simplify the properties of HP1 and SUV39h1 (histone modification binding and methylation of neighboring histones) into properties of the single transferase in our model. Thus, a methyltransferase with modification-induced recruitment properties can bind at any methylated nucleosome, and it can methylate both the occupied nucleosome and neighboring nucleosomes. However, the methyltransferase will not diffuse between the nucleosomes. In this sense the mechanism is similar to previously published models [30, 31]. In a model with combined diffusion and recruitment properties (Figure 1E) the transferase can perform all of these activities: diffuse to a neighboring nucleosome, bind to a methylated nucleosome, methylate the bound nucleosome, and methylate a neighboring nucleosome. Regardless of the propagation mechanism, the transferase can dissociate from a nucleosome at any time and position.
Overview of the kinetic parameters in the model
Influx of transferases at initiation site
2.4 (one enzyme models)
0.01 (two enzyme models)
Release rate of transferase from nucleosome
Modification rate of nucleosome
Modification rate of neighboring nucleosome
1D diffusion rate over the chromatin
Influx of transferases at modified nucleosome
Rate constant of demodification
1D diffusion and recruitment mechanisms show localized modification patterns
We introduce four characteristics to analyze how modification patterns are established, maintained, and removed with the three mechanisms of transferase movement. The first characteristic is the stationary methylation and methyltransferase occupancy pattern along the array in steady state. Such a simulated stationary pattern corresponds to a histone modification pattern as observed with ChIP-seq analysis. The second characteristic is the activation dynamics. This characteristic is plotted as the methyltransferase occupancy on the nucleosomal array starting from the fully acetylated array state. As third characteristic we compute the steady-state probability distribution of the number of methyltransferases and the three modification states on the array. Finally, we test the relaxation dynamics. This fourth characteristic determines the stability of a modification pattern. To test this we simulate the system from an initial state in which all nucleosomes are methylated and bound to a methyltransferase, subsequently, we determine the decline in the number of methylated nucleosomes. In this fourth test we do not allow any influx of new methyltransferases at the initiation site and we assume that the initiation site is inaccessible, for instance due to nucleosomal repositioning.
Strong synergism between 1D diffusion and recruitment leads to bistable behavior
Introduction of two opposite modification enzymes causes unstable boundary formation
Simulations performed with a single transferase indicated that an elevated influx of transferases at the initiation site causes a nearly permanent occupation of the initiation site. This situation leads to stable modification peaks around each initiation site (Figure 3 column i). Because we are interested in the origin of bistability we study conditions in which one modification overtakes the other. Therefore, we lower association rates at the initiation site, and take a rate constant of 0.01 s-1 for the initiation rate (kinitiation).
The recruitment mechanism with RE = 0.5 (so a low initiation rate) yields a modification pattern close to the initiation sites (Figure 6C). This is in agreement with Figure 3B in which each of the transferases cannot modify the nucleosomes further away from their initiation site. When the transferases are simulated with a higher efficiency recruitment mechanism (RE = 2.0) the opposing transferases form a ‘boundary’ in the center of the nucleosomal array (Figure 6D). In this case, the transferases can reach nucleosomes far from their initiation sites but since they are not able to pass each other, a clear boundary between acetylation and methylation states is generated. This boundary is on average in the middle of the array but its location is unstable. Moreover, methylation and acetylation can outcompete each other to give rise to a nucleosomal array that is either totally methylated or acetylated. The initiation site of the methyltransferase or acetyltransferase can be blocked by the presence of an antagonistic transferase on the initiation site. New influx of transferases is hampered if the krecruitment is higher than kinitiation and as a result, the entire system remains mostly in one modification state. When an initiation event occurs, in principle, the modification that is outnumbered is able to overtake the opposite modification but the probability of this event is relatively low. This simulation, therefore, indicates that spontaneous boundary formation by opposing modifications may eventually lead to one modification dictating the opposite modification. In Figure 6E, an indication of the stability of the boundary position is given by the difference between total methylation and total acetylation on the array. Within the 700 s of simulated time the recruitment mechanism rarely causes the formation of states where the majority of the array is methylated or acetylated. The combined recruitment and diffusion mechanism with RE = 0.2 (Figure 6F) also shows this boundary formation although in a more dynamic fashion.
Chromatin connectivity induces bistable dynamics in modification pattern formation
Without chromatin connectivity, we hardly find any modifications on the nucleosomal array in the simulations of the opposing transferases with the recruitment mechanism (RE =0.5; Figure 6C). In contrast, upon the introduction of 10 chromatin interaction sites, each with an average interaction frequency of one interaction per 100 s, the system forms highly stochastic modification dynamics with large variability in array-wide methylation and acetylation levels (Figure 7Ci). Around the individual chromatin interaction sites the low activity of recruitment causes dynamic and local methylation and acetylation regions. A higher interaction frequency (10 interaction sites, frequency =0.1 s-1) increases the modification extent of the array causing spontaneous state switching between an almost completely methylated and acetylated nucleosomal array (Figure 7Cii).
The introduction of chromatin connectivity in the recruitment mechanism of Figure 6E destabilizes the acetylation-methylation modification boundary. At an intermediate chromatin interaction frequency (at rate 0.1 s-1) and with five interaction sites, chromatin connectivity allows acetyl- and methyltransferases to pass over each other and enables the nucleosomal array to switch between alternative states (Figure 7Di). High frequency interactions (10 interaction sites, frequency 0.1 s-1) make the state switching more rapid and induce a more stable chromatin composition (Figure 7Dii). In the combined mechanism (Figure 7E) of propagation and low recruitment, the effect of the chromatin connectivity is similar to the results found in Figure 7D. Intermediate connectivity enables faster switching between opposing modifications states and an increased interaction frequency causes switching between nucleosome array states that are almost completely methylated or acetylated (Figure 7Eii). Both in the recruitment mechanism and the combined mechanism, bistability is introduced through chromatin connectivity.
In this study we present a stochastic mathematical model that describes the mechanisms of histone modification spreading. This model incorporates current biological understanding of the processes involved in the spreading of histone modifications, such as diffusion and modification induced recruitment of methyl- and acetyltransferases. We include parameter values that match realistic values enabling the model to describe processes of both dynamic and stable epigenetic pattern formation on both local chromatin sites and larger more widespread chromatin regions.
Previous model efforts have described the dynamics of histone modification state switching. Dodd et al.  developed a model of the behavior of epigenetic state switching at the mating locus of Schizosaccharomyces pombe. In this model of Dodd et al., nucleosomes and their modifications are simplified to three nucleosomal states with state conversions that depend on positive feedback interactions with other nucleosomes in the nucleosomal array. They simplified the model kinetics providing only two possible reactions for a nucleosome: (i) active conversion, driven by positive feedback between two selected nucleosomes, and (ii) noisy conversion of a nucleosome state into a different state. The Dodd model indicates the importance of the feedback-to-noise ratio in the histone modification reactions and the necessity for cooperativity to establish a bistable epigenetic system. Angel et al.  implemented the Dodd approach to develop a model that explains the Arabidopsis epigenetic switch for flowering after an environmental alteration. This Arabidopsis model is composed of a two-compartment system in which the nucleation site is modeled as a more interactive region than the surrounding regions. The model allowed the authors to fit the experimentally measured spreading of histone H3 lysine 27 methylation (H3K27me) on the flowering locus. Satake and Iwasa  also used a model to simulate Arabidopsis during environmental changes and to explain stability of the switch in a stochastic environment. In their study, the rates of modification were formalized as reactions that are dependent on the total amount of modified histones. Thus, in this model a positive feedback-loop was implemented independent of the position of the histones or the mechanism of feedback regulation. In these previous models the molecular interactions of the regulatory proteins (transferases) producing histone modifications were not the object of interest and therefore not mechanically incorporated. Our model extends those studies and considers the kinetic descriptions of nucleosome modifications, that is, diffusive propagation and recruitment of modification enzymes, and chromatin connectivity. As a result, our simulations concern physically realistic parameters, concentrations of modifications enzymes, and time-scales of kinetic and diffusive processes.
All our investigated mechanisms involving transferase propagation (that is, separate or combined diffusion and recruitment as well as chromatin connectivity induced propagation) show that stable patterns of histone modifications can be formed in a dynamic environment. The different mechanisms result in different types of modification patterns. For instance, the diffusion mechanism produces localized patterns around a DNA binding site. These patterns resemble the enrichment of transcription factors commonly found at promoter and enhancer sites of active genes [2, 9]. Additionally, the modification induced recruitment mechanism with a low recruitment-efficiency gives rise to patterns that are similar to patterns generated by the diffusion mechanism. However when the recruitment efficiency is increased, the patterns are widespread and stable. This increased recruitment efficiency corresponds with measured lifetimes of repressive marks that are on average longer than activating histone modifications , since long-lived modifications increase the RE of the modification. These widespread patterns are also observed when the mechanisms for diffusion and recruitment are combined together even at low recruitment efficiencies. The latter patterns are reminiscent of histone modifications commonly found at genomic regions consisting of transcriptionally silent genes in which widespread patches of defined histone modification patterns such as H3K9me3 or H3K27me are found [2, 5]. Remarkably, the dynamics of a single transferase in our model is fast, while the stability of the pattern is large. This corresponds with the findings of Cheutin et al.  showing that even though HP1 exhibits relatively fast dynamic binding kinetics it is involved in creating stable modification patterns.
An important factor for the stability of global modification patterns in the model is the recruitment efficiency. Here, we define this efficiency in terms of two parameters, the lifetime of the modification and the time before recruitment (1/kdemodification and 1/krecruitment). A change in these parameters can cause large differences in the stability of the modification patterns on the chromatin. Biological alterations in these parameters can be caused by changes in the concentrations of the involved enzymes or changes in their binding affinities. In addition, a change in the binding affinity could be caused by, for example, the presence of another histone modification at the same nucleosome, which could influence the stability of the interaction of the transferase with the nucleosome. Furthermore, we showed that the combined diffusion and recruitment mechanism has a high capability of initiation from a small amount of seeds. These experiments were shown in a more or less neutral background, in which acetylation is not explicitly modeled. This means that although a high seeding capacity is shown, the effects of a heterochromatin protein in an actively acetylated gene could have a different outcome. However, in a favorable environment a fast switch from active to inactive chromatin can be established by the activity of only a few methyltransferases.
Surprisingly, the rate of spreading that we observe in the widespread modification patterns is faster than we expected. In our study, we observe single gene level histone modification pattern formation at a minute scale; around 1 min for the combined mechanism (Figure 3C) and within 5 min in the case of high efficiency recruitment (Figure 4A). Our obtained time scales agree with diffusion times of molecules over chromatin [40, 41]. It should be noted that the mechanics used for this study are simplified in such manner that all protein complex formation and successive binding reactions are not explicitly included. Therefore, the model presented here shows the high end of the possible spreading rates. However, incorporation of complex formation or subsequent binding effects will not cause our model timescale to change from seconds to hours. Hathaway et al.  developed a model to mechanistically simulate histone modification spreading. They assumed a straightforward neighbor-neighbor interaction of the nucleosomes. With the model assumptions that each nucleosome can modify its neighbor at a certain rate and that each histone can lose its modification at another rate. This model specifically creates a way to study histone modification spreading and pattern formation. The spreading kinetics noticed by Hathaway et al. is fitted to their ChIP data of histone modification patterns and provides insight into the average rate of H3K9me3 spreading over the chromatin. The follow-up of this article by Hodges et al.  showed that the extent of histone modification spreading was restricted intrinsically by the properties of the spreading mechanism. The measurements in the Hathaway  article are based on cell population dynamics of histone modifications measured in ChIP experiments, whereas our simulations concern histone modification pattern formation at the level of a single gene. Our study provides insight in histone modification spreading along single genes. The experimental testing of this single gene behavior is an important next step to confirm its biological relevance.
The simulations with two opposing transferases show the capacity of the system to form a natural boundary between an epigenetic active and an inactive stretch of nucleosomes. This is, however, a boundary with an unstable position, indicating that additional control factors are necessary to stabilize the position of the boundary. Two model assumptions are decisive for this (unstable) boundary formation: (i) opposing transferases cannot pass each other on the chromatin causing confinement of the transferases, and (ii) the acetyltransferase and the methyltransferase have the same kinetic properties, which generates a balanced competition of the transferases for the nucleosomes. Interestingly, the system as we show it would be a plausible explanation for early stages in the differential expression of genes that are affected by Position Effect Variegation (PEV). In PEV, genes that are close to an inactive region of the chromatin (mostly telomeric or centromeric) can be differentially expressed in a subpopulation of cells dependent on the dynamic position of the boundary of the heterochromatin region [11, 37].
Similar to computational efforts on histone modification spreading introduced by other authors , we find that chromatin connectivity is a mechanism for bistable pattern formation. We show that the introduction of connectivity is necessary for the transferases to pass the boundary and create stretches of fully methylated or fully acetylated chromatin. In nucleosome arrays with two or three interaction sites this process is still quite inefficient but with more interactions, five or 10, and a high interaction-frequency the stability of an active or silenced chromatin state is increased, but also the switch process from one state to the other is faster. During embryonic development gene-expression patterns occur that are cell type specific. Bistability could represent advantageous system behavior, enabling differentiating cells to switch the expression activity of several genes by allowing the chromatin to change its conformation rapidly over a long length scale.
Our study shows that histone modification patterns can be established within minutes and that nearly irreversible transitions can result, which provides stable epigenetic patterns. Introduction of two opposing enzymes causes dynamic histone modification boundaries whereas chromatin connectivity can introduce both bistable epigenetic state switching and stable histone modification patterns. We demonstrate in our model that with biologically plausible parameter regimes both epigenetic stochastic switching and stable pattern formation are noticed to provide cell-to-cell heterogeneity and robustness in gene expression.
Here is the apparent first-order association rate constant for the binding of modification enzymes to a DNA site, kon is the second-order diffusion limited rate constant, nE is the number of modification enzymes per nucleus with volume Vnucleus, DE and DS are the diffusion coefficients of the enzyme and the DNA site, and rE and rS are the radii of the enzyme and the DNA site. Since we fix the total number of modification enzymes the relevant rate constant to consider is .
The diffusion constant of enzymes (DE) in a cell are between 0.5 and 5 μm2s-1 . The diffusion constant of the DNA binding site (D S ) is expected to be a few orders of magnitude lower than the diffusion constant of the enzyme, therefore we neglect this constant. We assume a transferase of about 5 nm in diameter and a DNA binding site of 15 bp, which measures 5 nm. The volume of an average human nucleus is 500 pL. The concentration of specific histone modification transferases in the nucleus has not been reported so far. Transcription factors are estimated to be present in a wide copy range from tens to thousands per cell . Here we assume that 5,000 copies of modification and demodification enzymes are present per cell. Together these values give an apparent diffusion limited association rate of 2.4 s-1. This value is used for the binding at the initiation site, the rate of demodification and the rate of modification for the acetyltransferase in the first section of the results (parameters kon, kdemodification, and kmodification).
The transferase nucleosomal release rate is based on FRAP analysis of transcription factors and histone modification transferases . The residence time of different factors ranges from milliseconds to several minutes. HP1β is shown to have a chromatin binding residence time of 11.3 s in the fast fraction, covering 88% of the HP1β population. We consider that the transferase has an average residence time of 10 s on the nucleosome array. Thus, the parameter for the transferase detaching from the nucleosome (koff) is 0.1 s-1.
Gorman et al. [40, 41] studied 1D sliding of proteins on the DNA. Several different proteins have been studied in an in vitro setting and their diffusion coefficients have been measured. Measured diffusion coefficients range from 2 × 10-4 μm2s-1 to 0.5 μm2s-1. For our model we assume the lower end of the range because most of the measurements were performed on stretched, nucleosome free DNA while the measurements on chromatin gave varying results, depending on the measured transferase . The diffusion parameter is defined as the rate of movement over the approximately 25 nm distance from nucleosome to nucleosome, resulting in a diffusion parameter (kslide) of 0.6 s-1.
We assume that the transferase modification reaction is fast once the transferase is bound to the nucleosome (ktransferase =1,000 s-1). The transferase neighbor modification reactions however are much slower due to the persistence length of the DNA over short distances (kneighbor =0.2 s-1). From these assumptions we calculated the model parameters (Table 1). Although we made assumptions in the parameter set, our model is able to describe the dynamics of the system based on biological data in a qualitative manner. A description for the NucleosomeModelBuilder, that is part of the NucleosomeTool, is included in Additional file 1: S1. The reactions used to describe the model are summarized in Additional file 10: S5.
StochPy data analysis and probability distributions
The model with different mechanisms and parameter settings was simulated with the NucleosomeTool plug-in designed for StochPy . StochPy uses the Gillespie algorithm for the stochastic simulations. We selected the next reaction method for all simulations in this study. StochPy and the NucleosomeTool are available from: http://stochpy.sourceforge.net/. Simulation data are saved and imported into Mathematica (Wolfram Research) for further analysis.Simulations with one transferase are run for at least 1,000 time trajectories of at least 400 s simulated time. This accounts for 34 to 89 million data points per simulated mechanism. Dynamics plots (Figure 3 column ii and iv) are calculated by binning of data into 1 s bins and then calculation of lowest value, highest value, median, and 0.25 and 0.75 quartile range of the data points per bin. These are plotted as error bars to the median value.To create the plots of steady-state probabilities of the model species the simulations are repeated 1,000 times with an initial state taken from an average steady-state situation. Pattern distributions (Figure 3 column i) are calculated from the final 60% of the steady-state data to ensure steady state are reached. The nucleosomes in each state are multiplied by the time spent in this state divided by the total time to calculate the probability of each state at each nucleosome position. Probability distribution (Figure 3 column iii) of the steady-state simulations is obtained by binning the total number of nucleosomes in each state (M, A, or U) and calculating the time spent in each total amount divided by the total time. Of the steady-state simulations only the final 60% of the simulated data are used to ensure steady-state values are calculated.
In the simulations with two opposing enzymes, lowering the initiation parameter from 2.4 s-1 to 0.01 s-1 created the opportunity for one transferase to completely occupy the array. Because the computational limits of our system are reached at approximately 600 to 800 model seconds the simulations are run for 500 s. Figures 6 and 7 contain simulations of 500 s that are placed adjacently to give an illustration of long-term model behavior.
PJV (Associate Professor and Group leader) investigates functional and dynamic organization of the mammalian genome exploring the use of engineering (synthetic) mammalian cell systems allowing to modulate the epigenetic state performing quantitative single-cell measurements and computational stochastic modeling. Her main research focus is on epigenetic state switching and transcriptional dynamics/stochasticity also in the context of deregulated behavior such as UV-induced DNA damage repair and age-related diseases.
FJB (Associate and Extraordinary Professor) uses theory and experiments to study the mechanisms and principles of cellular adaptation to new conditions developing simulation software for combined simulation of stochasticity of molecular reactions and cell growth. Main research topics include metabolic regulation by gene expression, mixing of transcription stochasticity and cell-growth stochasticity, and genome-wide adaptation of metabolism.
LCMA and TRM are PhD students using systems biology approaches to understand complex biological behavior.
Acetylated or active nucleosome state
Chromatin ImmunoPrecipitation followed by sequencing of the obtained DNA fragments
Histone H3 lysine 27 methylation
Histone H3 trimethylated at lysine 9 HP1, Heterochromatin Protein 1
Methylated or silenced nucleosome state
Position effect variegation
The mammalian ortholog of Drosophila suppressor of variegation 3-9 (Su(var)3-9)
Unmodified or neutral nucleosome state.
This work was supported in part by grants from the Netherlands Organization for Scientific Research (ALW-NWO-Meervoud 836.07.001 to PJV, NWO-VIDI 864.11.011 to FJB) and from the Netherlands Institute for Systems Biology (NISB) (Glue-project 2008-III) and the Centre for Mathematics and Computer Science (CWI). We thank Mannus Kempe, Prof. HW Westerhoff and the other members of the SSB-NOG group (SILS, UvA, Amsterdam, the Netherlands) and Prof. AB Houtsmuller (ErasmusMC, Rotterdam, The Netherlands) for critical discussions. We thank SURFsara (http://www.surfsara.nl) for the support in using the Lisa Compute Cluster.
- Kouzarides T: Chromatin modifications and their function. Cell. 2007, 128: 693-705. 10.1016/j.cell.2007.02.005.View ArticlePubMedGoogle Scholar
- Barski A, Cuddapah S, Cui K, Roh T-Y, Schones DE, Wang Z, Wei G, Chepelev I, Zhao K: High-resolution profiling of histone methylations in the human genome. Cell. 2007, 129: 823-837. 10.1016/j.cell.2007.05.009.View ArticlePubMedGoogle Scholar
- Turner BM: The adjustable nucleosome: an epigenetic signaling module. Trends Genet. 2012, 28: 436-444. 10.1016/j.tig.2012.04.003.View ArticlePubMedGoogle Scholar
- Karlić R, Chung H-R, Lasserre J, Vlahovicek K, Vingron M: Histone modification levels are predictive for gene expression. Proc Natl Acad Sci U S A. 2010, 107: 2926-2931. 10.1073/pnas.0909344107.PubMed CentralView ArticlePubMedGoogle Scholar
- Ernst J, Kheradpour P, Mikkelsen TS, Shoresh N, Ward LD, Epstein CB, Zhang X, Wang L, Issner R, Coyne M, Ku M, Durham T, Kellis M, Bernstein BE: Mapping and analysis of chromatin state dynamics in nine human cell types. Nature. 2011, 473: 43-49. 10.1038/nature09906.PubMed CentralView ArticlePubMedGoogle Scholar
- Cosgrove MS, Wolberger C: How does the histone code work?. Biochem Cell Biol. 2005, 83: 468-476. 10.1139/o05-137.View ArticlePubMedGoogle Scholar
- Barth TK, Imhof A: Fast signals and slow marks: the dynamics of histone modifications. Trends Biochem Sci. 2010, 35: 618-626. 10.1016/j.tibs.2010.05.006.View ArticlePubMedGoogle Scholar
- Lee TI, Jenner RG, Boyer LA, Guenther MG, Levine SS, Kumar RM, Chevalier B, Johnstone SE, Cole MF, Isono K, Koseki H, Fuchikami T, Abe K, Murray HL, Zucker JP, Yuan B, Bell GW, Herbolsheimer E, Hannett NM, Sun K, Odom DT, Otte AP, Volkert TL, Bartel DP, Melton DA, Gifford DK, Jaenisch R, Young RA: Control of developmental regulators by Polycomb in human embryonic stem cells. Cell. 2006, 125: 301-313. 10.1016/j.cell.2006.02.043.PubMed CentralView ArticlePubMedGoogle Scholar
- Li B, Carey M, Workman JL: The role of chromatin during transcription. Cell. 2007, 128: 707-719. 10.1016/j.cell.2007.01.015.View ArticlePubMedGoogle Scholar
- Berger SL: The complex language of chromatin regulation during transcription. Nature. 2007, 447: 407-412. 10.1038/nature05915.View ArticlePubMedGoogle Scholar
- Grewal SIS, Jia S: Heterochromatin revisited. Nat Rev Genet. 2007, 8: 35-46. 10.1038/nrg2008.View ArticlePubMedGoogle Scholar
- Weake VM, Workman JL: Inducible gene expression: diverse regulatory mechanisms. Nat Rev Genet. 2010, 11: 426-437. 10.1038/nrg2781.View ArticlePubMedGoogle Scholar
- Eissenberg J, Shilatifard A: Leaving a mark: the many footprints of the elongating RNA polymerase II. Curr Opin Genet Dev. 2006, 16: 184-190. 10.1016/j.gde.2006.02.004.View ArticlePubMedGoogle Scholar
- Saunders A, Core LJ, Lis JT: Breaking barriers to transcription elongation. Nat Rev Mol Cell Biol. 2006, 7: 557-567. 10.1038/nrm1981.View ArticlePubMedGoogle Scholar
- Cosgrove MS, Boeke JD, Wolberger C: Regulated nucleosome mobility and the histone code. Nat Struct Mol Biol. 2004, 11: 1037-1043. 10.1038/nsmb851.View ArticlePubMedGoogle Scholar
- Suganuma T, Workman JL: Signals and combinatorial functions of histone modifications. Annu Rev Biochem. 2011, 80: 473-499. 10.1146/annurev-biochem-061809-175347.View ArticlePubMedGoogle Scholar
- Gaszner M, Felsenfeld G: Insulators: exploiting transcriptional and epigenetic mechanisms. Nat Rev Genet. 2006, 7: 703-713.View ArticlePubMedGoogle Scholar
- Grewal SI, Elgin SC: Heterochromatin: new possibilities for the inheritance of structure. Curr Opin Genet Dev. 2002, 12: 178-187. 10.1016/S0959-437X(02)00284-8.View ArticlePubMedGoogle Scholar
- Bryant GO, Prabhu V, Floer M, Wang X, Spagna D, Schreiber D, Ptashne M: Activator control of nucleosome occupancy in activation and repression of transcription. PLoS Biol. 2008, 6: 2928-2939.View ArticlePubMedGoogle Scholar
- Heintzman ND, Stuart RK, Hon G, Fu Y, Ching CW, Hawkins RD, Barrera LO, Van Calcar S, Qu C, Ching KA, Wang W, Weng Z, Green RD, Crawford GE, Ren B: Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat Genet. 2007, 39: 311-318. 10.1038/ng1966.View ArticlePubMedGoogle Scholar
- Mao C, Brown CR, Falkovskaia E, Dong S, Hrabeta-Robinson E, Wenger L, Boeger H: Quantitative analysis of the transcription control mechanism. Mol Syst Biol. 2010, 6: 431.PubMed CentralView ArticlePubMedGoogle Scholar
- Kim HD, O’Shea EK: A quantitative model of transcription factor-activated gene expression. Nat Struct Mol Biol. 2008, 15: 1192-1198. 10.1038/nsmb.1500.PubMed CentralView ArticlePubMedGoogle Scholar
- Hager GL, McNally JG, Misteli T: Transcription dynamics. Mol Cell. 2009, 35: 741-753. 10.1016/j.molcel.2009.09.005.View ArticlePubMedGoogle Scholar
- Raj A, Peskin CS, Tranchina D, Vargas DY, Tyagi S: Stochastic mRNA synthesis in mammalian cells. PLoS Biol. 2006, 4: e309-10.1371/journal.pbio.0040309.PubMed CentralView ArticlePubMedGoogle Scholar
- Raj A, van Oudenaarden A: Single-molecule approaches to stochastic gene expression. Annu Rev Biophys. 2009, 38: 255-270. 10.1146/annurev.biophys.37.032807.125928.PubMed CentralView ArticlePubMedGoogle Scholar
- Blake WJ, Balázsi G, Kohanski MA, Isaacs FJ, Murphy KF, Kuang Y, Cantor CR, Walt DR, Collins JJ: Phenotypic consequences of promoter-mediated transcriptional noise. Mol Cell. 2006, 24: 853-865. 10.1016/j.molcel.2006.11.003.View ArticlePubMedGoogle Scholar
- Chubb JR, Trcek T, Shenoy SM, Singer RH: Transcriptional pulsing of a developmental gene. Curr Biol. 2006, 16: 1018-1025. 10.1016/j.cub.2006.03.092.View ArticlePubMedGoogle Scholar
- Raser JM, O’Shea EK: Noise in gene expression: origins, consequences, and control. Science. 2005, 309: 2010-2013. 10.1126/science.1105891.PubMed CentralView ArticlePubMedGoogle Scholar
- Mogno I, Vallania F, Mitra RD, Cohen BA: TATA is a modular component of synthetic promoters. Genome Res. 2010, 20: 1391-1397. 10.1101/gr.106732.110.PubMed CentralView ArticlePubMedGoogle Scholar
- Dodd IB, Micheelsen MA, Sneppen K, Thon G: Theoretical analysis of epigenetic cell memory by nucleosome modification. Cell. 2007, 129: 813-822. 10.1016/j.cell.2007.02.053.View ArticlePubMedGoogle Scholar
- Hathaway NA, Bell O, Hodges C, Miller EL, Neel DS, Crabtree GR: Dynamics and memory of heterochromatin in living cells. Cell. 2012, 149: 1447-1460. 10.1016/j.cell.2012.03.052.PubMed CentralView ArticlePubMedGoogle Scholar
- Hodges C, Crabtree GR: Dynamics of inherently bounded histone modification domains. Proc Natl Acad Sci U S A. 2012, 109: 13296-13301. 10.1073/pnas.1211172109.PubMed CentralView ArticlePubMedGoogle Scholar
- Satake A, Iwasa Y: A stochastic model of chromatin modification: cell population coding of winter memory in plants. J Theor Biol. 2012, 302: 6-17.View ArticlePubMedGoogle Scholar
- Sedighi M, Sengupta A: Epigenetic chromatin silencing: bistability and front propagation. Phys Biol. 2007, 4: 246-255. 10.1088/1478-3975/4/4/002.PubMed CentralView ArticlePubMedGoogle Scholar
- Angel A, Song J, Dean C, Howard M: A Polycomb-based switch underlying quantitative epigenetic memory. Nature. 2011, 476: 105-108. 10.1038/nature10241.View ArticlePubMedGoogle Scholar
- Mukhopadhyay S, Sengupta AM: The role of multiple marks in epigenetic silencing and the emergence of a stable bivalent chromatin state. PLoS Comput Biol. 2013, 9: e1003121-10.1371/journal.pcbi.1003121.PubMed CentralView ArticlePubMedGoogle Scholar
- Talbert PB, Henikoff S: Spreading of silent chromatin: inaction at a distance. Nat Rev Genet. 2006, 7: 793-803.View ArticlePubMedGoogle Scholar
- Lachner M, O’Carroll D, Rea S, Mechtler K, Jenuwein T: Methylation of histone H3 lysine 9 creates a binding site for HP1 proteins. Nature. 2001, 410: 116-120. 10.1038/35065132.View ArticlePubMedGoogle Scholar
- Canzio D, Chang EY, Shankar S, Kuchenbecker KM, Simon MD, Madhani HD, Narlikar GJ, Al-Sady B: Chromodomain-mediated oligomerization of HP1 suggests a nucleosome-bridging mechanism for heterochromatin assembly. Mol Cell. 2011, 41: 67-81. 10.1016/j.molcel.2010.12.016.PubMed CentralView ArticlePubMedGoogle Scholar
- Gorman J, Greene EC: Visualizing one-dimensional diffusion of proteins along DNA. Nat Struct Mol Biol. 2008, 15: 768-774. 10.1038/nsmb.1441.View ArticlePubMedGoogle Scholar
- Gorman J, Plys AJ, Visnapuu M-L, Alani E, Greene EC: Visualizing one-dimensional diffusion of eukaryotic DNA repair factors along a chromatin lattice. Nat Struct Mol Biol. 2010, 17: 932-938. 10.1038/nsmb.1858.PubMed CentralView ArticlePubMedGoogle Scholar
- Stanford NP, Szczelkun MD, Marko JF, Halford SE: One- and three-dimensional pathways for proteins to reach specific DNA sites. EMBO J. 2000, 19: 6546-6557. 10.1093/emboj/19.23.6546.PubMed CentralView ArticlePubMedGoogle Scholar
- Maarleveld TR, Olivier BG, Bruggeman FJ: StochPy: a comprehensive, user-friendly tool for simulating stochastic biological processes. PLoS One. 2013, 8: e79345-10.1371/journal.pone.0079345.PubMed CentralView ArticlePubMedGoogle Scholar
- Lee W, Tillo D, Bray N, Morse RH, Davis RW, Hughes TR, Nislow C: A high-resolution atlas of nucleosome occupancy in yeast. Nat Genet. 2007, 39: 1235-1244. 10.1038/ng2117.View ArticlePubMedGoogle Scholar
- Lam F, Steger D, O’Shea E: Chromatin decouples promoter threshold from dynamic range. Nature. 2008, 453: 246-250. 10.1038/nature06867.PubMed CentralView ArticlePubMedGoogle Scholar
- Blainey PC, Luo G, Kou SC, Mangel WF, Verdine GL, Bagchi B, Xie XS: Nonspecifically bound proteins spin while diffusing along DNA. Nat Struct Mol Biol. 2009, 16: 1224-1229. 10.1038/nsmb.1716.PubMed CentralView ArticlePubMedGoogle Scholar
- Tafvizi A, Mirny LA, van Oijen AM: Dancing on DNA: kinetic aspects of search processes on DNA. ChemPhysChem. 2011, 12: 1481-1489. 10.1002/cphc.201100112.PubMed CentralView ArticlePubMedGoogle Scholar
- Mirny L, Slutsky M, Wunderlich Z, Tafvizi A, Leith J, Kosmrlj A: How a protein searches for its site on DNA: the mechanism of facilitated diffusion. J Phys A Math Theor. 2009, 42: 434013-10.1088/1751-8113/42/43/434013.View ArticleGoogle Scholar
- Bonnet I, Biebricher A, Porté P-L, Loverdo C, Bénichou O, Voituriez R, Escudé C, Wende W, Pingoud A, Desbiolles P: Sliding and jumping of single EcoRV restriction enzymes on non-cognate DNA. Nucleic Acids Res. 2008, 36: 4118-4127. 10.1093/nar/gkn376.PubMed CentralView ArticlePubMedGoogle Scholar
- Maison C, Almouzni G: HP1 and the dynamics of heterochromatin maintenance. Nat Rev Mol Cell Biol. 2004, 5: 296-304. 10.1038/nrm1355.View ArticlePubMedGoogle Scholar
- Bohn M, Heermann DW: Diffusion-driven looping provides a consistent framework for chromatin organization. PLoS One. 2010, 5: e12218-10.1371/journal.pone.0012218.PubMed CentralView ArticlePubMedGoogle Scholar
- Zee BM, Levin RS, Xu B, LeRoy G, Wingreen NS, Garcia BA: In vivo residue-specific histone methylation dynamics. J Biol Chem. 2010, 285: 3341-3350. 10.1074/jbc.M109.063784.PubMed CentralView ArticlePubMedGoogle Scholar
- Cheutin T, McNairn AJ, Jenuwein T, Gilbert DM, Singh PB, Misteli T: Maintenance of stable heterochromatin domains by dynamic HP1 binding. Science. 2003, 299: 721-725. 10.1126/science.1078572.View ArticlePubMedGoogle Scholar
- Gorski SA, Dundr M, Misteli T: The road much traveled: trafficking in the cell nucleus. Curr Opin Cell Biol. 2006, 18: 284-290. 10.1016/j.ceb.2006.03.002.View ArticlePubMedGoogle Scholar
- Phair R, Scaffidi P, Elbi C, Vecerová J, Dey A, Ozato K, Brown DT, Hager G, Bustin M, Misteli T: Global nature of dynamic protein-chromatin interactions in vivo: three-dimensional genome scanning and dynamic interaction networks of chromatin proteins. Mol Cell Biol. 2004, 24: 6393-6402. 10.1128/MCB.24.14.6393-6402.2004.PubMed CentralView ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.