Epigenetic Transitions and Knotted Solitons in Stretched Chromatin

TL;DR

This study explores how mechanical stretching and topology of chromatin influence epigenetic patterning, revealing sharp transitions and knotted solitons that couple 3D structure with epigenetic states, suggesting new experimental avenues.

Contribution

It introduces a novel in silico model linking chromatin stretching, topology, and epigenetic dynamics, highlighting the role of knots and forces in epigenetic pattern formation.

Findings

Chromatin undergoes a sharp transition between disordered and coherent epigenetic states under force.

Knots along chromatin enhance local epigenetic order, creating traveling epigenetic solitons.

Coupling between 3D topology and epigenetics offers new insights into chromatin regulation.

Abstract

The spreading and regulation of epigenetic marks on chromosomes is crucial to establish and maintain cellular identity. Nonetheless, the dynamical mechanism leading to the establishment and maintenance of a given, cell-line specific, epigenetic pattern is still poorly understood. In this work we propose, and investigate in silico, a possible experimental strategy to illuminate the interplay between 3D chromatin structure and epigenetic dynamics. We consider a set-up where a reconstituted chromatin fibre is stretched at its two ends (e.g., by laser tweezers), while epigenetic enzymes (writers) and chromatin-binding proteins (readers) are flooded into the system. We show that, by tuning the stretching force and the binding affinity of the readers for chromatin, the fibre undergoes a sharp transition between a stretched, epigenetically disordered, state and a crumpled, epigenetically coherent, one. We further investigate the case in which a knot is tied along the chromatin fibre, and find that the knotted segment enhances local epigenetic order, giving rise to "epigenetic solitons" which travel and diffuse along chromatin. Our results point to an intriguing coupling between 3D chromatin topology and epigenetic dynamics, which may be investigated via single molecule experiments.

Figures

Equations

$$U_{LJ}^{ab}(x)=\frac{4{ϵ}_{ab}}{\mathcal{N}}\left[{\left(\frac{\sigma }{x}\right)}^{12}-{\left(\frac{\sigma }{x}\right)}^6-{\left(\frac{\sigma }{x_c^{q_a{q}_b}}\right)}^{12}+{\left(\frac{\sigma }{x_c^{q_a{q}_b}}\right)}^6\right].$$

