---
abstract: '**Understanding the spatial organisation of the genome in the cell nucleus is one of the current grand challenges in biophysics. Certain biochemical – or epigenetic – marks that are deposited along the genome are thought to play an important, yet poorly understood, role in determining genome organisation and cell identity. The physical principles underlying the interplay between epigenetic dynamics and genome folding remain elusive. Here we propose and study a theory that assumes a coupling between epigenetic mark and genome densities, at the scale of the whole nucleus. We show that equilibrium models are not compatible with experiments and a qualitative agreement is recovered by accounting for non-equilibrium processes which can stabilise . We finally discuss the potential biophysical origin of these terms.**'
author:
- Davide Michieletto
- Davide Colì
- Davide Marenduzzo
- Enzo Orlandini
bibliography:
- 'library.bib'
- 'nota1.bib'
- 'nota2.bib'
- 'notaglotzer.bib'
title: |
    Nonequilibrium Theory of Epigenomic\
    Microphase Separation in the Cell Nucleus
---

[^1]

[^2]

Establishing distinct and inheritable cellular identities in different tissues is crucial to the existence of complex multi-cellular organisms. Because all cells contain the same DNA, cellular fate cannot be directed by genetic cues alone. Instead, tissue-specific 3D genome organisation [@Cavalli2013] and biochemical (also called epigenetic) patterns [@Cortini2015] are thought to be key regulators of cellular fate [@Stadhouders2019]. Epigenetic patterns are composed of biochemical marks on DNA and histones – the proteins which package DNA into chromatin [@Alberts2014] – and can be deposited or removed by a range of specialised proteins [@Cortini2015] often recruited by complexes such as RNA polymerase [@Ng2003] and Polycomb Repressive Complex (PRC) [@Zentner2013].

Investigating the interplay between genome organisation and epigenetic patterns can shed light into the mechanism underlying cell fate decision [@Pombo2015; @Stadhouders2018]. To study the dynamics of epigenetic patterns, Ising-like models have been proposed in the biophysics literature [@Dodd2007; @Jost2014pre; @Teif2015; @Berry2017; @Erdel2016; @Erdel2018a]. Yet, the genome is poorly represented by a 1D array of Ising spins. Instead, it is a fluctuating polymer which can assume distinct spatial organisations in 3D and also in response to stimuli [@Zirkel2017; @Buckle2018] and epigenetic cues [@Stadhouders2019].

In this Letter, and the companion paper [@PRE], we aim to analyse the link between chromatin and epigenetic dynamics by using models for chromatin folding inspired by the physics of magnetic polymers [@Garel1999a]. We study a Landau-Ginzburg field theory where the dynamics of epigenetic marks is linked to that of genome folding within the nucleus. Our theory considers a 1D chain of Potts-like spins which is allowed to fluctuate in 3D. It markedly departs from previous works [@Michieletto2016prx; @Michieletto2018nar] as we formalise and study a theory for the epigenomic organisation of a whole eukaryotic nucleus and investigate how non-equilibrium processes affect its thermodynamics and kinetics.

We show that while an equilibrium theory for epigenomic organisation captures some features seen [*in vivo*]{}, such as segregation of different epigenetic marks, it fails to explain the experimentally observed coexistence of diverse epigenetic and genomic domains in eukaryotic nuclei [@Cremer2015; @Cortini2015]. In light of this we then propose a non-equilibrium field theory accounting for generic energy-consuming biochemical and biophysical processes in the nucleus. We discover that a simple first-order reaction leads to arrested (and tunable) phase separation of epigenomic domains, in qualitative agreement with experiments [@Cremer2015]. Alongside our theory, we show results from large-scale 3D Brownian Dynamics (BD) simulations of (non-equilibrium) magnetic polymer melts. Besides validating our field theory, they also represent a key step towards a more realistic non-equilibrium polymer model for genome organisation with dynamic epigenetic marks.

#### Model – 

