1. Introduction

Chromatin is the scene of a vast set of epigenetic modifications such as histone marks and DNA methylation. The epigenomic state of chromatin is locally defined along the genome, and this epigenomic state contributes to transcription regulation by modulating regulatory sequence accessibility, transcription factor binding, and the propensity of other genetic processes [1]. Recent initiatives successfully determined complete epigenomes for various human cell lines [2,3,4]. The consensual understanding of the functional role of epigenomic domains is that, in the course of development, [5,6,7], genomic loci are provided with various combinations of histone marks so as to selectively enable tissue-specific profiles of gene expression.

In parallel with the characterization of epigenomic domains, Chromosome Conformation Capture techniques (3C and HiC [8,9]) showed that chromosomes exhibit complex spatial architectures. In human cells, chromosomes are spatially segregated in territories which are themselves composed of topologically associating domains (TAD) at the kilobase-to-megabase scale [10,11]. TADs are consecutive sequences (from a few dozen to hundreds of kilobase pairs) of chromatin that preferentially form contacts with each other, and do not colocalize with sequences from other loci. Chromatin topology is linked to transcription control, for instance by selectively enabling regulatory sequence interaction with their target locus [12,13,14], by synchronizing expressions of genes located on the same TAD, thus acting as regulons [15], or by selectively localizing genes to transcription factories [16].

Although experimental data accumulates regarding both epigenomes and chromosome architecture thanks to advances in molecular biology and bioinformatics, a concise, clear-cut understanding of how epigenetic marks contribute to gene expression regulation through chromatin spatial reorganization remains elusive. The vast combinatoric complexity of epigenetic marks led to simplified models of epigenomic states where epigenetic marks are combined into a reduced set of so-called chromatin “colors” [4,17,18]. In this view, the genome is composed of a succession of blocks with various colors, with a typical size of a few hundred kilobase pairs. This also revealed that TADs and epigenomic color domains seem to roughly overlap, suggesting that epigenomic domains not only modulate transcription by acting as signalling “tags”, but also contribute to specify the spatial architecture of chromatin.

The modelling community basing the research on polymer physics seized this perspective to propose physically-driven models of chromosome conformation that would explain how epigenomic domains determine the architecture of the genome. These modelling studies are all based on copolymer physics principles, since chromosomes are copolymers and not just homopolymers, precisely because of epigenetic marking [19]. In this respect, there are two main modelling approaches: the String & Binders Switch (SBS) models, where chromosomes are flexible polymer combined with molecular actors that specifically bind chromatin loci together to form loops and other structures [20,21] ; and the interacting copolymer models in which monomers of chromatin have specific interaction energies according to their epigenomic state [22,23]. In [23], Jost et al. noticed that the overlap between TADs and epigenomic domains was strongly akin to block copolymers, and therefore introduced an interacting 2-color block copolymer model which provided a surprisingly rich conformational space and fairly reproduced contact maps observed experimentally using HiC assays. This suggests that the wealth of genetic expression patterns observed in human cell types is made possible by the richness of the chromatin conformational space, where combinations of TADs can be selectively opened/closed and/or expressed based on their conformation.

In parallel, we demonstrated recently that homogeneous polymers of different sizes would undergo a coil-globule transition at different temperatures [24]: for a given temperature, small polymers tend to adopt a coil conformation at equilibrium whereas larger polymers will adopt a globule conformation. In this study, we aim at demonstrating that this size-dependent coil-globule transition also affects the conformation of epigenomic domains and TADs. We show that the expression for the probability density of polymer conformations at equilibrium we proposed in [24], which takes into account the polymer size and its interaction energy parameters, also correctly describes the conformations observed when simulating a block copolymer with various block sizes. For a 2-color block copolymer, the conformation of epigenomic blocks is size-dependent and obeys the same coil-globule transition as homoplymers. Our main conclusion is that not only the strength of the interaction of epigenomic domains (both intra-domain and inter-domain interactions) affects chromatin conformations, but also the size of epigenomic domains physically regulates gene expression. This opens a new perspective in which simply changing the span of epigenomic domains (i.e. by spreading / erasing epigenomic marks) is sufficient to produce a vast set of selective chromatin folding patterns, hence a vast set of gene expression patterns.

2. Material and Methods