$$x_c^{q_a{q}_b}=\{\begin{array}{l}1.8\sigma \text{ if }{q}_a=q_b\ne 0\\ 2^{1/6}\sigma \text{ if }{q}_a\ne q_b\text{ or }{q}_a=0\text{ or }{q}_b=0.\end{array}$$

$$\Delta U=\sum _b{U}_{LJ}^{ab}(x).$$

$$p=\exp \left(\frac{-\Delta U}{k_B{T}_R}\right),$$

$$m\equiv \frac{\left|n_{\mathrm{red}}-n_{\mathrm{blue}}\right|}{L}$$

$$R_g^2=\frac{1}{L}\sum _{i=0}^L{\left[{\mathbf{r}}_i-{\mathbf{r}}_{\mathrm{CM}}\right]}^2,$$

$$⟨R_g⟩\sim L^{1/3}\sigma ,⟨m⟩>0$$

$$⟨R_g⟩\sim L\sigma ,⟨m⟩\simeq 0.$$

$${\delta }^{2}x(t)=\frac{1}{T-t}{\int }_0^{T-t}{\left[x(t_0+t)-x(t_0)\right]}^{2}d{t}_0$$

$$D=\underset{t\to \infty }{\lim }\frac{{\delta }^{2}x(t)}{2t},$$

$$Ua(x)=\sum _b\left[U_{LJ}^{ab}(x)+U_H^{ab}(x)+U_{KP}^{ab}(x)\right].$$

$$m\frac{d^2{\mathbf{r}}_a}{d{t}^2}=-\gamma \frac{d{\mathbf{r}}_a}{dt}-\nabla U^a(x)+{\mathbf{\xi }}^a(t)$$

$${\tau }_{Br}=\frac{{\sigma }^2}{D_{\mathrm{self}}}=\frac{3\pi \eta {\sigma }^3}{k_B{T}_L}\simeq 0.01\text{s},$$

$$\overline{m}(x)=\frac{1}{\Delta \tau }{\int }_{{\tau }_1}^{{\tau }_1+\Delta \tau }m(l,t)\delta (l-(c_{\mathcal{K}}+x))dt,$$

Full text

I Introduction

Cell-line-specific features in multi-cellular organisms are achieved by regulation of “epigenetic marks”, biochemical modifications of DNA and histone octamers which do not affect the underlying genomic sequence. Thus, the average pattern of epigenetic marks in a given-cell line well correlates with the pattern of genes which are transcriptionally active or inactive within that cell-line Alberts2014 ; Waddington1968 ; Waddington1942 ; Turner2002 ; Henikoff2016 ; Zhang2015a . Dissecting the biophysical mechanisms leading to the de novo establishment, spreading and maintainance of epigenetic marks is consequently a key step towards a better understanding of the dynamic organisation of genomes and of chromosomal re-arrangement throughout the cell cycle Naumova2013 , in cellular ageing Pal2016 and pluripotency Yamanaka2010 . Notwithstanding their pivotal role, these mechanisms are still poorly understood Cortini2016 .

The epigenetic patterning of chromatin – the fibre made of DNA wrapped around histone proteins Alberts2014 – has been shown to strongly correlate with the three-dimensional (3D) nuclear organisation of interphase chromosomes Rao2014 ; Beagrie2017 ; Barbieri2012 ; Brackley2016nucleus ; Jost2014 ; Dixon2012 ; Brackley2016nar . For instance, transcriptionally active regions can be co-localised in multi-enhancer hubs Beagrie2017 or transcription factories Cook2001book , whereas transcriptionally inactive regions may form large heterochromatic foci Gilbert2004 , mega-base (Mb) size lamin-associated domains Guelen2008 or Barr Pinter2012 and Polycomb bodies Wani2016 . This intimate connection is also further supported by mean replication timing data Yaffe2010 ; Baker2012 ; Ryba2010 ; Boulos2015 ; Boulos2014 ; Julienne2013 ; Boulos2013 .

On the other hand, it is important to realise that the establishment of epigenetic patterns is fundamentally a dynamic process, where biochemical tags are constantly deposited, removed and degraded on histones, which can themselves be displaced during transcription or replaced after replication Ramachandran2015 ; Zentner2013 ; Probst2009 ; Saksouk2015 ; Zhang2015a . For this reason, simple models where epigenetic marks are stably deposited along chromatin Sexton2012 ; Filion2010 ; Rao2014 ; Jost2014 ; Barbieri2012 ; Brackley2016nucleus are only crude approximations of a much more complex and dynamic scenario. Crucially, these “static” models fail to address key questions such as, how epigenetic patterns are first established along chromosomes, and how these change, for instance, with cellular ageing Pal2016 or during disruptive events in the cell cycle Zentner2013 . In addition, an understanding of cell-to-cell variability in genome organisation Nagano2013 ; Nagano2016 and of efficiency of the induced-pluripotency pathway Yamanaka2007 may be achieved only through models which can take into account the plasticity of the epigenetic landscape.

In recent years, wide-spread technological advances in the field of molecular biology allowed the biophysical community to identify some of the key players in the dynamics of epigenetic patterning Henikoff2016 ; Talbert2006 ; Zentner2013 ; Pal2016 . At the heart of this process are proteins which can “read” and “write” biochemical tags along chromatin Zentner2013 ; Zhang2015a ; Willyard2017 . Importantly, some of these proteins are found in the same complex Zentner2013 or are known to directly recruit one another Muir2016 . For instance, the heterochromatin binding protein HP1 (a reader) possesses a chromodomain recognising tri-methylation of Lysine 9 on histone 3 (H3K9me3) Lachner2001 ; Hiragami2016 ; Kilic2015 ; at the same time, it can recruit the methyltransferase complex SUV39h1 Zentner2013 (a writer), therefore engaging a positive feedback loop which self-perpetuates this repressive mark Muir2016 ; Narlikar1 ; Narlikar2 .

Computational models of this positive feedback loop in 1D Dodd2007 ; Dodd2011 have shown that accounting for long-range contacts is necessary to allow spreading of repressive marks. On the other hand, 3D models coupling epigenetic and chromatin dynamics Michieletto2016prx have shown that spreading of silencing marks on a mobile 3D fibre can be viewed as an “all-or-none” transition, where a local fluctuation triggers an epigenetic wave of repressive marks which stabilises a compact globular state Michieletto2016prx .

Existing experimental and computational studies suggest that the feedback between 3D chromatin structure and epigenetic dynamics along the chromosome may be a key potential mechanism underlying the establishment of epigenetic patterns which can display memory of their landscape. However, direct experimental observations of dynamical epigenetic and chromatin conformational changes in the nucleus are largely missing, in view of the difficulty to achieve enough spatio-temporal resolution in vivo. In light of this, here we propose, and investigate in silico, a novel avenue to directly observe the coupling between epigenetic dynamics and 3D chromatin folding.

The system we consider can be recreated in vitro via single-molecule experiments: it is inspired both by experiments studying the stretching of reconstituted chromatin Kruithof2009 ; Meng2015 ; Cui2000 , and by investigations of protein-DNA and protein-chromatin interactions in vitro Kanke2016 ; Stigler2016 . Specifically, we envisage a set-up where reconstituted chromatin is attached at both ends to macroscopic beads so as to exert stretching forces on the strand via optical or magnetic tweezers Kruithof2009 ; Cui2000 . Further, we imagine that the chromatin fibre is embedded in a solution where reader and writer proteins are added and activated. In practice, one may achieve this by including “readers” such as heterochromatin HP1 proteins Canzio2013 (which are known to be multivalent, hence act as bridges which fold chromatin), and “writers” such as SUV39h Zentner2013 . These proteins can respectively read and write the repressive H3K9me3 mark Canzio2013 ; Muir2016 . An alternative option is to use Polycomb PRC complexes and Enhancer-of-zeste (E(z)) proteins which are respectively able to bind and deposit the H3K27me3 mark on histones Talbert2006 ; Angel2011 ; Zentner2013 ; Laprell2017 ; Ciabrelli2017 .

Our results show that by tuning the strength of the external stretching force $f$ and the binding affinity of reader proteins $\epsilon$, the system can display either a stretched, epigenetically disordered, state (SD) or a compact, epigenetically ordered, one (CO). The two regimes are separated by an abrupt transition line which we quantitatively locate in the $(f,\epsilon)$ parameter space. We further show that when knots are tied along the chromatin, these can readily be identified by looking at the local order of epigenetic patterning. We argue that the topology of the knots localise and enhance 3D chromatin interactions, thereby “protecting” the epigenetically ordered region, which remains localised within the knotted arc. We dub these novel and remarkable states “epigenetic knotted solitons”.

We finally envisage that our findings will inform the design of novel single-molecule experiments in vitro and illuminate the interplay between epigenetic dynamics and chromatin topology in vivo.

I.1 Chromatin Model and Experimental Set-up

We model the chromatin as a “recolourable” bead-spring polymer (for further details see Methods). Each bead has a size $\sigma$ – for definiteness, we take $\sigma=3$ kb or $30$ nm Thoma1979 ; Widom1985 ; Robinson2006 , but our results do not depend on this precise choice – and it bears a “colour” representing a specific biochemical tag (see Fig. 1). For instance, one may think of blue beads as representing chromatin regions with an excess of heterochromatin mark, H3K9me3, and red beads as ones with excess of Polycomb mark, H3K27me3. We further include grey beads which represent unmarked regions of chromatin. For simplicity, in this work we consider only these three colours, in agreement with the generic experimental observation that few epigenetic “states” are sufficient to well capture the overall epigenetic landscape in several organisms Roudier2011 ; Liu2011 ; Sexton2012 ; Julienne2013 ; Julienne2013a . 111Extending our model to account for more colours (epigenetic states) do not change the qualitative behaviour of the system.

Beads interact through a Weeks-Chandler-Andersen potential, which strongly suppresses any overlap between beads. This interaction is further modified to account for attractive interactions (with affinity $\epsilon$) between beads bearing the same colour. This effectively mimics the action of reading proteins, or “bridges”, which can bring together distant chromatin segments bearing the same epigenetic mark. Unmarked (grey) or differently marked beads are therefore considered to be solely sterically interacting. The truncated Lennard-Jones (LJ) potential can be written as a function of distance $x$ between, and the colour $q$ of, the interacting beads $a$ and $b$ as follows

[TABLE]

for $x\leq x_{c}^{q_{a}q_{b}}$, whereas $U^{ab}_{LJ}(x)=0$ otherwise ($\mathcal{N}$ is a normalisation constant, see Materials and Methods). In Eq. (1) $x_{c}^{q_{a}q_{b}}$ is the cut-off distance between beads $a$ and $b$ with colors $q_{a}$ and $q_{b}$ respectively. This potential is a computationally efficient way to model repulsive and attractive interactions between beads. Specifically, it can display an attractive (negative) region by setting the cut-off larger than the minimum of the potential or can yield pure repulsion otherwise. In practice, we use the knowledge of the colours $q_{a}$ and $q_{b}$ to define a colour-dependent cut-off as

[TABLE]

The shift (last two terms in eq. (1)) ensure that there is no discontinuity in the potential. The binding affinity $\epsilon_{ab}=\epsilon$ when $q_{a}=q_{b}\neq 0$ and is a free parameter of our model, otherwise $\epsilon_{ab}=k_{B}T_{L}$ when $q_{a}\neq q_{b}$ or one of the two is unmarked.

The “recolouring” process is modelled via a Monte Carlo procedure which occurs at an inverse rate of $k_{R}^{-1}=\tau_{R}=10^{3}\tau_{Br}\simeq 10$ s (see Materials and Methods). Every time a bead ($a$) is selected for a recolouring attempt, its colour is randomly changed into one of the two remaining colours and the new energy is computed. Because the only colour-dependent potential employed in the simulations is the one written in Eq. (1), the difference in energy is given by

[TABLE]

The recolouring attempt is then accepted with a Metropolis probability

[TABLE]

where $T_{R}$ is a generic temperature that can, in principle, be tuned according to the efficiency of the writing process and is therefore independent on the solution temperature $T_{L}$, which regulates the stochastic dynamics of the beads in 3D Michieletto2016prx (see Materials and Methods). The recolouring process and the motion of the chromatin in 3D are engaged in a positive feedback loop: beads that are near each other in space are more likely to become equally coloured and beads bearing the same colour are more likely to stick together. Thus, both species (red and blue) compete over the chromatin strand in a similar way that “up” and “down” spins compete over a mobile string, or a “magnetic polymer” Garel1999a .

Here, we study the behaviour of the system upon tuning the stretching force $f$ and the attraction strength $\epsilon$, since this protocol may be realised in vitro by using reconstituted chromatin Kruithof2009 ; Cui2000 and proteins such as HP1 and SUV39 Narlikar1 ; Narlikar2 , as previously mentioned. For simplicity, we limit to the case $T_{L}=T_{R}$ which ensures that the epigenetic read-write mechanism and the chromatin folding are governed by transition rules between different microstates that obey detailed balance and that can be described in terms of an effective free energy. Considering out-of-equilibrium conditions Michieletto2016prx leads to transitions between states with similar epigenetic patterns, i.e. from swollen ordered (or disordered) to compact ordered (or disordered), which can be understood as the homopolymer (or heteropolymer) limit of our system. Because these transitions do not shed light onto the interplay between epigenetics and chromatin conformations, we here decide not to pursue them.

In order to quantify the behaviour of the system we perform Brownian Dynamics simulations for a typical runtime of $5$ $10^{6}$ $\tau_{Br}$ – corresponding to $1000$ $\tau_{R}$, or Monte Carlo sweeps – of a chain $L=1000\simeq 3$ Mb beads long (if not specified otherwise). Although chromatin has not yet been reconstituted to such a length, the trends we uncover are generic and also hold for smaller values of $L$.

II Results

II.1 Epigenetic and Conformational Transitions of Stretched Chromatin

We initialise the systems by assuming that no reader or writer protein is present at $t<0$ and by letting the chromatin equilibrate while subject to the stretching force $f$ and in a bath at temperature $T_{L}$. In practice, we do this by lowering the cut-off for the Weeks-Chandler-Andersen potential written in Eq. (1) to $2^{1/6}\sigma$ for any pair of beads so that the interactions between beads are purely repulsive.

At time $t=0$ the system is flooded with reader and writer proteins which can bridge beads with the same epigenetic mark (i.e., colour) and to attempt recolouring of the beads (i.e., writing) at a rate $\tau^{-1}_{R}$. To visualise the epigenetic evolution of the system, in Fig. 2(A-B) we report typical “kymographs”, which show the colour ($q$) of each bead along the polymer at a given time $t$. Kymographs readily capture the epigenetic “ordering” of the chromatin strand as a function of time for a specific choice of $\epsilon$ and $f$. These two latter quantities are the main free parameters of the system: the stretching force $f$ can be controlled in vitro through either optical or magnetic tweezers, while the binding affinity between readers and epigenetic marks may be varied by considering mutants of bridging proteins such as HP1 Narlikar1 ; Canzio2013 .

In order to further quantify the behaviour of the system in parameter space, we systematically vary $f$ and $\epsilon$ and quantify the equilibrium states by measuring the values of the absolute “epigenetic magnetisation”

[TABLE]

and of the average chromatin extension via its radius of gyration

[TABLE]

where $\bm{r}_{i}$ and $\bm{r}_{\rm CM}$ are the positions of segment $i$ and of the centre of mass of the chain, respectively. In order to construct an equilibrium phase diagram of the system, we estimate the average $\langle R_{g}\rangle\equiv\langle R_{g}^{2}\rangle^{1/2}$ and $\langle m\rangle$, first by time-averaging over a trajectory at steady state and then averaging the results across several trajectories. By measuring these two observables for different choices of $\epsilon$ and $f$ we report the phase diagram shown Fig. 2(C). One can readily notice that the typical equilibrium configurations can be grouped into two distinct phases: one is compact and epigenetically ordered (CO), i.e. with

[TABLE]

whereas the other is stretched and epigenetically disordered (SD), i.e. with

[TABLE]

These two regimes are separated by a transition line and one can readily appreciate from the profiles of $\langle R_{g}\rangle$ and $\langle m\rangle$ (Figs. 2(D-E)) that this transition is abrupt. In the limit of stretching forces $f\rightarrow 0$, we retrieve the first-order transition observed in Ref. Michieletto2016prx for un-stretched chromatin. When larger stretching forces are applied, the system therefore displays a force-dependent critical line $\epsilon_{c}(f)$ (or $f_{c}(\epsilon)$) which retains the first-order-like features observed at $f=0$.

Within the CO region, there is a parameter range where a multi-domain epigenetic structure emerges dynamically. This can be readily seen from the profile of the epigenetic magnetisation $\langle m\rangle$ in Fig. 2(E): near the transition line (from the compact ordered side) this quantity displays a sharp peak which then decreases for smaller values of the force. In the region labelled as “multi” in Fig. 2(E), the magnetisation is lower than unity because multiple macroscopic ordered regions populated by different epigenetic marks compete with one another, thereby lowering the overall magnetisation. These domains must be metastable (though long-lived) as in steady state a single domain is preferred since it removes domain walls, which have a free energy cost. The existence of multi-domain patterns arises because, at stretching forces far from critical $f_{c}(\epsilon)$, there may be multiple nucleation points along the fibre for the spreading of epigenetic marks. These nucleation points are generated by a local increase in chromatin density, for instance through the transient formation of loops and 3D interactions, which then trigger local spreading of different epigenetic marks (see Fig. 3). As the stretching force gets closer to the critical value $f_{c}(\epsilon)$, the nucleation probability declines, hence only one epigenetic mark is able to take over the whole chromatin strand, thereby enhancing its overall epigenetic magnetisation.

It is useful to compare the epigenetic transition between the CO and the SD state in Fig. 2 with the equilibrium transition between a compact and a stretched state which can be observed in a stretched homopolymer with self-attractive interactions (i.e., in a poor solvent). While the transition is first order in the homopolymer case as well homopolymer ; homopolymer2 , here the epigenetic degrees of freedom increase the entropy of the disordered phase, thereby the value of the critical force, for a given $\epsilon$, is lower. Another key difference is that multi-pearl structures are only observed transiently close to the transition for the homopolymer case homopolymer , whereas in the epigenetic case multi-domain states arise far from the transition, and they are long-lived. Both these differences should be experimentally detectable, as the homopolymer case may be recreated by reconstituting chromatin fibres with histone octamers with controlled biochemical tags or in the presence of linker histones (H1), while avoiding the action of writers.

II.2 Long-Lived Epigenetic Domains and Boundary Diffusion

The competition of epigenetic marks over a chromatin segment is a phenomenon relevant to many biological systems. For example, genes that are positioned near telomeres may switch between transcriptionally active and silenced states from one generation to the next Gottschling1990 . This phenomenon, generally referred to as “telomeric position effect”, is particularly relevant in yeast Gottschling1990 but it has also been observed in human cells Baur2001 . It occurs because heterochromatin marks largely populating the telomeric gene-poor regions of the genome invade gene-rich stretches of chromatin, effectively switching off gene transcription.

In our proposed single molecule set-up, the emergence of multiple competing epigenetic domains in reconstituted chromatin under (small) tension allows a quantitative study of stability and diffusive dynamics of boundaries between neighbouring epigenetic domains, processes which may be relevant to understand the mechanisms underlying competition and spreading of epigenetic marks in vivo.

In practice, an epigenetic domain in the stretched chromatin fibre is defined as a macroscopic region (e.g., consisting of $w\geq 50$ beads) over which more than {\color[rgb]{0,0,0}\theta}=90\% of the beads are homogeneously coloured. A boundary between domains is found if the beads on either side of the boundary have opposite colours (i.e., red and blue) and the signed magnetisation difference $\Delta\tilde{m}$ of the left and right domain is such that \left|\Delta\tilde{m}\right|>2{\color[rgb]{0,0,0}\theta} (i.e. the domains have opposite signed magnetisation).

By tracking the position of the domain boundaries over time, one can measure their mean square displacement

[TABLE]

where $x(t)$ is the position of a given boundary at time $t$ and $T$ the total measurement time. Plots of $\delta^{2}x(t)$ for different replicas and two choices of interaction energies $\epsilon$ are reported in Fig. 3(C-D). Data points are grouped according to the applied stretching force $f$, showing that the dependence on $f$ is much weaker compared to that on $\epsilon$. We extract the diffusion coefficient $D$ as the long time behaviour of $\delta^{2}x(t)$, i.e.

[TABLE]

which is reported in Fig. 3(F) and shows an exponential decay of $D$ as a function of $\epsilon$.

In our stretched chromatin assay, the diffusive dynamics of epigenetic boundaries is therefore mainly controlled by the strength of the attractive interactions between beads bearing the same epigenetic mark. In vivo, another possible important factors affecting diffusivity of epigenetic boundaries may be the presence of insulators and architectural proteins such as CTCF and cohesins Alberts2014 ; Stigler2016 ; Kanke2016 . These features may be included in future studies focused on understanding the actions of these architectural elements.

II.3 Epigenetic Knotted Solitons

The abrupt, first order, phase transition between the compact-ordered and the stretched-disordered state shown in Figure 2 is here due to the coupling between global 3D structure and 1D epigenetic dynamics Michieletto2016prx , which gives rise to a positive feedback loop where compaction aids spreading, which leads to further crumpling at small enough $f$, or sufficiently large $\epsilon$.

It is interesting to ask whether one may design an experimentally realisable situation whereby this positive feedback is only realised locally, rather than globally as done in Ref. Michieletto2016prx . To do so, we consider a knotted chromatin fibre: the idea behind this construct is that a polymeric knot tightly localises upon stretching Farago2002 ; Caraglio2015 , and that the knotted region is more compact, thereby providing a natural nucleation point for epigenetic ordering. In the single molecule experiment set-up we consider (Fig. 2), a physical knot may be tied within chromatin by using, for instance, micro-manipulating techniques used to generate knots on DNA arai1999tying ; Bao2003 ; Liu2016 . Alternatively, a knotted chromatin fibre may be self-assembled by first knotting the naked DNA strand, and later allowing nucleosome formation by adding histone octamers to the solution Brackley2015nar .

In Figure 4 (A-B) we show the kymograph for a chromatin strand along which a figure-of-eight ($4_{1}$) knot is tied. As one can see, even though we choose values of $\epsilon$ and $f$ so that the system is in the stretched-disordered region of the phase diagram (see Fig. 2), the kymograph clearly displays a localised epigenetically ordered region. To confirm whether this region corresponds to the knotted arc, tightened upon stretching and self-attraction, we pinpoint and monitor the time evolution of the knotted region by resorting to a top-down search of the smallest portion of chromatin that yields, upon suitable closure, a ring with the same topology of the whole chain Tubiana2011 . Remarkably, the knotted portion found with this well-established algorithm perfectly matches the epigenetic domain in the kymograph (black lines in Fig. 4 identify the boundary of the knotted arc, see also Methods).

We call this remarkable configuration an “epigenetic knotted soliton”, because, as a soliton, it travels along the fibre by keeping near-constant shape, hence displaying particle-like behaviour. Similarly to solitons recently observed in meta-materials Chen2014 and minimal surfaces Machon2016 , it is the non-trivial topology of the localised knot which keeps the structure together. The epigenetic state of the knotted soliton can “switch” from red to blue, or vice versa – switching events occur, in Figure 4, at about $2$ $10^{5}\tau_{Br}$ and $3.5$ $10^{5}\tau_{Br}$). Switching between the two stable states is observed because the soliton has a finite size, hence the effective free energy barrier separating the red and blue states can be occasionally overtaken by fluctuations.