Epigenetic marks are biochemical (e.g., methyl) groups, that are transiently deposited on histones or DNA. Albeit a slew of modifications exists [@Alberts2014], it is typical to consider a generic situation in which two classes of epigenetic states may be present on a given chromatin segment, marking either transcriptionally active or silent DNA [@Cortini2015].

Specific protein complexes can bind and bridge genomic segments bearing the same marks [@Barbieri2012; @Brackley2013pnas; @Jost2014B] while others are known to deposit biochemical groups [@Hathaway2012]. Importantly, evidence suggest that complexes depositing a given epigenetic mark also contain protein domains which can bind the same mark [@Dodd2007; @Cortini2015] . This read-write mechanism sets up a positive feedback loop which is here captured by a “ferromagnetic” interaction tending to align 3D proximal Potts spins – here playing the role of epigenetic marks – or to bring them together when already aligned [@PRE].

We describe the system with two fields $n({\boldsymbol{\mathrm{x}}},t)$ and $m({\boldsymbol{\mathrm{x}}}, t)$ – the density and “magnetisation” fields. The former records the local genomic density at a given spatio-temporal location, the latter epigenetic state. These fields should be interpreted as local averages over a volume element centred in ${\boldsymbol{\mathrm{x}}}$, large enough to smooth out microscopic fluctuations. Since the total DNA mass is constant during most of the cell cycle, we assume that the mean density $n_0 = \int d{\boldsymbol{\mathrm{x}}} \, n({\boldsymbol{\mathrm{x}}},t)/V$ is conserved. A generic Landau-Ginzburg free energy density describing this model is $$\beta f_{nuc} =  a m^2 + b m^4 + c n^2 + d n^3 -\chi m^2 n \, .
\label{eq:FMelt}$$ The rationale behind Eq.  is that it must respect the $\mathbb{Z}_2$ symmetry of the magnetisation field (active/inactive chromatin [@Lieberman-Aiden2009]) and it must describe the behaviour of the density field via a standard virial expansion for non-ideal gases [@ChaikinLubensky] so that, for $m=0$, the ground state is a nucleus homogeneously filled with DNA. Hereafter we set $a=b=c=d=1$ unless otherwise stated, thereby suppressing phase transitions driven by any of the two uncoupled fields. The key coupling term, $\chi m^2 n$, is chosen to model the feedback between chromatin folding ($n>0$) and epigenetic ordering ($m^2 > 0$). In this respect, $\chi$ can be thought of as parametrising the self-attraction of equal epigenetic marks, which is mediated by bridging proteins [@Barbieri2012; @Brackley2013pnas] such as HP1 [@Larson2017; @Hiragami-Hamada2016], polymerases [@Cook2002] and transcription factors [@Sabari2018].

![**Phase Diagram at the Nuclear Scale.** **A** Sketch of the epigenetic and folding dynamics of magnetic polymers. **B** Equilibrium phase diagram of the free energy density (Eq. ). The three equilibrium phases are: (UD) uniform ($n=n_0$) and epigenetically disordered ($m^2=0$); (UO) uniform ($n=n_0$) and epigenetically ordered ($m^2>0$) and (DO) demixed and epigenetically ordered ($n=n_+$, $m^2>0$ and $n=n_-$,$m^2=0$). The dotted line marks the critical value of the coupling $\chi_c(n_0)=a/n_0$, the solid lines identify the boundaries of the coexistence region (binodals) and the dashed lines identify the spinodal region. Insets report representative snapshots from BD simulations performed at monomer density $\rho$ and self-affinity $\epsilon$ (see SI and Ref. [@PRE] for details). The value of the parameters used for the snapshots are: UD = ($\rho \sigma^3=0.5$, $\epsilon/k_BT=0.7$); UO = ($\rho \sigma^3=0.8$, $\epsilon/k_BT=1$); DO = ($\rho \sigma^3=0.1$, $\epsilon/k_BT=0.9$). See also Supp. Movies **M1-M2** and **BD M1-M3**. ](panel2.png){width="44.00000%"}