We present in this section the theoretical model and computational method we used for the finitesize scaling analysis of the conformational space of epigenomic domains. We explain the order parameters we retained for the theoretical study of the finite-size coil-globule transition for a homogeneous (single color) polymer, and the expected behaviors when varying the polymer size. We briefly introduce the block copolymer model akin to the one presented in [23] that we used for modelling chromatin domains. For the sake of simplicity, we chose a 2-color model (domains of A and domains of B). We then present the Langevin dynamics simulation software used for our numerical results.

2.1. Chromatin as a block copolymer

Our framework uses the same fundamental assumptions as presented in [23] : a chromatin fiber is modelled as a bead-spring chain (see fig.1). Each bead is assigned one color (either A or B). The color c_i of monomer i determines its interaction energy U_ij with monomer j of color c_j as follows:

Figure 1. Copolymer model. (A) Snapshot of a simulation showing an intermediate coilglobule conformation. monomers of type A are in red, monomers of type B are in green. (B) Principles of the bead-spring copolymer simulation model. Beads are linked by harmonic springs with stiffness k_s. Two monomers of type A (respectively B) have an interaction energy U_A (respectively U_B). Monomers of different colors do not interact except by volume exclusion. (C) The monomer sequences of simulated copolymers were chosen so that they are composed of 3 tandem repeats of a sub-sequence (AN_ABN_B) with N_A + N_B = 60.

DownLoad: Full-Size ImgPowerPoint

The interaction model is a cut-off Morse potential whose depth was set to U_s, with maximum range r_max, potential well width a and minimal energy distance r₀ (see table 1). U_s is the specific interaction energy between monomers of the same color, with U_S = U_A (respectively U_B) if c_i = c_j = A (respectively B) in Eq.1. We neglect non-specific interactions between monomers. The specific interaction energy represents the effect of epigenomic marks on chromatin monomer physico-chemical properties and more generally on their propensity to bind through other, possibly unknown, molecular mechanisms (insulator binding, loop-mediating proteins, etc.). Each bead represents 10kb of chromatin. The stiffness of the connecting springs is set to $k_s = \frac{3k_{B}T}{b^2}$, where k_B is Boltzmann’s constant, T the temperature and b is both the effective Kuhn length of the polymer and the spring relaxed extension.

Epigenomic domains are introduced by setting the color of monomers along the polymer in the following way : an alternance of N_A monomers of type A and N_B monomers of type B, repeated N_R times (denoted (A_NAB_NB)_NR ). Unless otherwise mentioned, N_A+N_B = 60 and N_R = 3 for all simulations. To simulate a homopolymer, we set N_A = N; N_B = 0; N_R = 1. We will call N the polymer size of which the conformation is observed, so that for homopolymers, N is the total polymer size, and for block copolymers, N is the size of the observed epigenomic domain (N_A or N_B respectively). The premise of the model is therefore similar to [23], except that we use a global-scheme Langevin dynamics engine (detailed further) which allows us to efficiently sample the conformational space, and we vary the domain sizes as well as the interaction energy to elucidate the size dependence of conformations.

2.2. Finite-size scaling of the coil-globule transition

A flexible homopolymer composed of N identical monomers can adopt a conformation at equilibrium that is either coil (self avoiding random walk) or globule (compact random walk) [25]. The polymer can transit from a coil to a globule as the temperature decreases. The typical physical situation where such a transition is observed occurs for a polymer in a bad solvent : as the temperature decreases, the interaction energy between monomers favors a more compact, dense conformation. In our study, the temperature is kept fixed at T = 300K and the interaction energy U_S is varied, which achieves the same effect. The same transition from coil to globule is therefore expected as the interaction energy increases at fixed temperature.

The conformation of the chain can be characterized by the probability density of the radius of gyration r at equilibrium:

where r_i is the position of the ith monomer and R_G is the position of the center of mass of the polymer. At a given chain size N, the larger r, the more relaxed the conformation, and the smaller r, the more globular, the denser the conformation. A polymer of size N can have one of three principal types of conformations, each having its specific scaling : globules for which < r > $\propto N^{1/d}$ (d is the spatial dimension), coils for which < r > $\propto N^\nu$ ($\nu$ is the Flory exponent), and stretchs for which < r > $\propto N$. These scalings are however only valid in the thermodynamic limit ($N \to \infty$). In previous works [26,24], we proposed an expression for the finite-size Boltzmann-Gibbs distribution of the chain conformations by combining these three scalings into a single distribution by means of finitesize scaling analysis:

The coefficients a_i depend on the monomer interaction energy U_S and the temperature T and are related to the global shape of the distribution.