The physical properties of the soliton are tunable: for instance, by increasing the stretching force the knot tightens and becomes more compact, so that the soliton shrinks in size but becomes more stable due to the increase of 3D interactions which locally “protect” the epigenetically ordered state. By increasing values of $\epsilon$, such that the system remains overall disordered ($\epsilon<\epsilon_{c}(f)$) the switching rate decreases (i.e., the epigenetic ordering is more robust, see Fig. 4(D)) – the diffusion coefficient, however, does not depend on $\epsilon$ appreciably (unlike for the case of epigenetic domains, see Fig. 3(E-F)). It would be of interest to quantify how the soliton diffusion coefficient depends on knot type, and whether, as in knots in swollen polymers, twist knots and more complex knots are less mobile Adams1994 ; Bao2003 ; diStefano2014 ; matthews2010 .

For values of $\epsilon>\epsilon_{c}(f)$, we observe the same SD-CO phase transition reported in Figure 2, where the localised knot now acts as a nucleation point for the spreading of its epigenetic mark.

Finally, we investigate how knotted solitons interact on a chromatin strand. We consider a system where two knots, a pentafoil $5_{1}$ and a trefoil $3_{1}$, are initially tied along the chromatin. As the system evolves each knot forms its independent epigenetic soliton. As shown by the kymographs in Figure 5(A-D), both solitons diffuse along the chain Bao2003 ; they also interact and may merge to form a single domain, and later on split again (Fig. 5(A)), or even cross each other Trefz2014 (Fig. 5(D)). Following every collision or interaction, each soliton retains its overall structure and re-establishes its locally coherent epigenetic mark.