\[fig:pd\_melt\]

#### Equilibrium Thermodynamics – 

We seek to thermodynamically characterise the theory in Eq.  first by finding an $m$ minimising the free energy for a given $n$, and then by using the common tangent construction to assess whether the resulting uniform state is unstable to phase separation [@Matsuyama2002; @Fosado2017]. The optimal $m^*$ is $m^* = \pm\sqrt{ (\chi n({\boldsymbol{\mathrm{x}}},t) - a)/ 2b}$ if $n({\boldsymbol{\mathrm{x}}},t) > a/\chi$ and 0 otherwise. This solution indicates a second-order phase transition between an epigenetically disordered ($m^2=0$) and ordered ($m^2>0$) phase at the critical line which is shown as a dotted line in Fig. \[fig:pd\_melt\].

Plugging $m^*$ into Eq.  we obtain the minimised $f^*=f(m^*,n)$, which we need to further minimise with respect to $n$, subject to the constraint that $ \int_V n({\boldsymbol{\mathrm{x}}}) \,\mathrm{d} {\boldsymbol{\mathrm{x}}}/V = n_0$. This can be done via the common tangent construction [@Fosado2017; @Matsuyama2002; @ChaikinLubensky], by finding points in phase space where pressure $P =f^\star - n \partial f^\star/\partial n$ and chemical potential $\mu=\partial f^\star/\partial n$ of the two phases are equal, while having a lower free energy than the mixed/uniform phase. Graphically, these conditions can be solved by finding the points $(n_-, f^*(n_-))$ and $(n_+, f^*(n_+))$ at which the tangents to $f^\star$ have the same slope $\mu$ and the same intercept $P$ [@nota1]. By repeating this procedure for different $\chi$ one finds the so-called “binodals”, $n_-(\chi)$ and $n_+(\chi)$, which are plotted as thick lines in Fig. \[fig:pd\_melt\]. Binodal lines delimit the region of phase separation but contain regions in which the mixed state is metastable. The “spinodal” region of linear instability, where the system spontaneously demixes into low ($n_-$) and high density ($n_+$) phases, is that where $\partial^2f^\star/\partial n^2 < 0$ yielding $a/\chi < n_0 < (\chi^2-4 b c)/(12 b d)$ shown in Fig. \[fig:pd\_melt\] as dashed lines.

This construction allows us to discover three possible equilibrium phases for our theory: (i) for $\chi < \chi_c(n_0) = a/n_0$ the system is in a uniform ($n=n_0$) and disordered ($m=0$) phase (UD); (ii) at $\chi=a/n$ we find a second order phase transition to a uniform state ($n({\boldsymbol{\mathrm{x}}},t)=n_0$) with ordered epigenetic field $(m^2>0)$ (UO); (iii) for $\chi > \chi_c (n_+)$ or $\chi > \chi_c (n_-)$, we observe a phase separated state that we call demixed ordered (DO).

In a biological context, our theory suggests that certain perturbations can trigger nuclear macroscopic phase separation driven by epigenetics. This is achieved either by moving the system into the spinodal region or by overcoming the metastability of the UO phase outside this region.

![**Dynamics and Reorganisation of Epigenomic Domains.** **A** Snapshots following the quench $Q_1$ from Uniform Disordered (UD) to Uniform Ordered (UO) ($\chi=0 \rightarrow \chi=1$ at $n_0=2$). **B** Snapshots following the quench $Q_2$ from Uniform Disordered (UD) to Demixed Ordered (DO) ($\chi=0 \rightarrow \chi=3.5$ at $n_0=0.5$). Top rows in **A-B** show the evolution of Eq.  on a $100\times 100$ grid. Bottom rows show the evolution through BD simulations (UD$\to$UO, $\rho \sigma^3=0.8$, $\epsilon/k_BT=1$; UD$\to$DO, $\rho \sigma^3=0.1$, $\epsilon/k_BT=0.9$). Ref. [@PRE] gives more details on the BD simulations. **C** Scaling of the mean epigenetic domain (units of lattice size $dl$) computed as the average size of clusters of connected lattice sites with magnetisation $|m|>1$ for $Q_{1,2}$, and $|m|>2$ for $Q_3$ ($\chi=0 \rightarrow \chi=6$ at $n_0=1$). See also Supp. Movies **M1-M2** and **BD M1-M3**. The other parameters are set to $\kappa_m=\kappa_n=5$ and $\Gamma_m=\Gamma_m=0.1$. ](panel2dynamics.png){width="42.00000%"}

