Magnetic Polymer Models for Epigenomic Organisation and Phase Separation

The genetic instructions stored in the genome require an additional layer of information to robustly determine cell fate. This additional regulation is provided by the interplay between chromosome-patterning biochemical (“epigenetic”) marks and threedimensional genome folding. Yet, the physical principles underlying the dynamical coupling between three-dimensional genomic organisation and one-dimensional epigenetic patterns remain elusive. To shed light on this issue, here we study by mean field theory and Brownian dynamics simulations a magnetic polymer model for chromosomes, where each monomer carries a dynamic epigenetic mark. At the single chromosome level, we show that a first order transition describes the unlimited spreading of epigenetic marks, a phenomenon that is often observed in vivo. At the level of the whole nucleus, experiments suggest chromosomes form micro-phase separated compartments with distinct epigenetic marks. We here discover that for a melt of magnetic polymers such a morphology is thermodynamically unstable, but can be stabilised by a nonequilibrium and ATP-mediated epigenetic switch between different monomer states.

Introduction Each cell in our body contains the same DNA and hence it carries the same genetic information; yet, cells in different tissues possess distinct identities that are robustly inherited following multiple rounds of cell division [1,2].Thus, cellular fate cannot be directed by genetic cues alone and it requires an additional layer of information involving 3D genome organisation [3][4][5] and tissue-specific "epigenetic" patterns [6][7][8][9][10][11].The latter consist of biochemical tags that are deposited along the genome and on histones -the proteins in charge of packaging DNA into chromatin [1,2].The interplay between spatial genome organisation and epigenetic patterns guides the tissue-specific selection of which genes will be translated into proteins, in turn determining cellular identity [3,5,[12][13][14].One of the outstanding problems in biophysics is to understand how genome organisation and epigenetic patterns are linked to each other dynamically and what are the physical principles through which they regulate genome functionality and cellular memory [15][16][17][18][19][20][21].
To shed light on this issue here we introduce and study, both analytically and numerically, models of the genome where its 3D spatial organisation is coupled to a dynamically evolving epigenetic field [22].These models describe each chromosome as a magnetic polymer whose monomers encode (epigenetic) states which can change over time; they are therefore in the same universality class of annealed copolymers without global conservation laws [23].This model is markedly different from previous works on annealed copolymers with conserved number of elements in each state [24][25][26], and can be seen as a generalisation of the 1D Ising (or Potts) system where the substrate is allowed to diffuse in 3D space [15,16,27].We combine analytical mean-field theories with Brownian Dynamics (BD) simulations to simultaneously map the distribution of epigenetic marks and the 3D genomic arrangement within the cell nucleus.Together they describe the nuclear "epigenomic" organisation that can be directly compared with experiments [3,28].
At the single chromosome level, our magnetic polymer undergoes a first order transition between a swollen, epigenetically disordered fibre and a compact, epigenetically ordered one.Dynamically, the model generically predicts uncontrolled growth of the dominating epigenetic mark, reminiscent of the process through which transcriptionally repressed chromatin is often seen to spread in vivo, e.g., in X-chromosome inactivation [29,30] or positioneffect-variegation [6,31].At the whole nucleus level, a melt of magnetic polymers can initially phase separate into multiple thermodynamically metastable epigenetic domains; these though evolve into a single domain at large times.Introducing a local non-equilibrium epigenetic switch between an epigenetically active and an inactive state -mimicking ATP-dependent chromatin remodelling processes which modify chromatin accessibility to the deposition of biochemical marks -arrests the phase separation and ordering kinetics and yields microphase separation of the genome into multiple epigenetic domains, reminiscent of those observed in the cell nucleus [3,28,[32][33][34].
Single Chromosome To describe the equilibrium properties of a single chromosome fibre with a fluctuating epigenetic profile we consider an N -step self-avoiding walk (SAW) on a lattice with coordination number z where each vertex displays an epigenetic state q.The partition function of the model reads ∆ ri,rj J(q i , q j )   (1) arXiv:1807.11101v1[cond-mat.soft]29 Jul 2018 where 1/β = k B T and ∆ r,r = 1 if (r, r ) are nearestneighbour on the lattice (and 0 otherwise) thus restricting the interaction to 3D proximal segments.
For simplicity, we limit our model to three possible epigenetic states [9] (q i ∈ {−1, 0, 1}) and define J(q i , q j ) = − if q i = q j = ±1 and 0 otherwise.With this choice we implicitly assume that two marks are self-attractive (q i = ±1) while the third (q i = 0) is neutral (or unmarked) [9,15].This choice is also motivated by biological consideration, as we can assume that the non-neutral polymer states are associated with read-write protein complexes that can bridge polymer segments bearing the same epigenetic mark while "infecting" spatially neighbouring segments with the same mark [9,10].Both processes are captured by the same energy term and are akin to ferromagnetic interactions that align 3D proximal Ising or Potts spins, or bring them together when already aligned [23].
Eq. ( 1) can be solved within a mean field approximation [23,35] for an Ising-like model on a SAW (see SI for details).This approximation leads to the free energy density where α ≡ β z is the interaction parameter strength, ρ ≡ N/V the chromosome density (as the chain is confined into a box of volume V ) and φ is an epigenetic field (here modelling global epigenetic ordering).In analogy to ferromagnetic systems [23], we can identify φ as the average magnetisation of the system (see SI).By minimizing Eq. ( 2) with respect to ρ and φ one obtains the equilibrium phase diagram (see Fig. 1) where we distinguish two phases.At low α the system is in a swollen-disordered phase (SD): the chain is extended (ρ = 0) and heterogeneously coloured (φ 2 = 0).At large α we find a compact-ordered phase (CO) where the chain is crumpled (ρ = 0) and nearly uniformly coloured (φ 2 > 0).The discontinuous jumps of the order parameters ρ and φ at the transition point (α c 3.96) signal a first order transition between these two regimes [15,36].In Fig. 1 we also report snapshots of representative configurations from BD simulations of a corresponding polymer model where the Langevin dynamics of the polymer backbone is coupled to a Monte-Carlo annealing procedure that evolves the states of the polymer beads (see [15] and SI for details).The first order nature of the transition, as noted in [15], provides a mechanism to endow memory to a global epigenetic state.
A relaxation dynamics for ρ and φ can be written down starting from Eq. (2) in terms of two coupled "Model A" equations [37] as both fields are not conserved at the single chromosome level.[Here ρ should be understood as the density of beads within the smallest box containing the polymer chain].Such equations read where Γ ρ/φ and κ ρ/φ are mobilities and surface tensionlike coefficients, respectively.In Eqs.
(3) we decouple α into two independent parameters affecting the dynamics of the polymer (α ρ ) and of the epigenetic field (α φ ) separately.Note that the case α ρ = α φ leads to nonequilibrium dynamics as these equations no longer derive from a free energy.By numerically integrating Eqs.(3) we obtain the non-equilibrium phase diagram shown in Figure 1(b).We discover a new phase that is absent in equilibrium (α ρ = α φ ), featuring a crumpled and epigenetically disordered (CD) polymer.Yet, within the mean field approximation, we do not observe the swollen-ordered (SO) phase seen in BD [15,36].A single magnetic "epigenetic" polymer therefore exists in one of three phases in steady state, each reminiscent of a biologically relevant configuration.The SD phase models the conformation of a chromosome exiting mitosis, when epigenetic patterns and 3D folding are not yet established [3,38].The CO phase resembles the "Barr body" into which the inactive X-chromosome folds in female mammalian cells [29].This is a dense globular structure which is homogeneously marked with a repressive epigenetic state [30].Finally, the CD phase is akin to inert chromatin which experimental contact maps suggest is compact [32,39], yet has no clear epigenetic signature [39][40][41].
Our theory also offers a framework within which to understand the spreading of repressive marks (heterochromatin) in X-chromosome inactivation or in other position-effect-variegation where a transcriptionally silent domains spreads onto a nearby gene, switching off its expression [31,42].In our model the spreading occurs via a t 1/2 growth [37] when both epigenetic states are equally likely, whereas if one is favoured we expect linear Fisher-like growth [43].
Whole nucleus At the scale of the entire nucleus (volume V ) we assume that chromosomes are initially homogeneously filling the space.The overall density of chromatin n 0 ≡ N/V (where N is now the total length of the genome) is conserved, and the system can be described as a melt of magnetic polymers.
A minimal free-energy density describing the equilibrium properties of this model is where the fields n(x, t) and m(x, t) are the local chromosome density distribution and the average epigenetic marks respectively.The terms in Eq. ( 19) can be justified as follows: (i) the magnetisation field should not explicitly break its intrinsic Z 2 symmetry (if both marks are equally likely); (ii) the density field should be described by a standard virial expansion for non-ideal gases; (iii) the minimal coupling χm 2 n should capture the interplay between chromatin folding (n > 0) and epigenetic ordering (m 2 > 0).For convenience and without lack of generality we set a > 0, b > 0, c > 0, d > 0 and χ = χ(T ) > 0.
The equilibrium phase diagram (see Fig. 2) is obtained by first minimising Eq. ( 19) with respect to the non conserved field m, i.e. ∂f /∂m| m * = 0, and then by analysing the resulting f (m * , n) as a function of the conserved field n, via a common tangent construction [44,45].
For small values of χ the system is in a uniform (n = n 0 ) and epigenetic disordered phase (m = 0) (UD) (no epigenomic domains).Upon increasing the overall density n 0 (keeping χ ≤ χ c (n 0 ) fixed) we find a second order phase transition to a uniform state with ordered epigenetic field (m(x) 2 > 0) (UO) (see SI).The dynamics of the UD-UO transition is characterised by long-lived bicontinous spanning domains with alternated epigenetic marks (see Suppl.Movies), similar to growing magnetic domains in Ising systems [37,46].At large times these domains coalesce into a single system-spanning epigenetic domain (see Figs. S1,S2 and inset in Fig. 2 from the BD simulations).Finally, for χ > χ c (n 0 ), we observe that the uniform state is unstable and the system phase separates into high (n + ) and low (n − ) density regions forming a demixed-ordered (DO) phase.The high-density regions are associated with strong epigenetic domains (m 2 > 0) whereas the low-density regions with neutral epigenetic signature (m 2 = 0) (see Fig. 2).This phase is contained within the binodal curves which are determined using a common tangent construction [37,44,45] (see SI).We also mention that close to the critical point, where the binodal lines meet, the DO phase displays weaker variations of density throughout the system, i.e. the low density phase is strictly non-zero (0 < n − < n + ).We call this regime partially demixed ordered (PDO) phase (see Fig. 2).Pleasingly, the equilibrium phases obtained from the mean-field free-energy (19) are confirmed by BD simulations of a more realistic model in which the genome is described as a dense solution of magnetic polymers (see insets of Fig. 2 and SI).Some of the observed phases are reminiscent of the epigenomic organisation seen in experiments.The UD phase (as the SD phase for a single polymer) may represent a genomic configuration upon exit from mitosis, when spatial structure and epigenetic patterns are yet to be established (although our model does not account for mitotic chromosome structure).The (P)DO phase may be associated to strongly phase-separated nuclei, for instance in retinal [47] or senescent [48,49] cells.When quenching from the UD phase into the DO region, which may model the mitosis-interphase transition, the system phase separates into competing epigenomic domains which slowly evolve into homogeneously marked systems.These transient states display epigenomic organisations that are reminiscent of typical cell nuclei [28].Yet, the long-time steady state lacks epigenetic state coexistence and is fully phase separated, so is qualitatively different from typical nuclear organisation.The metastable multidomain state can be stabilised though, by driving the system away from equilibrium as detailed below.
Non-Equilibrium Epigenomic Organisation We now propose a non-equilibrium model for epigenomic organisation that can be derived starting from the free energy in Eq. (19).We consider its "model C" equations [37,46] and add two kinetic terms that dynamically convert the chromosomal density field from an "active" state, which can be biochemically marked (n a ) to an "inactive" one that is refractory to biochemical modification (n i ), and vice versa.This switch is inspired by the process of ATPdependent chromatin remodelling which changes local fibre structure and is coupled to histone modification [1].Note that now it is only the sum of the two density fields needs to be conserved at all times, i.e. n i + n a ≡ n 0 .The modified equations read where the parameters σ a/i describe the rates at which chromatin is activated/inactivated.In general we will consider σ i = σ a (see SI).We numerically evolve Eqs.(30) starting from the UD phase.Importantly, we find that the presence of non equilibrium switching terms now lead to arrest of both density phase separation and epigenetic ordering [50].The system stabilises into coexisting domains with high local density and non-zero epigenetic signature separated by regions with low active density (see Fig. 3).Large-scale BD simulations of magnetic polymer melts in which beads are switched from a passive "non-magnetisable" state to an active "magnetisable" one at rate κ confirm this phenomenology (see Fig. 3, and SI).

Conclusions
We have proposed and solved models of magnetic polymers that can be used to describe the coupling between epigenetic patterns and genome organisation both at the single chromosome and at the whole nucleus scale.
For a single chromosome, our magnetic polymer model can be solved at the mean field level [23] and displays three possible phases in steady state.The phase diagram is in agreement with that found from BD simulations [15], and the dynamics of the model generically entails uncontrolled spreading of the dominant epigenetic mark, which is reminiscent of epigenetic silencing dynamics in vivo [29].At the whole nucleus scale, we consider a Landau free energy density to describe the coupling between epigenetic states and chromosomal density.By combining dynamical mean field theory based on this free energy and direct BD simulations, we find that the model now leads to growth of many epigenetic domains with different marks, as found experimentally.In equilibrium one epigenetic domain eventually takes over the whole nucleus by spontaneous symmetry breaking.Unlimited spreading can though be contrasted by a non-equilibrium switching mechanism motivated by the phenomenon of ATP-dependent chromatin remodelling where each genomic segment can switch between a state in which it can be epigenetically marked and an inert one in which it cannot.
Our magnetic polymer model for epigenomic ordering can be extended in a number of ways.One is by introducing genomic bookmarking to seed domain formation [16].Another interesting avenue to explore would be to pursue a spin-glass model [51] instead of a Potts model for the underlying polymeric ordering.In this case, the rough free energy landscape of spin-glasses [52] might provide another avenue to stabilise a genome with micro-phase separated epigenetic domains.

Supplementary Information Single Chromosome Model
Here we obtain the mean field approximation presented in the main text to describe the thermodynamics of a single chromosome fiber with epigenetic marks.
Following Ref. [53], we describe the chromosome fiber as a N -steps self-avoiding walk (SAW) on a lattice with coordination number z.Each vertex of the walk carries an epigenetic state q that can assume three possible values (q ∈ {−1, 0, 1}).
Any pair of neighbouring (but non consecutive) vertices interact with each other via a contact potential that depends on their q-value.More precisely, if the i-th and the j-th vertices are nearest neighbours on the lattice, their contact energy J(q i , q j ) is with > 0. Note that the mark q = 0 does not contribute to this configurational energy and we will define it as a neutral mark.The equilibrium properties of this system is described by the following partition function where 1/β = k B T .The sums SAW and {q} run over the set of all N -steps SAWs and all the possible epigenetic states respectively.The matrix ∆ rirj is the adjacency matrix associated to a given SAW and is given by Notice that the partition function in Eq. ( 7) presents a clear Z 2 symmetry as J(−q i , −q j ) = J(q i , q j ).
Since we are here interested in the critical properties of the system, we can restrict the phase space of the epigenetic variables, q, to the case where the abundance of the state q = 0 is equal to the one of q = −1.With this restriction the system can be faithfully described by a two-valued spin variable S = {1, − 1  2 } where S = 1 corresponds to the mark q = 1, while the values S = − 1 2 has multiplicity 2 as it corresponds both to q = 0 and q = −1 [54].
By using the spin variable S, Eq. ( 6) becomes which can be re-written as and the partition function in Eq. ( 7) is then recast into Let us first evaluate, at a fixed γ ∈ SAW, the term: By using an Hubbard-Stratonovich transformation Eq. ( 12) becomes where dφ = N i=1 dφ i .By summing over all possible spin configurations we get This integral can be evaluated through an homogeneous saddle point approximation and by assuming the translational invariance of the field φ.This gives In general, the term N i,j=1 ∆ −1 rirj , depends on the given SAW and it is not easy to compute.However, it can be estimated if we restrict the set of SAWs to the ones that are almost space filling, i.e. ones that can be approximated as Hamiltonian walks [53].
An Hamiltonian walk is a path that visits each vertex of a lattice embedded in a volume V exactly once and have been used to study equilibrium properties of highly compact polymers [55,56].For an Hamiltonian walk, the adjacency matrix of the SAW ∆ takes the same form of the adjacency matrix of the underlying lattice and it is characterised by the coordination number z.Hence, N i,j=1 ∆ −1 rirj = N z .Here, we consider N -steps configurations that, similarly to Hamiltonian walks, are contained in a volume V but may in principle display a lower mean number of nearest neighbours, i.e. ρz instead of z.With this approximation N i,j=1 Notice that for generic SAWs with low ρ values Eq. 15 is not exact but is an upper bound.Finally, we evaluate the last term in Eq. ( 11), i.e.
By following the approach described in Ref. [57] we can approximate F SAW as By collecting all the terms and taking f = − T N log Z, we obtain the following mean-field free energy density where α ≡ β z.The equilibrium properties of the model are then obtained by minimizing Eq. ( 18) with respect to both, magnetisation φ and density ρ.As stated in the main text, this mean field approximation gives two possible equilibrium phases.For large values of α we find a compact-ordered phase (CO) where the chain is globular (ρ = 0) and nearly uniformly coloured (φ > 0).For small values of α the system is instead in a swollen-disordered phase (SD) where the chain is extended in space (ρ = 0) and heterogeneously coloured (φ = 0).At the transition point (α c 3.96) we observe a discontinuous jump of the parameters ρ and φ, proving the existence of a first order transition between the two phases [58,59].

Genome-wide Model
Here, we discuss the model we introduced in the main text to describe the equilibrium properties of epigenomic organisation at the scale of the full genome.By assuming that chromosomes fill a fixed volume V , we can define a conserved mean density n 0 ≡ N/V , where N is the total length of the genome.The equilibrium properties can be described by the following free-energy density where the fields n(x, t) and m(x, t) correspond to the local chromatin density distribution and the average epigenetic marks (or magnetisation) respectively.The phenomenological parameters of the uncoupled system are constant and set to be a > 0, b > 0, c > 0, d > 0.
The parameter χ > 0, governing the coupling between the epigenetic profile and the chromatin organisation, is temperature dependent.Since V is fixed, the local density n obeys the following constraint: The equilibrium properties are found by minimizing the functional F = V f (x) dx with the constraint in Eq. ( 20), i.e. constant n 0 .This is equivalent to find the minima of the functional G = F − µ V [n(x, t) − n 0 ] dx i.e. to solve the set of equations: where f [m(x), n(x)] denotes the free energy functional and δf /δm the functional derivative.By finding the solution to the first equation, i.e.
we restrict the problem to the effective free-energy den- ] that depends only on the conserved field n and reads: This procedure simplifies Eqs. ( 21) to the set of equations which is satisfied by the trivial uniform solution n(x) ≡ n 0 .In the non-trivial solution of these equations, instead, we find that in the system there is a cohexistence between two density phases n(x) = n − and n(x) = n + , have the same pressure P = f − n δf δn and chemical potentials µ, and are found via the so called common tangent construction [37].
Finally, the (spinodal) region in which the homogeneous solution n(x) = n 0 is unstable is characterised by < 0, which leads to Inside this region of values the homogeneous solution is linearly unstable and the system spontaneously demixes into low density (n − ) and high density (n + ) phases.By applying this procedure to the free-energy in Eq.( 19) we obtain the equilibrium phase-diagram as a function of the coupling parameter χ and genome density n 0 (see main text and Fig. 2).

Nature of the Phase Transitions
We now discuss the nature of the lines of phase transitions found in the equilibrium phase diagram: First, from Eq. ( 22) one can notice that the order parameter m goes continuously to zero.This strongly suggests that the transition from UD to UO is second order.Second, if a system is driven from the homogeneous phase, where n(x) = n 0 ∀x, to a region in which this solution becomes unstable, then it must cross a binodal line.At this point the pressure is the critical one P and we find that either lim {χ,n0}→P Similarly, if the system is driven from one demixed region (e.g.PDO with 0 < n − < n 0 < n + ) to another (e.g.DO with 0 = n − < n 0 < n + ), then the system must cross another point P where lim {χ,n0}→P n − = n − and lim {χ,n0}→P In light of this, and of the fact that the magnetic order parameter is continuous, we can conclude that every transition line in the phase diagram is continuous.Below we will focus in more detail on the transition from the homogeneous to the demixed phase, but a similar argument can be used for phase transitions between demixed phases.
A homogeneous phase displays F (eq) . By looking at the first derivatives of the free energy in Eq. ( 19), computed at the equilibrium, and using the above conditions we get This equality, together with the common tangent construction which gives the constraints ∂f , leads to: Since the order parameter n is continuous, if we drive the system from the homogeneous phase, we expect that either α → 0, Similarly, one can show that: Therefore, as the system pass from an homogeneous phase, to a demixed one, the first derivatives of the free energy are continuous.

Dynamical Scaling
Here we characterize the dynamical evolution of the system when it is quenched from a point within the UD phase into one within either the UO or the DO phase.This analysis may provide insights into the dynamics of the genome-wide spatial re-organisation and epigenetic recolouring, for example, at the beginning of interphase.In this model the density is a conserved order parameter while the magnetisation needs not be conserved.Hence, the dynamics of the system can be described by "model C" [37] equations: by using the free energy in Eq. ( 19) these become We numerically solve Eqs. ( 29) and monitor the time evolution of the density and magnetisation fields during several possible quenching trajectories in the phase space.We start from the UD phase and perform four representative quenches: (small n 0 ) (see inset of Fig. 4).
Following Q 1 , we observe that the density remains uniform while the epigenetic field coarsens into clusters of coherent colours which slowly evolve into one systemspanning domain through spontaneous symmetry breaking (see movie M1).The scaling of the typical epigenetic domain size grows as L(t) ∼ t α where α = 0.46 is compatible with Model A dynamics [46] (see Fig. 4).This is expected since the density field remains uniform.
We also observe that the other three quenches evolve on slower timescales as both fields need to be re-organised since we drive a transition from a homogeneous system to a demixed one (see Fig. 5).Specifically, for quenches Q 2 , Q 3 , and Q 4 , L m (t) ∼ t β with β 0.25 in agreement with previous results on Model C dynamics [46] (see Movie M2, M3, M4).

Non-Equilibrium Epigenetic Switching
Here we present the details of the non-equilibrium model for genome organisation with epigenetic switching.As reported in the main text, the dynamical equations are Eqs. (30) describe the dynamics of a "model C" [37] with two additional kinetic terms that dynamically convert the density fields from one that can be epigenetically marked (or active, n a ) to the one that is uncoupled from the epigenetic field (or inactive, n i ).As discussed in the main text, these terms may effectively account for the nonequilibrium action of so-called chromatin remodelling complexes [1] that render a local genomic region available for, or refractory to, epigenetic marking at a certain time.The amplitudes of σ a/i describe the rates at which the density fields (n a , n i ) are activated/inactivated, i.e. the rates at which chromatin remodelling factors act on the genome.One should notice that in this case the total density n a + n i must be conserved, i.e.
whereas n a and n i need not to be individually conserved.Nevertheless, since ∂(n a + n i )/∂t can be written as the divergence of a certain quantity, equation Eq. ( 31) is always satisfied.We also mention that by imposing a free energy of the form: where G(x, t) is a function such that: then Eqs.(30) can be derived from an effective free energy only if σ a = σ i .In this case G takes the form where G(x) is the Green function that solves the equation ∇ 2 G(x) = δ(x) and it depends on the system dimension.Note that, if the switching rates are equal, σ i = σ a , then the dynamical equations of the system can be understood as underlying an effective free energy, thus entailing that the system is in equilibrium.On the other hand, the general condition that σ a = σ i , entails that Eqs.(30) describe a purely non-equilibrium system.

Steady States of the Switching Model
We now study the dynamics and the steady states of the model described by Eqs.(30) varying the values of σ i Figure 5. Representative snapshots of the system at different times of its evolution for two quenches.On the left-hand side, the quench Q1 brings the system from an UD phase to the UO phase: while the density field remains homogeneous, the magnetisation shows the formation of clusters which slowly evolve into a uniformly coloured system.Here the dynamics is akin to the Model A one.On the right-hand side, Q3 brings the system from the UD phase to the DO phase: in this case both magnetisation and density fields needs to re-organize and the system evolves on slower time-scales.and σ a .We keep the phase diagram of the system (Fig. 2 of main text) as a reference and fix the values of n 0 and χ such that a phase in the limit of negligible density of inactive marks, i.e. σ a σ i can be observed.If we quench the system either into the Uniform Ordered or the Demixed Ordered phases, then by varying σ i and σ a leads to a new stationary state similar to the Partially Demixed Ordered phase, i.e. one characterised by weak variations of the total density (n − > 0) and denoted by a non-null magnetisation m 2 > 0 (see movie MS1, movie MS2, movie MS3, movie MS4).
The non-equilibrium phase diagrams of the model at fixed n 0 and χ as a function of the two kinetic rates σ a/i are shown in Fig. S6.In most of the cases, these pictures show that the ordered phases arise when the fraction σ i /σ a is lower than a certain critical ratio r which can be estimated as follows: in steady state, Eqs.(30) predict a mean active density On the other hand, in Eq. ( 25) we have shown that the ordered states are stable only if the active density n a > a χ .Thus, one can conclude that the Ordered phases (Uniform or Demixed) are strongly favoured if in very good quantitative agreement with the observations from the numerical evolution of the system (see Fig. S6 black lines).
Finally, the free parameter ε(q i , q j ) is: and N is a parameter which ensures that the minimum of the attractive part is − k B T L .
3. U H describes the connection between consecutive beads along the chain: where k h models the connectivity strength and it is set ot kH = 200 .
We then use these potentials to evolve the equations of motion for each bead in the system using a fixedvolume and constant-temperature molecular dynamics (MD) simulations (NVT ensemble).The simulations are within the LAMMPS engine [64] and the equations of motion are integrated using a velocity Verlet algorithm, in which all beads are weakly coupled to a Langevin heat bath with friction γ = τ −1 B where τ B = 3πησ 3 /k B T is the self-diffusion (Brownian) time of a bead of size σ moving in a solution with viscosity η (which we consider water, i.e. η = 1cP , for the mapping to real units).Finally, the integration time step is set to ∆τ = 0.01 τ B .
As mentioned before, "recolouring" steps are performed every τ R = 100τ B and in each step we attempt a number of moves equal to the number of beads in the system.In each move, we randomly select a bead and randomly change its colour to a different one.If the move lowers the energy of the system we accept it, otherwise we assign an acceptance probability p = e −∆E/k B T R where ∆E is the change in system energy after and before the move.
In this scheme, it is straightforward to implement nonequilibrium switching by defining a fourth bead type (or q = 2) which does not participate to the recolouring dynamics, i.e. beads bearing q = 2 are excluded from the recolouring moves.Then, at rate σ i , beads bearing q = {−1, 0, 1} are randomly converted into q = 2 and viceversa at rate σ a .

Single Chromosomes
We employ single chromosome BD simulations of this model to confirm the results obtained through our continuum model in Fig. 1 of the main text.In the equilibrium case (T L = T R ) the main parameter that is varied to confirm the phase diagram reported in Fig. 1a is α = ε/k B T L .In the non-equilibrium case, we break detailed balance and independently vary T L and T R while maintaining ε = 1.Our results are robust with respect to the choice of recolouring rate τ −1 R and initial conditions.

Full Nucleus
To model the whole nucleus we perform simulations of a melt of annealed polymers at different monomer densities ρ = N/V and ε/k B T L .We here consider N = 50 polymers with M = 256 beads each and the range of parameters employed are ρ = 0.1 -0.8 σ −3 and ε/k B T L = 0.75-1.1.The insets of Fig. 2 in the main text are obtained using the following parameters:

Captions of Supplementary Movies
• Movie M1: Time evolution of the system described by eqs.( 29), after quench and n 0 = 2).The system is initialised in a UD phase (homogeneous density, incoherent magnetisation), and evolves towards a UO phase, where the system is still homogeneous, but the magnetisation is organised in big clusters of coherent magnetisation.
• Movie M2: Time evolution of the system described by eqs.( 29), after quench Q 2 (Γ m = Γ n = D m = D n = 1, χ = 3.5, and n 0 = 0.5).Here the system is initialised in a UD phase (homogeneous density, incoherent magnetisation), and evolves towards a PDO phase, where the system organizes in clusters, and it is characterised by weak density variations.
• Movie M3: Time evolution of the system described by eqs.(29), after quench Q 3 (Γ m = Γ n = D m = D n = 1, χ = 6, and n 0 = 2).The system is initialised in a UD phase (homogeneous density, incoherent colouring), and evolves towards a DO phase, where the system organizes in clusters and it is characterised by strong density variations.
• Movie M4: Time evolution of the system described by eqs.(29), following the quench Q 4 (Γ m = Γ n = D m = D n = 1, χ = 6, and n 0 = 0.2).The system is initialised in a UD phase (homogeneous density, incoherent magnetisation), and evolves towards a DO phase, where the system organizes in clusters and it is characterised by strong density variations.Compared to the Movie M3, the clusters appear to be smaller.
• Movie M5: BD simulations of a melt of magnetic annealed polymers with monomer density ρ = 0.5σ −3 and ε/k B T L = 0.7.The systems is initialised with a random colouring and it evolves towards a uniform ordered state where the large majority of beads are red via spontaneous symmetry breaking.
• Movie MS1: Numerical integration of eqs.(30), with parameters Γ m = Γ n = D m = D n = 1, χ = 4, n 0 = 2, σ a = 0.1, and σ i = 0.5.The system is initialised in a UD phase (homogeneous density, incoherent magnetisation).We observe that both the active and inactive densities organize in patterns similar to the ones observed in the PDO phase at the equilibrium.Remarkably, while the magnetisation in the equilibrium PDO phases was negligible in the low density regions, here it assumes a positive (or negative) value that is consistent with the neighbouring high-density areas.
• Movie MS2: Numerical integration of eqs.(30), with parameters Γ m = Γ n = D m = D n = 1, χ = 6, n = 2, σ a = 10, and σ i = 2.The system is initialised in a UD phase (homogeneous density, incoherent magnetisation).Both the density fields, and the magnetisation field show a behaviour similar to the one observed in the equilibrium DO phase.
• Movie MS5: BD simulations of a melt of magnetic annealed polymers with monomer density ρ = 0.8σ −3 and ε/k B T L = 0.9 and switching rate κ = 10 −4 τ B .This Movie shows that the evolution towards a uniformly coloured state is arrested and epigenomic (epigenetic and density) domains appear.For simplicity we only show the beads that are either red or blue (q = −1, 1) and not the neutral or inactive types.
• Movie MS6: BD simulations of a melt of magnetic annealed polymers with monomer density ρ = 0.8σ −3 and ε/k B T L = 0.9 and switching rate κ = 10 −5 τ B .Compared with Movie MS5, the domains appear larger.For simplicity we only show the beads that are either red or blue (q = −1, 1) and not the neutral or inactive types.