III Discussion and conclusions

In several areas of chromosome biology, researchers assume the existence of a tight coupling between 3D chromatin structure and epigenetic dynamics, as this coupling provides an appealing mechanism for the de novo establishement and maintainance of epigenetic patterns Dodd2007 ; Michieletto2016prx ; Canzio2013 ; Angel2011 ; Ciabrelli2017 . However, as both epigenetic spreading and chromatin dynamics may occur on comparable and relatively fast timescales in vivo (minutes to hours Zentner2013 ; Barth2010 ), it is extremely difficult to design an experiment to demonstrate this coupling dynamically inside the cell. Here, instead, we proposed, and investigated in silico, a set-up for an experiment to test and quantify the coupling directly in vitro. The set-up involves a reconstituted chromatin fibre Kruithof2009 ; Cui2000 , which is stretched (e.g., by laser tweezers) and interacts with an ensemble of reader and writer proteins – such as heterochromatin HP1 and SUV39h1 complexes Narlikar1 ; Canzio2013 or polycomb Angel2011 ; Ciabrelli2017 enzymes able to bind to and deposit their respective marks.

We showed that by varying, for instance, the stretching force, one can trigger a phase transition between a compact-ordered phase, where a single epigenetic mark invades the whole chromatin fibre, and a stretched-disordered phase, where no single mark can take over the system. It is important that this transition may be observed by varying only the stretching force, as this is simpler to fine tune continuously in vitro with respect to temperature (which can inactivate readers or writers altogether) and effective self-attraction between chromatin segments (which may be achieved through the use of mutant reader proteins).