\[fig:dynamics\]

#### Equilibrium Kinetics – 

. It is therefore important to characterise such kinetics within our theory. To this end we now consider the time-Dependent Landau-Ginzburg equations derived from Eq. . These are obtained considering the steepest descent to the free energy minimum for the magnetisation field, $\partial_t m \sim -\delta \mathcal{H}/\delta m$, and a diffusive dynamics of the conserved density, $\partial_t n \sim \nabla^2 \delta \mathcal{H}/\delta n$, with $\mathcal{H} = \int d{\boldsymbol{\mathrm{x}}} \left[ f + \kappa_n (\nabla n)^2 + \kappa_m (\nabla m)^2 \right]$ [@ChaikinLubensky]. These “Model C” [@ChaikinLubensky] equations can be explicitly written as $$\begin{aligned}
&  \dot{m} = \Gamma_m \left(2 \chi m n  -2 a m - 4 b m^3 + \kappa_m \nabla^2 m \right)  \notag \\
&  \dot{n} = \Gamma_n \nabla^2 \left( 2 c n + 3 d n^2  - \chi m^2 - \kappa_n \nabla^2 n\right) \, . 
\label{eq:modelC}\end{aligned}$$ We numerically evolve these coupled equations on a $100 \times 100$ grid initialising the system in a UD state with some small fluctuations in both magnetisation and density. From here, a quench $Q$:$\chi \rightarrow \chi^\prime$ is performed and the evolution towards a new stable state monitored.

![image](panel4-new.png){width="85.00000%"}

\[fig:Switching\]

Interestingly, the quench into the UO state ($Q_1$ in Fig. \[fig:dynamics\]) follows a different scaling with respect to all the others. In $Q_1$ the density field remains uniform while the epigenetic domains form bicontinuous spanning clusters evolving into one percolating domain (see [**Movie M1**]{}). The scaling of the typical domain size grows as $L(t) \sim t^{\alpha}$ where $\alpha = 0.46$ is compatible with the coarsening of a non-conserved field such as the magnetisation in the Ising model [@Kockelkoren2002]. In all other quenches there is reorganisation of both epigenetic and density fields and our theory yields a slower growth with $\alpha \simeq 0.25 - 0.3$. The lower end is interpreted as a transient behaviour which would eventually lead to the $1/3$ exponent expected in density-conserving Model B kinetics [@Kockelkoren2002] (see **Movie M2**). It would be interesting to follow the growth of epigenomic domains in live cells, e.g. in oncogene-triggered senescence [@Boumendil2019; @Chiang2018]

. Their coarsening is arrested and does not proceed indefinitely as predicted by Eqs. ; for this reason, we now aim to improve our model and account for non-equilibrium processes.

#### Non-equilibrium Model – 

. Some of them can affect the local chromatin state and its accessibility to proteins [@Brangwynne2015] and under certain conditions, they may destabilise or hinder the deposition of histones with a given epigenetic mark [@Henikoff2008].