The complete derivation of Eq.3 is available in [26]. For the sake of clarity and completeness, we give here the basic elements of this derivation. First, Eq.3 is the probability density of the order parameter $\hat{t}$, which is a proxy for the polymer conformation derived from the radius of gyration r and the density $\rho = N/r^{3}$ :

Hence, the larger $\hat{t}$, the denser the conformation. Then, the two main classes of conformations in the thermodynamic limit (coils and globules) appear in Eq.3 as follows. The globule limit of the distribution is found in the factor :

while the coil limit of the distribution is found in the factor :

with $q=\frac{\nu - 1/3}{1-\nu} \sim 0.617$ [27] and $\nu \sim 0.588$ is the Flory exponent. The final contribution to the whole distribution normalization is the factor $\hat{t}^c\hat{\mathcal{C}}^{(\beta)}_{N}$, where the exponent c ～ -1.13 arises from the cardinality of the space of self-avoiding walks, and ${\mathcal{C}}^{(\beta)}_{N}$ N is a normalization constant depending on the polymer size and the interaction energy US .

Eq.3 for $P_N(\hat{t})$ is relevant to polymers with both excluded-volume and attractive interactions between monomers. This expression predicts that, at a given temperature, small polymers adopt a coil conformation while large polymers adopt a globule conformation, as observed in simulations [28,27,29] and experimentally [30]. In this work, we demonstrate that it also correctly describes the conformation of individual blocks of a copolymer. One notable result of our previous finite-size scaling analysis is that not only the average < $\hat{t}$ > scales differently with respect to the size (i.e. the radius of gyration), but the whole shape of the distribution changes due to finite-size effects [24]. This is shown in Eq.3 as $-N(a_{1}\hat{t} + a_{2}\hat{t}^{2})$ becomes dominant for large N compared to $-a_{3}(N\hat{t})^{-q}$. For a polymer in a coil conformation (high temperature or low attractive interactions), the distribution typically exhibits an exponential peak at $\hat{t} \sim 0$ whereas for a globule conformation (low temperature or high attractive interactions), the distribution exhibits a bell shape around a central value $\hat{t} \gg 0$. At a given temperature T near the so called $\theta$-point (or at intermediate attractive interactions), coil and globule conformations coexist.

We have checked the validity of Eq.3 for homopolymers by fitting it to empirical distributions of $\hat{t}$ obtained from simulations, for various N and various interaction energies. Then we showed that empirical distributions of $\hat{t}$ taken from the conformations of epigenomic domains of simulated copolymers also satisfy Eq.3, demonstrating that each epigenomic domain undergoes its own size-dependent coil-globule transition when its size changes and its interaction energy remains fixed.

2.3. Langevin Dynamics simulation

We used a software library developped in-house, CGMDODE [31], to build a simulator of copolymer dynamics. The simulator builds a chain of N beads from a specified sequence of the form (A_NAB_NB)_NR . We specify the physical parameters of each bead (mass, radius, ...) as well as the specific interaction energy based on its color. For simplicity and to reduce the parameter space, we used the same physical properties for the two colors except for the interaction energy (see parameter table 1). Each bead represents 10kb of chromatin (i.e. ～ 60 nucleosomes). Domains spanned 5 to 55 beads in size, corresponding to actual chromatin sizes of 50 to 550 kb, i.e. the typical range of observed human epigenomic domain sizes. We used the built-in physics engine collision detection algorithm to enforce volume exclusion. A-A (respectively B-B) interactions were implemented by applying a force on each pair of monomers derived according to Eq.1.

We implemented the global Langevin dynamics scheme [32] into our simulator, as recently developped for physics engines [33,34]. The global Langevin dynamics allows for a very fast convergence of conformations towards their equilibrium distribution as well as an efficient sampling of the conformational space at equilibrium, while preserving the accuracy of the ensemble statistics of dynamical properties [33].

Table 1. Langevin dynamics simualtions parameters for beads of 10kb.

description	notation	value
thermostat coupling frequency	$\gamma$	$10^8\text{s}^{-1}$
timestep	$\Delta t$	$5\times10^{-9}$s
total simulation time	$t_{sim}$	$5\times10^{-2}$
Morse potential max. range	$r_{max}$	$7.5 \times 10^{-7}$m
Morse potential width factor	$a$	$10^{7}\text{m}^{-1}$
Morse potential energy minimum distance	$r_0$	$6 \times 10^{-8}$m
bead mass	$m_{b}$	$2.2 \times10^{-20}$kg
bead radius	$r_{b}$	$6 \times 10^{-8}$m
Kuhn length	$b$	$1.5 \times 10^{-7}$m
Temperature	T	300K

 | Show Table