By using single-molecule imaging Yardimci2012 ; Kanke2016 ; Stigler2016 or super-resolution techniques Huang2008 one may envisage to employ antibodies recognising H3K9me3 and H3K27me3 epigenetic states to separately stain the marks and thereby discriminate between an epigenetically ordered and an epigenetically disordered state. In principle, one may even be able to observe the transition from one to the other by increasing the stretching force in the single molecule assay. Another possibility to assess the degree of epigenetic order is to use single-cell chromatin immunoprecipitation Rotem2015 (ChIP) on the reconstituted chromatin, to quantify histone modifications along the fibres. Note that, for this approach to be viable, one would have to use genomic DNA to reconstitute chromatin (rather than repeating, such as 601, sequences), so that locations along the DNA can be mapped uniquely. Therefore, whilst an experimental investigation of the compact-ordered to stretched-disordered transition is certainly challenging, and will require state-of-the-art experimental techniques, it is in principle feasible, and would constitute the first direct measurement of the so far elusive coupling between chromatin structure and epigenetic dynamics. A more experimentally accessible insight (although less informative regarding the epigenetic state of the system) may also be achieved via optical-tweezers by measuring force-extension curves Kruithof2009 ; Cui2000 ; VanDenBroek2010 . These can discriminate between a (partially) collapsed and a stretched coil VanDenBroek2010 , while cannot determine whether the chromatin is epigenetically ordered or disordered. On the other hand, the typical force-extension curves obtained with this set-up may display sensitivity on the concentration or efficiency of epigenetic readers and writers thereby offering indirect quantification of the recolouring process.