The dynamical equations with epigenetic switching are $$\begin{aligned}
&  \dot{m} = \Gamma_m \left(2 \chi m n_a -2 a m - 4 b m^3 + \kappa_m \nabla^2 m \right)  \notag \\
&  \dot{n}_a = \Gamma_n \nabla^2 \left(g(n_a) - \chi m^2 - \kappa_n \nabla^2 n_a\right)  + \sigma_{a} n_i - \sigma_{i} n_a  \notag \\
& \dot{n}_i = \Gamma_n \nabla^2 \left(g(n_i) -\kappa_n \nabla^2 n_i\right) -\sigma_{a} n_i + \sigma_{i} n_a
\label{eq:modelC+switch}\end{aligned}$$ where $g(x)=2 c x + 3 d x^2$ and the parameters $\sigma_{a/i}$ describe the rates at which chromatin is activated/inactivated (in general $\sigma_i \neq \sigma_a$). Note that we now require the [*sum*]{} of the two density fields to be conserved [@notaglotzer]. We numerically evolve Eqs.  on a $100 \times 100$ grid starting from the UD phase (other parameters are given in Fig. \[fig:Switching\]).

The rates are set so that any region of the system with $m\ne 0$ can become ; yet, because the field $n_i$ decouples from $m$, occurs through a neutral ($m=0$) state. This sets up a current in epigenetic states which breaks detailed balance and time-reversal symmetry (Figs. \[fig:Switching\]A,B), similarly to what occurs in theories for scalar active matter [@ElsenPRX]. The current is visible in our BD simulations, when following the evolution of a chromatin bead with a given epigenetic mark (Fig. \[fig:Switching\]B, inset).

The non-equilibrium switching terms in Eqs.  yield an arrested phase separation with multiple epigenomic domains in steady state, or microphase separation, as in the majority of cell types [@Cremer2015]. We note there is comparatively much less density variations with respect to equilibrium (Figs. \[fig:Switching\]B,C). For large $\sigma_{a/i}$, chromatin droplets of the same magnetisation form “super-beads” which come together forming branched aggregates (Fig. \[fig:Switching\]D), as this configuration minimises interfaces in the fields $m$ and $n_{i}$ (see Movies **M6**-**M9**). The domain size in steady state is controlled by the switching rate, and we find it approximately decays as $\sigma^{-1/4}$ (Fig. \[fig:Switching\]E).

Large-scale BD simulations of magnetic polymer melts in which beads switch between a inert and magnetisable states at rate $\kappa_s$ confirm that chromatin switching can yield microphase separation of epigenomic domains (Fig. \[fig:Switching\]B and Ref. [@PRE]). Besides switching, other ingredients may be at play, such as effective emulsification by RNA [@Hilbert2018], RNA-binding gel-forming proteins [@Nozawa2017; @Michieletto2019rnareview] epigenetic bookmarking [@Michieletto2018nar], and differential underlying transcriptional activity [@Smrek2017]. We note that copolymers made by heterochromatin and euchromatin would also generically microphase separate [@Leibler1980b; @Jost2014B; @Brackley2013pnas; @Scolari2015] in equilibrium – yet the epigenetic marks and the emerging domains would be static, unlike those found here. .

#### Conclusions – 

In summary, we have proposed a Landau-Ginzburg theory to describe the coupling between the dynamics of epigenetic states and that of the genome within the cell nucleus. By combining dynamical mean field theory and BD simulations, we find that this theory leads to a growth of competing epigenetic domains, in agreement with the well-known spreading of certain histone modifications [@Hathaway2012]. Our theory predicts several dynamical exponents that can be tested experimentally by triggering epigenomic reorganisations.

In equilibrium, one dominant epigenetic domain eventually takes over the whole nucleus. . We found that a generic nonequilibrium “epigenetic switch” between modifiable and refractory chromatin states (or first-order biochemical reaction) can arrest the unlimited spreading of epigenetic marks and stabilise a microphase separated state with coexisting epigenomic domains. We hope that our predictions will be tested in vivo and in vitro by controlling the rates of post-translational modification of certain protein complexes that are known to bind and/or remodel chromatin.

#### Acknowledgements – 

We thank the European Research Council (ERC CoG 648050 THREEDCELLPHYSICS) for funding. DMi and EO would also like to acknowledge the networking support by EUTOPIA (CA17139).

[^1]: Joint first author

[^2]: Joint first author