DownLoad: CSV

3. Results

We simulated trajectories of a homopolymer to check Eq.3, then simulated 2-color block copolymers with various domain sizes and interaction energies. We computed the values of $\hat{t}$ from the trajectories for each simulation parameter set. Unless otherwise mentioned, we generated at least 12 independent runs for each simulation parameter set (20 for block copolymers), for a total simulation time of 5 × 10^-2s with a timestep of 5 × 10^-9s (a total of 10⁷ timesteps). We observed that conformations reached equilibrium after 10^-3s of simulation time. We exploited one trajectory frame every 10 000 frames, starting after 2 × 10^-3s of simulation time, to ensure that conformations were at equilibrium. The frame sampling ensures that sampled conformations are independent of each other, as we checked that the auto-correlation function of the snapshot radius of gyration r decreased exponentially with a characteristic time ～ 10^-5s (approximately 5 000 timesteps). We computed $\hat{t}$ for every frame sampled (a single value based on the total polymer conformation for homopolymers, or 3 independent values based on each of the 3 individual block conformations for block copolymers). This yielded, for every condition, at least 8 000 uncorrelated conformations for homopolymers, and at least 40 000 for block copolymer domains. We fitted Eq.3 on the obtained empirical probability densities. Fits were obtained using the Levenberg-Marquardt algorithm for non-linear regression implemented in the package minpack.lm of the R software.

3.1. Size-dependent conformations of a homopolymer

Empirical probability densities of conformations for a homopolymer (N monomers of the same color) are shown in fig.2, for different polymer sizes N and attractive interaction energies U_S (given in absolute value). At low attractive energy (fig.2-A), all polymer sizes yield a peaked distribution at $\hat{t}$ close to 0, characteristic of coil conformations. When attractive interactions increase, the distribution progressively slides towards larger $\hat{t}$ values, and long polymers undergo the transition to globule conformations sooner - i.e. at weaker attractive interactions - than short polymers (fig.2-B and -C). Globule densities stabilize around a central $\hat{t}$ value that depends on the strength of the attractive interaction U_S and the polymer size N. Eq.3 is in very good agreement with simulation data for all N and U_S (fig.2 solid lines).

Figure 2. Size-dependent coil-globule transition for a homopolymer. Empirical probability density for the order parameter $\hat{t}$ (symbols) of the conformation of a homopolymer obtained from simulations, at fixed temperature. N refers to the total size of the polymer. Solid lines are fits for P_N($\hat{t}$) using Eq 3.

DownLoad: Full-Size ImgPowerPoint

3.2. Size-dependent conformations for epigenomic domains

We simulated block copolymers (see fig.1-C for block sequences) with N_A, N_B between 5 and 55. We observed a similar size-dependent coil-globule transition for epigenomic domains of simulated block copolymers, as shown by probability densities of $\hat{t}$ at fixed energies and various domain sizes (fig.3), or at fixed domain size and varying energies (fig.3). In fig.3-A (US = 0.1k_BT), attractive interactions are too low and all domains are in a coil conformation. When attractive strength increases, large domains undergo a transition towards globule conformations, with larger domains transiting first. It appears in fig.3-B at U_S = 0:2k_BT: domains of N=55 monomers form globules (bell-shaped curve on the right of the x-axis) whereas smaller domains keep a coil conformation (peaked curve on the left of the x-axis), and in fig.3-C, domains of sizes between 30 and 55 have transited, domains of size 20 have mixed coil-globule conformations, and domains of size 10 are still mainly coils and just beginning to transit into globules. Eq.3 for P_N($\hat{t}$) was fitted on the densities obtained from block copolymer domain conformations with very good agreement (solid lines fig.3).

Figure 3. Size-dependent coil-globule transition of blocks in a copolymer. Empirical probability densities for the order parameter $\hat{t}$ (symbols) of the conformation of an epigenomic block obtained from simulations, at fixed temperature. N here refers to the size of the block. Solid lines are fits for P_N($\hat{t}$) using Eq.3.

DownLoad: Full-Size ImgPowerPoint