By tying a knot along the reconstituted chromatin fibre we observed that the system can harbour a novel structure, an “epigenetically knotted soliton” (see Figs. 4- 6). This is a tight knot which is locally ordered epigenetically and diffuses freely, within an epigenetically disordered background. The size of the soliton is tunable by varying the stretching force applied to the chromatin fibre, and different solitons interact in a variety of ways: they may bounce off one another or cross through each other when they collide. We also expect similar topological solitons to be universally found in knotted magnetic polymer, a new kind of topological soft matter which has not yet been realised experimentally Garel1999a .

The findings we have reported here may also be employed to detect knots in chromatin strands. Starting from an ordered and crumpled state (such as heterochromatin Alberts2014 ), one may imagine to apply and continuously increase an external stretching force $f$ and monitor the evolution of the system; whereas ordered unknotted regions will undergo an abrupt phase transition and become epigenetically disordered (see phase diagram in Fig. 2), knotted chromatin will instead preserve a localised epigenetically ordered region that can be identified with the knotted segment (Fig. 4).

Besides all this, the epigenetic solitons provide another observable consequence of the dynamic coupling between 3D chromatin structure and epigenetics which can be tested in future single-molecule experiments.

IV Methods

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

IV.1 Chromatin Model

Chromatin is modelled as a bead-spring polymer chain Kremer1990 where each bead has nominal size $\sigma$. Attractive and repulsive interactions are controlled by the potential described in the main text (Eq. (1)) and plotted in Fig. S1. This is a truncated-and-shifted Lennard-Jones potential which is broadly used in Molecular Dynamics FrenkelMD to model short-ranged interactions. We include a colour-dependence in order to model like-colour attraction and different-colour repulsion as described in eq. (1) and shown in Fig. S1. Harmonic springs between beads $a$ and $b$ are imposed as $U^{ab}_{H}(x)=(\delta_{b,a+1}+\delta_{b,a-1})k_{B}T_{L}\kappa(x-x_{0})/2$ with $x_{0}=1.1\sigma$ and $\kappa=100k_{B}T_{L}$ to ensure the connectivity of the backbone. The chain stiffness is regulated by a Kratky-Porod potential between triplets of beads forming an angle $\theta=\left(\bm{t}_{a}\cdot\bm{t}_{b}\right)/\left(\left|\bm{t}_{a}\right|\left|\bm{t}_{b}\right|\right)$ where $\bm{t}_{a}$ is the vector joining beads $a$ and $a+1$ as $U^{ab}_{KP}(\theta)=(\delta_{b,a+1}+\delta_{b,a-1})k_{B}T_{L}l_{p}\left(1+\cos{\theta}\right)/\sigma$ with $l_{p}=3\sigma$. The total potential $U^{a}(x)$ experienced by each bead is given by the sum over all the possible interacting pairs and triplets, i.e.