Figure 1 .
Figure 1.Phase Diagrams at the single chromosome scale.(a) Equilibrium phase diagram of a magnetic "epigenetic" polymer described by the free energy in Eq. (2).The system undergoes a first order transition (marked by a discontinuity in the order parameters) at α = /kBT 3.96 between a swollen-disordered (SD) and a compact-ordered (CO) phase.(b) Non-equilibrium phase diagram obtained integrating Eqs.(3) in the parameter space (αρ, α φ ).Insets: snapshots of representative configurations from BD simulations (see SI for details).

Figure 2 .
Figure 2. Phase Diagram at the Nuclear Scale Equilibrium phase diagram of a melt of magnetic polymers obtained by using the common tangent construction on the free energy Eq.(19).The three equilibrium phases are: (UD) uniform (n = n0) and epigenetically disordered (m 2 = 0); (UO) uniform (n = n0) and epigenetically ordered (m 2 > 0); (DO) demixed and epigenetically ordered (n = n+, m 2 > 0 and n = n−,m 2 = 0).A fourth partially-demixed ordered (PDO) phase is characterised by weaker variations in density (n− > 0) and denoted by a white shading within the DO phase.The dotted line marks the critical value of the coupling χc(n0), the solid lines identify the boundaries of the coexistence region (binodals) and the dashed lines identify the spinodal region where the uniform solution is linearly unstable[37,44].Insets report representative snapshots from Brownian Dynamics simulations of dense solution of magnetic polymers (see SI for details).

Figure 3 .
Figure 3. Non-Equilibrium Switching Drives Arrested Phase Separated Epigenomic Domains.(a) Snapshot from a BD simulation of a melt of polymers with switching rate κ = 10 −4 τB (τB is the diffusion time of a monomer), monomer density ρ = 0.8σ −3 and ε/kBTL = 0.9 (see SI for details).(b) Snapshot of a steady state configuration obtained evolving Eqs.(30) with parameters Γm = Γn = km = kn = 1, χ = 6, n0 = 0.5, σa = σi = 5.(c) Evolution of typical epigenetic domain size as a function of time and for different switching rates.These figures show that the ordering dynamics is arrested and domains with well-defined (self-limiting) size are formed when epigenetic switching is included in the model.

Figure 4 .
Figure 4. Growth of Epigenetic Domains.Evolution of typical epigenetic domain size L following a quench from the uniform disordered phase.The growth displays a power law that is compatible with Model A dynamics (α = 1/2) when the density field remains uniform (quench Q1), but it is significantly slowed down when both fields are re-organised (α = 1/4).In the inset we schematically show the paths of the quenches in the phase diagram.Here we evolved Eqs.(29) with fixed a = b = c = d = 1 and Γm = Γn = κm = κn = 1.