Figure 4. Block copolymer conformation samples. Snapshots of simulated conformations for a block copolymer of sequence (A₂₀B₄₀)₃. Left and middle : attractive interactions U_s = 0.225k_BT. Small domains (green) adopt a coil conformation while large domains (red) adopt a globule conformation. Size-dependent folding produces a vast set of conformation combination possibilities, such as string of globules (left) and central globule with extruding loops (middle). Right : example of microphase separation conformation (two separated globules) for a block copolymer (A₂₀B₄₀)₃, with attractive interactions US = 0.25k_BT.

DownLoad: Full-Size ImgPowerPoint

The main conclusions of this study are, for a given attractive interaction strength US :

· (i) if the temperature is high (or U_S is low), then the scaling of $\hat{t}$ is the same regardless of N : the average radius of gyration of short domains and long domains has the same scaling in N ($\sim N^{\nu}$ for coils, ～ N^1/3 for globules).

· (ii) if the temperature is near $\theta$ ($=\theta(N \to \infty)$) (or attraction strength is close to $1/\theta(N \to \infty)$), then the scaling becomes size-dependent : for a polymer of size N, if T > $\theta(N)$ then the polymer is a coil, whereas if T < $\theta(N)$, the polymer is a globule.

This holds for homopolymers as well as for blocks of a copolymer taken individually, and therefore for epigenomic domains. Our simulations for chromatin showed that a change of interaction energy of 0.1 k_bT is enough to change the scaling for domains of 550 kb (fig.3-A to B, N=55). In terms of size, at U_S = 0.2kT, going from domains of size N = 30 to N = 55 was enough to complete the transition from coil to globule. This suggests that small changes in domain sizes may have huge consequences in TADs and epigenomic domain foldings.

4. Discussion

This work was set in the context of a simple model of chromatin as a block copolymer, stemmed from the observation that finite polymers undergo a size-dependent coil-globule transition. Using Langevin dynamics simulations, we showed that, at fixed temperature, a polymer with attractive and excluded-volume interactions adopts a conformation that depends on its size: large polymers will have a globule conformation and small ones a coil conformation. We demonstrated that this still holds when considering blocks of a 2-color copolymer taken individually: the conformations of the blocks depend on their size. We provided an expression for the probability density of size-dependent conformations and verified its extremely good agreement with simulation results for both homopolymers and block copolymer domains.

Our model uses the same mechanical properties on both simulated epigenomic colors, as only their respective specific interaction energies were modulated. Yet, it is quite possible that mechanical properties of the chromatin fiber might also be affected by epigenetic marks. However experimental data are still lacking in mammalian cells concerning the persistence length of the chromatin fiber - the elusive 30nm fiber- in euchromatin and heterochromatin, note to mention that actual structures of the chromatin fiber are still currently debated [19,35]. Additionnally, the compaction rates (in bp per nm) of chromatin fibers in euchromatin and heterochromatin are most probably very different. When the structural and mechanical properties of various chromatin types will be better characterized experimentally, they could be implemented in a more refined version of our model.

5. Conclusion

The selective folding of domains depending on their size observed demonstrates that two epigenomic domains of the same epigenomic color taken anywhere on the genome will have different conformations if their sequence lengths are different, and notably, that the size difference is not required to be huge for a dramatic change in domain folding. This opens new perspectives for linking epigenomes and chromosome conformation data, as we showed that the span (in terms of sequence) of an epigenomic domain directly influences its folding. This size-dependent folding could link the molecular processes that write / erase epigenetic marks on the chromatin [36] during development and changes in chromosome architecture. During development, by removing the barriers and/or the biochemical cues that delimit domains, epigenetic marks are spread and/or erased by dynamical processes [37,38], leading to modified domain sizes, and thus alternative domain folding. Alternative folding then leads to specific gene expression patterns, as chromosome spatial organization contributes to transcription regulation. Hence, this simple scenario offers new perspectives on how cells can produce a large set of chromosome architectures using a very limited set of epigenomic states. This is why, in future chromatin modelling studies, the effect of the sizes of epigenomic domains should not be overlooked. For instance, combined statistical analysis of conformation and epigenomic data aimed at determining the interaction energy of chromatin colors will fail if the model does not take into account the principle presented in this study. Finally, single-cell HiC experiments [39] should provide invaluable insights regarding the importance of this effect, and its role in gene expression stochasticity and robustness.

Acknowledgments

We gratefully acknowledge the support of NVIDIA which provided computing resources used for this study. This work was funded by the French Agence Nationale de la Recherche, Grant No. ANR- 13-BSV5-0010-03 and by the French Institut National du Cancer, Grant No. INCa 5960.

Conflict of Interest

The authors declare no conflict of interest.