[TABLE]

The dynamics of each bead therefore obeys the Langevin equation

[TABLE]

where $\gamma=1$ (in dimensionless LJ units) is a friction coefficient and $\bm{\xi}^{a}(t)$ is stochastic delta-correlated noise which obeys the fluctuation dissipation relationship $\langle\xi_{\alpha,a}(t)\xi_{\beta,b}(t^{\prime})\rangle=2\gamma k_{B}T_{L}\delta(t-t^{\prime})\delta_{\alpha\beta}\delta_{ab}$ where the Latin indexes run over particles and Greek indexes over Cartesian components. The mass of the beads are taken to be unity in the dimensionless Lennard-Jones units FrenkelMD . The beads employed at the ends of the chromatin are five times larger than the beads forming the polymer and have the same mass. In order to simulate the stretching of the chromatin, we directly apply a force $f$ on these beads. Finally, Eq. (11) is integrated with a velocity-Verlet scheme within the LAMMPS Plimpton1995 engine.

Simulations units can be mapped to real ones by considering that the polymer beads can coarse-grain groups of nucleosomes. The nominal size of the beads can is here considered to be 30 nm, as capturing the thickness of reconstituted chromatin fibre Thoma1979 ; Widom1985 ; Robinson2006 ; Kruithof2009 . The typical timescales over which each bead diffuses its own size is therefore

[TABLE]

where we used the Einstein relation $D_{\rm self}=k_{B}T_{L}/3\pi\eta\sigma$ and we considered $\eta=150cP$ as the effective viscosity.

The recolouring process occurs at inverse rate $\tau_{R}=10^{3}\tau_{Br}$. In each attempt, a bead is selected at random and a colour change is proposed; the move is then accepted if satisfying the Metropolis criterion as described in the text. In other words, on average, all beads will have been attempted to change their colour every $\tau_{R}=10^{3}\tau_{Br}$ steps.

For simplicity, we model reader and writer proteins implicitly. This is effectively done by setting attractive pair-wise interactions between beads bearing the same colour and by assuming a uniform selection probability during the recolouring attempts. In other words, the physical proximity of reader and writers is not required for bridging segments and changing the beads colour. This model is therefore a good approximation for the case in which reader and writer proteins saturate the solution. On the other hand, lowering the concentration of reader and writer proteins would effectively slow down the kinetics towards the equilibrium state: this situation is more complicated and we defer it to future studies.

IV.2 Knotted Soliton

To test the stability of the knotted soliton we consider the time-averaged degree of epigenetic coherence $\overline{m}(x)$ at position $x$ relative to the centre of knotted arc, $c_{\mathcal{K}}$. This can be defined as

[TABLE]

where $\Delta\tau$ is the typical time between epigenetic “switching”, where the magnetisation of the whole knot changes sign (such as the one occurring in Fig. 4 at about $2$ $10^{5}\tau_{Br}$ and $3.5$ $10^{5}\tau_{Br}$) and $m(l,t)$ is the magnetisation at position $l$ and time $t$. In Fig. 4(C) we report $\overline{m}(x)$ for three system replicas and compare it with the distribution of knot sizes obtained from the 3D polymer configurations (inset). As one can notice, the length of the ordered domain agrees with the length of the knotted arc. Further, the shape of the magnetisation profile $\overline{m}(x)$ displays a broad peak, which is characteristic of soliton-type solutions.

IV.3 Acknowledgements

We acknowledge ERC for funding (Consolidator Grant THREEDCELLPHYSICS, Ref. 648050). The authors would like to thank Nick Gilbert and Jim Allan for insightful discussions and useful remarks on the manuscript.

IV.4 Author contributions statement

All authors conceived the research. D.Mi. performed the simulations and analysed the data. All authors wrote and reviewed the manuscript.

IV.5 Additional Information

The authors declare no conflict of interest.