In interphase nuclei, chromatin forms dense domains of characteristic sizes, but the influence of transcription and histone modifications on domain size is not understood. We present a theoretical model exploring this relationship, considering chromatin-chromatin interactions, histone modifications, and chromatin extrusion. We predict that the size of heterochromatic domains is governed by a balance among the diffusive flux of methylated histones sustaining them and the acetylation reactions in the domains and the process of loop extrusion via supercoiling by RNAPII at their periphery, which contributes to size reduction. Super-resolution and nano-imaging of five distinct cell lines confirm the predictions indicating that the absence of transcription leads to larger heterochromatin domains. Furthermore, the model accurately reproduces the findings regarding how transcription-mediated supercoiling loss can mitigate the impacts of excessive cohesin loading. Our findings shed light on the role of transcription in genome organization, offering insights into chromatin dynamics and potential therapeutic targets.

The three-dimensional organization of chromatin within the nucleus is key to understanding the biophysical origin of critical cellular activities ranging from cell fate decisions to migration, proliferation, and metabolism. The existence of a multiscale chromatin organization has been observed not only from sequencing and contact-mapping techniques 1,2 , but also super-resolution imaging [3][4][5][6][7][8] . At the microscale, chromatin is organized into transcriptionally distinct compartmentsa transcriptionally active, loosely packed euchromatin phase and a tightly packed, predominantly silent heterochromatin phase. Finer resolution of the chromatin conformation reveals the existence of a more detailed spatial organization ranging from self-interacting topologically associated domains (TADs) to chromatin loops, a feature of the chromatin polymer mediating direct contact between gene regulatory elements bound by the CCCTC-binding factor (CTCF) and the cohesin complex 1,2 . The chromatin fibers can be trapped in and pushed through cohesin rings via a process called loop extrusion, until either CTCF bound sites are encountered or cohesin is unloaded [9][10][11] . In addition to direct extrusion of DNA loops via cohesin motor activity [12][13][14][15][16][17][18][19] , RNA polymerase II (RNAP II), a protein complex essential for DNA transcription, has been identified to play a significant role in enabling the movement of the chromatin fiber resulting in chromatin loop extrusion through cohesin 19,20 . Specifically, by altering the DNA winding thereby supercoiling it, transcriptional activity has been proposed to play a role in in-vivo chromatin loop extrusion 19,20 . A recent experimental study combined super-resolution imaging of chromatin and single molecule tracking of cohesin with various biological perturbations, such as pharmacological and genetic inhibition of transcription, supercoiling, and loop extrusion. This approach provided compelling evidence that transcription-mediated supercoiling regulates loop extrusion, as well as the spatial organization of chromatin within the nucleus 19 . These observations present a promising avenue of crosstalk between chromatin's multiscale structural organization and its transcriptional status. This indicates that a bi-directional coupling exists, such that not only do the distinct phases of chromatin organization regulate transcription, but transcriptional activity can also affect genome organization via chromatin tethering, extrusion, and decompaction 19,21 . While the local microscopic effects of transcription on spatial DNA organization have been previously investigated, a fundamental quantitative understanding of the physical mechanisms involved in the global genomic organization, due to transcriptional and epigenetic regulation, is not yet fully understood.

Here, we propose a mesoscale coarse-grained, polymer physicsbased mathematical model to capture the formation of chromatin domains while incorporating the spatiotemporal role of transcriptiondriven chromatin extrusion kinetics. Chromatin-chromatin interactions establish an energy landscape which drives a separation of hetero-and euchromatin phases. The dynamics of this evolution are governed by the diffusion of nucleoplasm and epigenetic reactions. Such evolution leads to the formation of functionally distinct heterochromatin domains of characteristic sizes. Chromatin-lamina interactions along the nuclear periphery give rise to lamina associated heterochromatin domains. The supercoiling-driven chromatin loop extrusion through active transcription is captured via the conversion of inactive compacted heterochromatin into transcriptionally active euchromatin loops along the chromatin phase boundaries. Essential and unique to our model is the interplay of the epigenetic and transcriptional kinetics in governing meso-scale chromatin organizationincluding the size of heterochromatin domains and their spacing in the interior and periphery of the nucleus.

Using this model, we make quantitative predictions that offer a mechanistic explanation for the emergence of size scaling of compacted heterochromatin domains with the rate of supercoilingmediated loop extrusion at the domain interfaces. Importantly, by including the interactions of chromatin with the nuclear lamina, we show the quantitative dependence of the sizes of lamina-associated domains (LADs) as well as those of interior chromatin domains on the level of transcriptional activity. The predictions on the size scaling of heterochromatin domains made by the model are agnostic to specific interactions, and thus are not limited to a particular cell type. Indeed, the model predictions are qualitatively validated experimentally on five different cell lines and using two different nanoscopic imaging approaches. We used partial wave spectroscopy (PWS), which enables high-throughput, label-free, live cell imaging, in conjunction with scanning transmission electron microscopy tomography with Chro-mEM staining (ChromSTEM), which allows 3-dimensional high-resolution quantification of chromatin mass distribution, to quantify statistical domain properties upon inhibition of transcription. We, further, quantitatively validated our predictions by analyzing the length scales of compacted chromatin domains previously reported using stochastic optical reconstruction microscopy (STORM) imaging 19 . In conjunction with super-resolution microscopy and nanoimaging techniques, our model establishes a foundation for a predictive framework with broad implications for understanding the role of transcriptional and epigenetic crosstalk in defining mesoscale genome organization.

We have developed a mathematical model to capture dynamic chromatin organization in the nucleus, in terms of its compaction into the heterochromatic phase or decompaction into the euchromatic phase (Fig. 1a). We treat the meso-scale genomic organization as a dynamic, far-from-equilibrium process, governed by the energetics of phaseseparation in conjunction with the kinetics of epigenetic reactions and the formation of chromatin loops aided by supercoiled DNA extrusion through cohesin due to RNAPII-mediated transcription. The model ingredients are depicted schematically in Fig. 1a. We begin by defining the energetics of the chromatin distribution in terms of the entropicenthalpic balance of chromatin-chromatin interactions, the chromatinlamina interactions as well as the penalty on the formation of phase boundaries via Eq. ( 6) (refer Methods, and Supplementary Section S1.2 in the SI). The gradients in the free-energy landscape, defined as the chemical potential (refer Supplementary Eq. ( S3)), drive the dynamic evolution of chromatin towards the two energy wells corresponding to the euchromatin and heterochromatin phases via Eq. (7a, b) (refer Methods, Supplementary Section (S1.4) in the SI). Interconversion of the two phases of chromatin can occur via (a) epigenetic regulation of histone acetylation and methylation (Fig. 1b), and (b) supercoiling-driven extrusion of chromatin loops from heterochromatin into euchromatin along the phase boundaries (Eq. (7b)) as shown in Fig. 1c.

The process of phase separation is initiated by adding a random perturbation to the initially uniform chromatin configuration (as shown in Fig. 2a, left panel) which captures the intrinsic intranuclear heterogeneities. As the simulation progresses heterochromatin domains (in red, center panel of Fig. 2a) spontaneously nucleate and grow. The evolution ultimately stabilizes resulting in a steady state (right panel of Fig. 2a) with a quasi-periodic distribution of stable domains of heterochromatin rich phase (ϕ h = ϕ max h ) in red and euchromatin rich phase ϕ h = 0 À Á in blue. Each of these domains are nearly circular (see Supplementary Section S2 of SI for a discussion on non-circular lamellar domains) with characteristic sizes. Concomitantly, heterochromatin domains localized to the nuclear lamina (called LADs) of comparable sizes appear in our simulations (Fig. 2a).

The meso-scale distribution of chromatin throughout the nucleus predicted by the mathematical model presents a striking qualitative similarity with the experimentally observed distribution of DNA in the nucleus using ChromSTEM, and STORM as reported previously 19 (Fig. 2b). Domains of compacted chromatin with a characteristic size are observed via a high histone density distinguished from regions of low histone density (Fig. 2b). Lastly, the preferential accumulation of heterochromatin domains along the nuclear periphery seen via STORM imaging (Fig. 2b), again with similar size scale, is also in excellent agreement with the experiments.

When defining the free energy density of chromatin organization in the nucleus (see Supplementary Eq. (S1) in SI), we penalized the formation of sharp interfaces via an interface penalty η, defined as the energy cost associated with the formation of the interfaces between heterochromatin and euchromatin phases. As we show in the SI (Supplementary Section S1.5), the energy penalty η results in the formation of a smooth rather than a sharp interface between the heterochromatin and the euchromatin phases. Numerical simulations of chromatin organization exhibit such smooth interfaces around chromatin domains, as shown in the zoomed in image in Fig. 2c (right panel). The width of the interface δ is controlled by the competition between the interfacial and bulk energy contributions (refer Supplementary Section S1.5).

Smooth chromatin phase boundaries are indeed observed in-vivo via Chrom-STEM imaging (Supplementary Section S1.11). We characterized the 3D chromatin density around individual heterochromatin domains in a BJ fibroblast nucleus using Chrom-STEM (Fig. 2c, left panel; Supplementary Fig. S5). We estimated the average chromatin density within concentric circles emerging from the center of individual domains to the periphery (Fig. 2c, Supplementary Fig. S5). The chromatin density was highest at the core of the domain and dropped slowly from the center of the domain to the periphery. The smooth decrease in radial density indicates that the chromatin domain boundaries are not abrupt (Fig. 2c), in agreement with the numerical simulations.

We next investigate how the size scaling of the heterochromatin domains is regulated by the epigenetic reactionsacetylation and methylation of histonesand supercoiling-driven chromatin extrusion which together can lead to interconversion between heterochromatin and euchromatin. First, we see that in the absence of the epigenetic as chromatin-chromatin interactions, the chromatin-lamina interactions and epigenetic regulation. Additionally, extrusion of chromatin loops due to DNA supercoilingwhich is increased by transcriptional activityalso plays a role in mesoscale genomic organization. While this may occur within either chromatin phases (red circle), we further explore the role of chromatin loop extrusion at the heterochromatin-euchromatin interface (black circle). b The model captures the chromatin-chromatin interaction energetics via a double well free energy description as shown in the contour plot. The two wells correspond to the heterochromatin (red circle) and euchromatin phases (blue circle). Any initial configuration (light blue circle) spontaneously decomposes into these wells at steady state. The dynamics of this transition are governed by diffusion and reaction kinetics comprising of epigenetic regulation and kinetics of supercoiling-driven chromatin extrusion (red box inset). c Loading of cohesin assisted by NIPBL/MAU2 initiates the formation of chromatin loops. Cohesin can also be dynamically unloaded via unloading factors viz. WAPL/PDS5. Active processes such as RNAPII mediated transcription further drive the extrusion of trapped DNA, supercoiling it into chromatin loops.

reactions and chromatin extrusion multiple domains of a characteristic size are not obtained as shown in Supplementary Fig. S10 (detailed discussion in Supplementary Section S5). In this case, although nucleation of multiple heterochromatin domains occurs even without reactions (Supplementary Fig. S10a), all of them merge into a single large cluster driven by Ostwald ripening so as to minimize the interface formation.

The model also predicts that the size of the heterochromatin domains in the interior and periphery can be regulated by the epigenetic reaction rates of acetylation and methylation as shown in Supplementary Fig. S6 (Supplementary Section S2). We see that as methylation increases the size of the interior domains increases too. On the other hand, increase in acetylation results in the formation of smaller heterochromatin domains. The trends followed by the domains towards the interior of the nucleus are replicated by the LADs as well. Lastly, we identify that the size scales of the domainsthe domain radii in the interior of the nucleus and the LAD thickness along its peripherydepend on the level of transcription governed supercoiling-driven chromatin extrusion rate e Γ a (Fig. 2d, Supplementary Fig. S6). We note that, as the transcription ( e Γ a ) is increased, the sizes of the heterochromatin domains decrease, both in the interior as well as at the periphery. At the same time, we also note that as chromatin extrusion rate is increased, the average volume fraction of heterochromatin ϕ h À Á in the nucleus decreases, while that of euchromatin ϕ e À Á increases.

Next, we theoretically predict an explicit dependence of the sizes of interior heterochromatic domains and LADs on epigenetic and transcription reactions and the diffusion kinetics of the epigenetic marks.

Intuitively, in the presence of more repressive methylation the overall heterochromatin content in the nucleus should increase, while in higher histone acetylation conditions the overall euchromatin content will increase. Thus, the epigenetic reactions can independently determine the average volume fractions of each form of chromatin, thereby breaking the detailed balance condition where the free energies of each phase determine their relative abundance in a thermodynamic equilibrium. A mathematical relation between the average volume fraction of each chromatin phase and the epigenetic reaction parameters can be determined by averaging the chromatin evolution equation (Eq. (7b)) at a steady state (i.e. ∂ϕ d ∂ e t = 0). In the absence of transcription driven chromatin extrusion (i.e. e Γ a = 0), we see that the epigenetic kinetics regulates the average heterochromatin content of the nucleus as,

(Supplementary Eq. (S23), refer Supplementary Section S3 for more details). The presence of transcription-mediated loop extrusion kinetics (i.e., e Γ a ≠ 0 in Eq. (7b)) further augments the deviation from thermodynamic equilibrium (i.e., the breaking of detail balance) via surface reactions that actively extrude DNA at the interface of heterochromatic domains. In the presence of transcription, the average heterochromatin (and euchromatin) content in the nucleus becomes (refer Supplementary Eq. (S22)),

where κ is a function of ϕ max h , volume fraction change across the interface Δϕ, and the length of the interface between the two chromatin phases (refer Supplementary Section S3 for derivation). Since supercoiling-mediated chromatin extrusion converts the tightly packed heterochromatin into low density transcriptionally active euchromatin phase, as extrusion rate e Γ a increases, the average heterochromatin content decreases.

Thus, the overall mean chromatin composition of the nucleus ϕ h , ϕ e À Á is determined by the reaction kinetics of epigenetic regulation along with transcription. The reaction kinetics alone would drive a homogenous chromatin organization with ϕ h , ϕ e À Á . On the ϕ d ,ϕ n À Á phase space we see that the average composition (shown as a light blue circle in Fig. 1b) determined by reactions is energetically unfavorableit does not lie in the energy wellsand hence must evolve in time.

Next, we show that the average composition of the two chromatin phases, shown in Fig. 2e(i), plays a key role in the emergence of the characteristic sizes of the heterochromatin domains. To illustrate this, we first observe that the mean chromatin composition ϕ h , ϕ e À Á lies in neither of the energy wells as shown in Fig. 1b (light blue circle) and is thus energetically unfavorable. The need to reduce the total free energy in the nucleus drives the system to phase separate by nucleating heterochromatin domains (Fig. 2e(iii)) corresponding to the red energy well labeled heterochromatin in Fig. 1b surrounded by euchromatin domains corresponding to the dark blue energy well labeled euchromatin. The events entailing the individual steps in the nucleation and growth of a single droplet of heterochromatin due to phase separation, as shown in Fig. 2e, are as follows:

1. Due to phase separation, the heterochromatin volume fraction immediately outside the droplet is ϕ h = 0 corresponding to the euchromatic energy well. Far away from the droplet, the mean composition ϕ h , ϕ e À Á remains undisturbed. The resulting spatial gradient in the chromatin composition (blue curve in Fig. 2e(iv)) sets up a diffusive flux of heterochromatin into the droplet, allowing it to grow. 2. On the other hand, within the heterochromatin droplet (with ϕ h = ϕ max h ) histone acetylation reactions will allow conversion of heterochromatin inside the droplet into euchromatin outside. Active supercoiling-mediated chromatin loop extrusion further adds to the heterochromatin outflux. Together loop extrusion and acetylation oppose the diffusive influx of heterochromatin and thereby reduce the size of the droplet (Fig. 2e). 3. Based on the above observations, the rate at which the nucleated heterochromatin droplet grows can be written in terms of the balance of reaction-diffusion gradient driven influx and acetylation and transcription driven outflux of heterochromatin as (refer Supplementary Section S4, Supplementary Eq. (S25)),

4π e R 2 d d e R d d e t = 4π e R d ϕ h |fflfflfflffl ffl{zfflfflfflffl ffl} inwards diffusion À 4 3 π e R 3 d ϕ max h |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} Acetylation working against inwards diffusion À 4π e R 2 d δ 2 e Γ a ϕ max h |fflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflffl ffl} Chromatin extrusion working against inwards diffusion

where δ is the rescaled width of the interface, which is in turn related to the length scale obtained via the competition between the interfacial energy and chromatin-chromatin interaction (refer Supplementary 2 | Numerically predicted chromatin distribution in the nucleus captures the salient features of in-vivo chromatin organization. a Visualization of the chromatin organization obtained from the simulations. The initial chromatin organization is a homogenous distribution with a small perturbation added, resulting in nucleation of heterochromatin domains (center panel) which grow into heterochromatin domains of characteristic sizes at a steady state. b Superresolution visualizations of chromatin organization observed in-vivo via STORM imaging of HeLa nuclei (left panel, scale bar 3 μm, data previously reported in ref. 19, n = 19 nuclei) and ChromSTEM imaging of BJ fibroblast nuclei (right panel, scale bar 1 μm, n = 1 nucleus) show that chromatin organization in nucleus is characterized by interspersed heterochromatic domains of comparable sizes. c The smooth boundaries of the chromatin packing domains as seen in ChromSTEM observations are captured by the model. d Numerically predicted trend of sizes of heterochromatin domains as the transcription-mediated chromatin extrusion rate increases. e Schematic diagrams of the step-by step events (events 'i' through 'vi') involved in the nucleation, growth and stabilization of heterochromatin domains at a steady state. f Plot of theoretically evaluated growth rate of heterochromatin domains with (red) and without (blue) reactions. Reactions give rise to a stable domain radius. In the absence of reactions, no stable heterochromatin domain length scales are observed. g The evaluation of stable radius (blue) and stable LAD thickness (red) as transcription mediated surface reactions are changed. Here, the relative radius is defined as the steady state radius relative to its value when transcription is zero, i.e., relative radius = e R SS d = e R SS d jΓ a = 0. The relative LAD thickness is similarly defined. Section S1.5). The resulting evolution of the droplet growth rate d e R d =d e t as the radius of the droplet increases is shown in Fig. 2e. Notice the two fixed points (Fig. 2f, labeled critical and stable radius) where d e R d =d e t = 0. Beyond the critical radius the domains grow in size. 4. The second fixed point (stable radius) corresponds to the rescaled steady state (i.e., d e R d =d e t = 0) heterochromatin domain size as determined by the active epigenetic and the transcriptional regulation in tandem with passive diffusion, and can be written as (derivation shown in Supplementary Section S4, Supplementary Eq. (S27)), e R ss d = À 3 e Γ a δ 4 + ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 3 e Γ a δ 4 2 + 3 ϕ max h e Γ me 1À

From Eq. ( 3), we observe that the steady state droplet radius e R ss d

depends on both diffusion and reaction kinetics. With increase in methylation, e R ss d increases implying bigger heterochromatin domains. On the other hand, with increase in either the acetylation or transcription-mediated loop extrusion the steady state radius decreases. The quantitative dependence of the steady state radius on transcriptional kinetics is shown in Fig. 2g (blue solid line). Note that the steady state radius shown in Fig. 2g is normalized relative to the steady state radius with no transcription. Thus, our theory predicts an increase in the sizes of compacted chromatin domains in the interior of the nucleus upon inhibition of transcription.

The size dependence of chromatin domains along the nuclear periphery can be similarly determined by the balance of reaction, transcription, and diffusion kinetics for the LADs. The affinity of chromatin to the nuclear periphery due to the chromatin-lamina interactions in Eq. ( 6) induces a preferential nucleation of LADs. A schematic representation of heterochromatin compaction along the nuclear periphery resulting in LAD growth is shown in Fig. 2e. As with the interior heterochromatin droplet, phase-separation drives the heterochromatin compaction ϕ h = ϕ max h À Á within the LADs, while the chromatin immediately outside corresponds to the euchromatin energy minimal well ϕ h = 0 À Á . Far away from the peripheral LAD nucleation sites, the chromatin composition remains undisturbed at the average composition of ϕ h , ϕ e À Á . The variation of chromatin composition with distance from nuclear periphery is shown in Fig. 2e (blue line). Like in the case of the interior heterochromatin droplets, the heterochromatin composition gradient driven diffusive influx is balanced by the epigenetic and transcriptional regulated heterochromatin outflux, which determines the rescaled steady-state thickness of the LADs (refer to the Supplementary Section S7, Supplementary Eq. (S34)),

As with the interior domains, we observe that the LADs become thicker with increase in methylation, while they become thinner with increasing acetylation or chromatin extrusion rates. A quantitative dependence of steady state LAD thickness on transcription rate based on Eq. ( 4) is plotted in Fig. 2g (red dashed line). Our theory predicts an increase in the sizes of LADs along the nuclear periphery upon inhibition of transcription. While the theoretical analysis helps develop a fundamental biophysical understanding of the role of energetics and kinetics in chromatin phase separation, a nucleus-wide chromatin organization and its dynamic evolution can only be obtained numerically.

Next, we use the in-silico model to make testable quantitative predictions of the meso-scale chromatin organization in the nucleus. We also report the in-vivo nuclear chromatin reorganization upon transcription inhibition using complimentary STORM 19 and ChromSTEMon nuclei from multiple cell lines. The choice of the parameters for rates of acetylation e Γ ac , methylation e Γ me , and the strength of chromatinlamina interactions e V L , were held constant for all the following simulations, and the choice of the level of spatial noise is discussed in the Supplementary Section S8. We calibrate the active chromatin supercoiling-driven loop extrusion rate Γ a to obtain an in-silico change in the interior domain sizes quantitatively comparable to that observed upon transcriptional inhibition. The calibrated model is then used to predict the change in LAD thickness due to inhibition of transcription, which upon comparison with experimental images serves to validate the model. A schematic for the workflow utilized to calibrate and cross-validate the model predictions in the interior and along periphery of the nucleus is shown in Supplementary Fig. (S14) (Supplementary Section S8).

ChromSTEM was used to obtain super-resolution images in terms of statistical descriptions of chromatin packing domains for BJ fibroblasts. ChromSTEM allows the quantification of 3D chromatin conformation with high resolution 22 . ChromSTEM mass density tomograms were collected for BJ fibroblasts treated with Actinomycin D (ActD) (Fig. 3a, center) and compared to DMSO treated mock controls (Fig. 3a, left) to evaluate the average size and density of chromatin packing domains. We have previously demonstrated that chromatin forms spatially well-defined higher-order packing domains and that, within these domains, chromatin exhibits a polymeric power-law scaling behavior with radially decreasing mass density moving outwards from the center of the domain 23 . As the ChromSTEM intensity in the reconstructed tomogram is proportional to the chromatin mass density, we estimated the size of the domains based on where the chromatin mass scaling and the radial chromatin density deviate from their predicted behavior (discussed in Supplementary Section S1.11). Based on the statistical analysis of individual packing domains, in a single tomograph shown in Fig. 3a, we observed 71 domains in DMSO and 48 domains in the ActD-treated nucleus. Of the identified domains, the average domain radius ( ± S.E) of BJ cells treated with DMSO and ActD was estimated to be 103.5 ± 4.73 nm and 129.7 ± 6.78 nm, respectively (Fig. 3a, right panel), representing a 20.2% increase in size. Overall, fewer domain centers, and larger chromatin packing domains were experimentally observed upon ActD treatment compared to the control.

In addition to evaluating domain properties using ChromSTEM, we utilized live-cell partial wave spectroscopy (PWS) imaging to observe the change in chromatin organization after transcription inhibition in various cell lines (Fig. 3b). The PWS images demonstrate a significant reduction in average chromatin packing scaling upon ActD treatment in live cells across four different cell types. Next, the size of the domains is quantitatively approximated via polymer scaling relationships discussed in Supplementary Section S1.13 22,24 . The quantification of the domain sizes (boxplots in Fig. 3b) shows that, for all cell types studied, packing domains are larger for upon transcription inhibition with ActD treatmentin agreement with the ChromSTEM results on BJ fibroblasts.

Additionally, we have previously used STORM imaging to observe the nucleus wide changes in chromatin organization caused by transcription abrogation in HeLa nuclei after ActD treatment 19 . Heatmaps of chromatin density obtained via Voronoi tessellationbased color-coding of STORM images (see 19 for analysis) are shown in Fig. 3c. The zoomed in images of heatmaps of the chromatin cluster density (Fig. 3f) clearly show the increasing heterochromatin domain sizes when RNAPII activity is inhibited, in agreement with our theoretical and numerical predictions (Fig. 2d, e). Importantly, we see that the changes in chromatin organization occur not only in the interior domains of the nucleus but also along its periphery (Fig. 3f, g).

Altogether these complementary imaging techniques establish that nucleus wide increase in sizes of compacted chromatin domains occurs upon the loss of transcription in a wide range of cell lines.

The chromatin cluster density maps obtained from STORM imaging were further analyzed to quantify the sizes of heterochromatin domains after DMSO and ActD treatment. A density-based threshold was used to isolate the high-density heterochromatin regions, which were then clustered via a density based spatial clustering algorithm (see Supplementary Section S1.8) and further sub-classified into LADs and interior domains depending on the distance from nuclear

0 DMSO **** ActD (f) DMSO zoom ActD zoom D M S O A c t D Domain size (normalized) Control/DMSO ActD (e) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Heterochromatin Content (ϕ h ) **** * * * * Model D M S O A c t D DNA -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 Log 10 (Voronoi Area (nm 2 )) 1.0 1.2 1.4 1.6 1.8 1 2 3 4 5 (c) Experiment (g) (a) BJ Fibroblasts (ChromSTEM) Inner Domains LADs D M S O A c t D D M S O A c t D * * * * Inner Domains LADs Inner Domains LADs (b) DMSO ActD HeLa (STORM) Control ActD ** Domain size-scale Control ActD Domain size-scale Control ActD Control ActD Domain size-scale Domain size-scale HCT116 (PWS) A549 (PWS) U2OS (PWS) HeLa (PWS) 0.5 1.0 1.5 2.0 2.5 C o n t r o l A c t D 0.5 1.0 1.5 2.0 2.5 C o n t r o l A c t D **** 0.5 1.0 1.5 2.0 2.5 **** C o n t r o l A c t D 0.5 1.0 1.5 2.0 2.5 * C o n t r o l A c t D C o n t r o l A c t D Domain radius (nm) 0 100 200 300 ** DMSO ActD HC domain size (d) Transcription Chromatin loop extrusion Loss of extruded loops 0 DMSO **** ActD (f) DMSO zoom ActD zoom D M S O A c t D Domain size (normalized) Control/DMSO ActD (e) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Heterochromatin Content (ϕ h ) **** * * * * Model D M S O A c t D DNA -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 Log 10 (Voronoi Area (nm 2 )) 1.0 1.2 1.4 1.6 1.8 1 2 3 4 5 (c) Experiment (g) (a) BJ Fibroblasts (ChromSTEM) Inner Domains LADs D M S O A c t D D M S O A c t D * * * * Inner Domains LADs Inner Domains LADs (b) DMSO ActD HeLa (STORM) Control ActD ** Domain size-scale Control ActD Domain size-scale Control ActD Control ActD Domain size-scale Domain size-scale HCT116 (PWS) A549 (PWS) U2OS (PWS) HeLa (PWS) 0.5 1.0 1.5 2.0 2.5 C o n t r o l A c t D 0.5 1.0 1.5 2.0 2.5 C o n t r o l A c t D **** 0.5 1.0 1.5 2.0 2.5 **** C o n t r o l A c t D 0.5 1.0 1.5 2.0 2.5 * C o n t r o l A c t D C o n t r o l A c t D Domain radius (nm) 0 100 200 300 ** DMSO ActD HC domain size (d) Transcription Chromatin loop extrusion Loss of extruded loops

periphery (Supplementary Section S1.9). The quantitatively extracted distribution of interior heterochromatin domain radii for DMSO and ActD treated nuclei shows that their mean radius after transcription inhibition was nearly 1.61 times that in DMSO controls (Fig. 3g). Indeed, our model (Eq 3-4, Fig. 2d, g) predicts that loss of transcription results in increased heterochromatin domain size. This is because under control conditions, extrusion of heterochromatin phase into euchromatin occurs. We assume, based on previous experimental findings 19 , that the presence of RNAPII activity drives the supercoiling of the DNA loop, thereby extruding it from the heterochromatin phase into the euchromatin phase at the phase boundaries (Fig. 3c, left panel). However, when RNAPII is inhibited with ActD treatment (Fig. 3c, right panel), the absence of this driving force for supercoiling-mediated loop extrusion keeps more DNA in the heterochromatin phase thereby increasing the domain sizes. The in-silico chromatin distribution predicted under control (left panel) and transcription inhibited (Γ a = 0, right panel) conditions is shown in Fig. 3e. The phase separated heterochromatin domains ϕ h = ϕ max h À Á are shown in red in a loosely compacted euchromatin background (blue, ϕ h = 0). We quantify the change in the sizes of the heterochromatin domains predicted by the model as the active extrusion rate Γ a is parametrically varied. The value of Γ a under control conditions is chosen (Supplementary Table S2) such that the change in the interior domain sizes with respect to transcription inhibition (with Γ a = 0) is quantitatively the same as observed experimentally.

The model predicts changes in LAD thickness due to transcriptional inhibition with no additional parameters. Next, we quantitatively validate the choice of Γ a under control conditions by comparing the predicted change in LAD thickness against that quantified from the STORM images. Our theoretical predictions (Eq. ( 4)) show that the reduction in transcription increases the thickness of the LADs reflecting the behavior predicted in the interior of the nucleus (Fig. 2d, g). Our simulations of chromatin distribution in the nucleus (Fig. 3e) show that inhibition of transcription (Γ a = 0) results in thicker LADs. Of note, the chromatin-lamina interaction strength V L À Á stays unchanged between the two simulations. Yet, we see a higher association of chromatin with the periphery. Upon quantitative comparison (Fig. 3g, left panel) we see that the LADs grow approximately 1.37 times thicker upon loss of transcription.

To validate this prediction, we compare the predicted change in LAD thickness with that quantified from in-vivo STORM imaging. (Fig. 3g, refer to Supplementary Sections S1.8 and S1.9 for procedure). The quantified comparison of LAD thickness between DMSO and ActD nuclei (Fig. 3g) shows nearly 1.3 times increase upon ActD treatment, in close quantitative agreement with the model prediction. Overall, with both model predictions and cellular observations, our results suggest that impairment of transcription plays a significant role in determining the size scaling of the interior heterochromatin domains and LADs.

We next enquire how, in addition to altering the size of the compacted domains, abrogation of transcription changes the extent of DNA packing. For this we analyzed the chromatin distribution in HeLa nuclei under DMSO and ActD treatments from STORM images previously generated 19 . Under control conditions the distribution of DNA is qualitatively more homogenous while ActD treated nuclei exhibit more isolated distinct domains of compacted chromatin surrounded by region of very low chromatin density (Fig. 4a). For quantification, we plot the chromatin intensity along a horizontal line chosen to run across two heterochromatin domains with euchromatin between them (see zoomed images in Fig. 4b, blue and red horizontal line). The chromatin intensity, plotted in Fig. 4c (in blue) shows that even in the euchromatin region, the DNA presence is substantial. On the other hand, chromatin intensity across a horizontal line chosen across a heterochromatin domain in ActD nucleus (Fig. 4b, c; in red) shows a much steeper gradient outside the domain.

The increased presence of DNA in the euchromatic phase in presence of transcription as observed experimentally is captured by the simulations. The in-silico distribution of DNA (measured as the sum of volume fractions of the chromatin phases, ϕ e + ϕ h ) in a nuclear region far from LADs is plotted in Fig. 4d for control and transcription inhibited in-silico nuclei. We see that the euchromatic phase (outside white circles) is darker when transcription is inhibited, indicating the presence of much lesser DNA than in control euchromatin. A quantification of the total DNA along cut-lines chosen in the control and ActD in-silico nuclei confirm the observations (Fig. 4e).

Since the lack of transcription inhibits supercoiling-mediated chromatin loop extrusion from heterochromatin into euchromatin, we see a reduced density of DNA in the euchromatin phase of the nucleus under ActD conditions. Further, due to the lack of chromatin extrusion out of the heterochromatin domains when transcription is inhibited, we also observe that they are larger in size. Thus, transcription, via chromatin loop extrusion, results in removal of DNA from compacted heterochromatin region by converting it into active euchromatin form.

Taken together, our results suggest that transcription not only affects the scaling of the lengths (radius or thickness) of the heterochromatin domains, but also significantly changes the relative amounts of DNA in the euchromatin and heterochromatin phases.

We have established that change in transcription activity affects the global chromatin organization of the nucleus via altered supercoiling mediated loop extrusion. In turn, chromatin loop extrusion is initiated by the loading of cohesin onto DNA via a balance between cohesin loaders such as NIPBL and cohesin unloaders like WAPL (Fig. 1c 2,12,25,26 ). If the chromatin loop extrusion is responsible for the global chromatin reorganization, altering the cohesin loading/unloading balance must ActD treated nucleus. f Zoomed in views of DMSO and ActD treated nuclei localized to the nucleus interior (top panels) and the periphery (bottom panels). Red and blue boxes shown in c are zoomed into. All scale bars -1 μm. g Left: Simulations predict domains in ActD nuclei are on average 1.63 times larger than in DMSO nuclei (n = 127 (DMSO), 77 (ActD) unpaired two tail t-test, p = 0) while LADs are 1.37 times thicker (n = 38 (DMSO), 15 (ActD); unpaired two tail t-test, p = 0). Right: Domain radii observed experimentally in ActD treated nuclei (n = 3584 loci, 20 nuclei) are 1.61 times (unpaired two tail t-test, p = 0) larger than in DMSO nuclei (n = 5830 loci,19 nuclei), while LADs are 1.3 times thicker (n = 1082 loci (DMSO), 1015 loci (ActD), unpaired two tail t-test, p = 0.0006). All boxplots show the mean (cross), median (horizontal line), upper and bottom quartiles (box outlines) and the maximum and minimum non-outlier data points (whiskers). All source data are provided as a source data file. also result in chromatin reorganization. Thus, next, we study the chromatin arrangement in WAPL-deficient (WAPLΔ) nuclei marked by increased levels of loaded cohesin.

In vivo, WAPL depletion causes an accumulation of large amounts of cohesin on chromatin 27 . This results in a much more homogenous distribution of DNA, which was previously termed "blending" due to excessive extrusion of chromatin loops, as shown schematically in Fig. 5a 19 . In our mathematical model, WAPL deficiency is simulated as an increase in the rate of chromatin extrusion ðΓ a Þ. Based on the theoretical size scaling of the interior heterochromatin domains and LADs, as seen from Eq. ( 3) and Fig. 2g, our model predicts that increase in Γ a would result in a decrease in the radius of the steady state heterochromatin domains (Fig. 5b).

STORM images of HeLa nuclei without (labeled Cas9) and with WAPL-deficiency previously revealed genome-wide changes in the chromatin organization induced by excessive loading of cohesin (Fig. 5c, d) 19 . A visual comparison between representative zoomed-in regions (white boxes in Fig. 5c) demonstrates the reduction of heterochromatin domain sizes in the interior of the nuclei in WAPLΔ nuclei (Fig. 5d). Using clustering analysis (refer Supplementary Section S1.8 and S1.9), we quantify the altered chromatin domain sizes in control and WAPLΔ HeLa cell nuclei. We observe that WAPLΔ nuclei with increased chromatin blending have heterochromatin domains with a mean radius approximately 15% smaller than control nuclei (Fig. 5e).

In-silico, we parametrically vary the active chromatin extrusion rate Γ a above the control level (Supplementary Table S2, determined for control treatment). The value of Γ a for WAPLΔ nuclei is chosen (Supplementary Table S2) such that the decrease in the size of interior heterochromatin domains reduces by 15% (Fig. 5f) to agree with the experimental observation (Fig. 5e).

As discussed previously (Fig. 2g), the model predicts that the effects of chromatin extrusion observed in the interior domains of the

Total DNA (ϕ h +ϕ e ) 8 . 0 0 . 1 6 . 0 4 . 0 2 . 0 0 . 0 Normalized distance DMSO ActD (d) (e) DMSO ActD 0 0 . 2 0 . 4 0.8 DMSO ActD x -distance (μm) Chromatin Grayscale Intensity DMSO Zoom ActD Zoom DMSO ActD 3 μm 3 μm (b) (a) - ----------------------Experimental Analysis ------------------------- nucleus are replicated along the nuclear periphery. Simulation of nuclear chromatin organization (Fig. 5b) reveals that by changing only the rate of chromatin extrusion Γ a , keeping all other parameters including chromatin-lamina interaction potential V L constant, we see a reduction in the association of chromatin with the lamina. Specifically, a 2.5-fold increase in Γ a calibrated to occur due to WAPL-deficiency predicts a 51.2% decrease in the average LAD thickness, as shown in Fig. 5f.

The predicted change in LAD thickness is consistent with previous experimental observations and was further quantitatively validated by measuring the thickness of LADs in STORM images of control and WAPLΔ nuclei (Fig. 5e) 19 . A reduction in the sizes of domains, as seen in the nucleus interior, can also be observed at the nuclear periphery, as shown in a representative zoomed in region (white dashed boxes in Fig. 5c) in Fig. 5d. The mean thickness of the LADs at the nuclear periphery is approximately 20% smaller for WAPLΔ nuclei (Fig. 5h) as compared to the control-treated nuclei.

Together, these results confirm that the meso-scale spatial chromatin organization is strongly regulated by the chromatin loop (n = 2386 loci in 6 nuclei for Cas9-treatment and 2416 loci in 7 nuclei for WAPLΔ treatment). WAPLΔ treated nuclei exhibit a significantly lower ( ∼ 0.86 times) mean heterochromatin radius (unpaired two-tailed t-test, p = 6e-10). Quantification of LAD thickness along the periphery of Cas9-and WAPLΔ -treated nuclei. (n = 219 loci in 6 nuclei for Cas9-treatment and 169 loci in 7 nuclei for WAPLΔ treatment). WAPLΔ treated nuclei exhibit a significantly lower ( ∼ 0.43 times) mean LAD thickness (unpaired two-tailed t-test, p = 1e-13). f Boxplot in left panel shows the distribution of domain radii predicted numerically. WAPLΔ nuclei have a mean domain radius 0.8 times that of Cas9-treated nuclei (unpaired two-tailed t-test, p = 0). Boxplot in right panels shows the distribution of LAD thicknesses predicted numerically. WAPLΔ nuclei have a mean LAD thickness 0.82 times that of Cas9treated nuclei. All boxplots show the mean (cross), median (horizontal line), upper and bottom quartiles (box outlines) and the maximum and minimum non-outlier data points (whiskers) of the plotted distribution. All source data are provided as a source data file.

formation, and this effect can be modulated not only by the transcription activity, but also by altering the extent of loading or unloading of cohesin rings on the DNA. These results provide further evidence for the link between transcriptional regulation and nucleuswide chromatin distribution via transcription-driven supercoiling mediated chromatin loop extrusion.

Since we have established, via both quantitative analysis of experimental data and simulations, that extrusion of chromatin loops is governed by both cohesin loading/unloading balance and RNAPII mediated transcription, a question of their tandem role emerges.

To simulate the individual effects of cohesin loading and transcriptional activity, we decompose the overall active chromatin extrusion rate into its distinct constitutive steps. The individual steps involved in the process of supercoiling mediated chromatin loop extrusion from heterochromatin into euchromatin (as discussed previously in Section "Introduction") are shown in Fig. 6a. As a first step, a balance between the loading of cohesin via NIPBL/MAU2 25 on chromatin occurring at a rate Γ l and its unloading via by WAPL/PDS5 2,12,26 occurring at a rate Γ ul results in the association of cohesin rings with chromatin at an overall rate Γ coh = Γ l À Γ ul . In other words, Γ coh denotes the overall rate of cohesin loading on DNA. The entrapment of DNA by cohesin is followed by the extrusion of supercoiled loops of chromatin via DNA supercoiling by the RNAPII mediated transcription, at a rate denoted by Γ tr . Thus, as shown in Fig. 6a, by assuming a first-order As previously, ActD treated nuclei exhibited a significantly increased domain size (unpaired two-tailed t-test, p = 0) while WAPLΔ treated nuclei exhibit a significantly lower mean heterochromatin radius (unpaired tw-tailed t-test, p = 0). However, the differences between Cas9-ActD treated and WAPLΔ-ActD treated nuclei was insignificant (unpaired two-tailed t-test, p ∼ 0.9). All boxplots show the mean (cross), median (horizontal line), upper and bottom quartiles (box outlines) and the maximum and minimum non-outlier data points (whiskers) of the plotted distribution. All source data are provided as a source data file.

reaction kinetics for both steps, the overall rate of active chromatin extrusion Γ a at the interface of heterochromatin and euchromatin is proposed to be multiplicatively decomposed as,

In addition to the extrusion of loops via RNAPII mediated DNA supercoiling activity 12,13,19,[28][29][30] , in vitro experiments proposed that cohesin once transiently loaded onto DNA, could independently drive the formation of loops via its ATPase machinery 9,11,[31][32][33] . Cell based experiments demonstrated that in WAPLΔ cells, clusters of cohesin in WAPLΔ cells assemble together into vermicelli-like structures and these structures disappear upon transcription inhibition, but not upon partial loss of cohesin 19 . These results, taken together, present strong evidence for the important role of transcription in powering cohesin mediated loop extrusion. While the relative role of cohesin's motor activity and transcription in loop extrusion inside cells remains to be determined, here we focus on the latter given the previous in vivo experimental findings. We indeed show that a kinetic model captured by Eq. ( 5) sufficiently explains the effect of extrusion of the specific chromatin loops extending from transcriptionally silenced heterochromatin into genetically active euchromatin on determining the meso-scale chromatin domain sizes.

The chromatin organization is simulated in a nucleus under control and transcription inhibition treatments for nuclei with and without WAPL deficiency. The chromatin organization in a control nucleus (labeled Cas9-DMSO), simulated via parameters listed in Supplementary Table S1 is shown in Fig. 6b, top-left panel. The individual inhibition of transcriptional activity without affecting the cohesin loading (Cas9-ActD) results in a chromatin organization with increased heterochromatin domains sizes and LAD thickness, as shown in Fig. 6b, top-right panel. On the other hand, the simulation of chromatin distribution in nucleus with depleted cohesin unloading, without disturbing the transcriptional activity, (WAPLΔ-DMSO) is shown in Fig. 6b, bottom-left panel. Finally, the chromatin distribution predicted in a WAPLΔ nucleus with inhibited transcription (WAPLΔ-DMSO-treatment) is shown in Fig. 6b, bottom-right panel. As shown in Fig. 3e and Fig. 3g, ActD (mathematically, Γ tr = 0 in Eq. ( 5)) results in larger heterochromatin domains and thicker LADs, while WAPLΔ nuclei (increased cohesin loading; mathematically, Γ ul =Γ l increases in Eq. ( 5)) show the opposite effect with smaller heterochromatin domains and LADs. For a WAPLΔ nuclei in which transcription is inhibited (WAPLΔ -ActD; mathematically, Γ tr = 0 and Γ ul =Γ l increases in Eq. ( 3)), the model predicts that inhibition of transcription returns the chromatin organization to the control (Cas9-ActD) levels. Transcription inhibition thus blocks the reduction in chromatin domain sizes induced due to WAPL deficiency due to lack of impetus for chromatin supercoiling.

To quantitatively validate the model predictions, we investigate the in-vivo chromatin organization under individual and tandem changes in transcription and cohesin unloading by re-analyzing previously reported super-resolution images shown as heatmap density plots in Fig. 6c 19 . Visual inspection of this data agrees with the model predictions that transcriptional inhibition counteracts the chromatin blending observed in DMSO treated WAPLΔ nuclei, which was also previously reported 19 . We thus focused on extracting the radius of heterochromatin domains and LAD thickness to further validate the model results quantitatively (Fig. 6d). Cas9 -ActD treated nuclei show an increased heterochromatin domain radius compared to control while WAPLΔ nuclei show a significant reduction in domain radius and LAD thickness (Fig. 6d). However, WAPLΔ -ActD treated nuclei show no significant difference in comparison to Cas9 -ActD treated nuclei (Fig. 6d), in quantitative agreement with the numerical predictions.

These results further confirm that the effect of transcription on global chromatin distribution occurs via supercoiling mediated chromatin loop extrusion, especially at the interface of heterochromatin and euchromatin phases. Furthermore, these results also present a significant validation of the mathematical phase-field model of chromatin organization in the nucleus.

Significant inroads into mechanistic modeling of chromatin organization as physically and functionally distinct states with finer architectural sub-features such as topologically associated domains (TADs) and chromatin loops have been made from a polymer physics perspective. Such models were developed with different levels of finegraining to capture biophysics of chromatin organization at different length-scales spanning single or multiple nucleosomes [34][35][36] , multiple nucleosome clutches [37][38][39][40][41][42] , single and multiple TADs with salient sub-TAD features [43][44][45][46] , single and multiple chromosomes 45,[47][48][49][50][51][52] and the whole genome [53][54][55][56] . Depending on the focus on the chromatin functionalities or structure being simulated, any of these models can be adopted. For instance, first-principles thermodynamics driven approach may capture chromatin as a copolymer with two states, whereas a data-driven approach trained on conformation capture (e.g., Hi-C) or sequencing (e.g., CHIP-seq) data may incorporate over 50 states spanning the entire genome.

The experimentally observed role of RNAPII-mediated transcription in DNA supercoiling and subsequent loop extrusion 2,12,19,20,33,[57][58][59][60] has also been studied using molecular dynamics simulations and polymer physics-based models at nanoscale 20,[61][62][63][64][65][66] . However, quantitative predictions of sizes of heterochromatin domains which organize at a nucleus-wide meso-scale level are beyond the purview of such models. Furthermore, to the best of our knowledge, polymer models lack the far from equilibrium kinetic considerations of active epigenetic regulation, chromatin extrusion and diffusion kinetics, which we find are intricately involved in the spatiotemporal regulation of heterochromatin domain sizes. The current study, incorporating coarsegrained continuum model of chromatin organization at a mesoscale, presents the following advantages over previous polymer-based models:

a. Nucleus-wide characteristic size distribution of heterochromatic domains. b. LADs of finite thickness co-existing with interior heterochromatin domains and their dynamic size-regulation. c. The kinetic interplay of diffusion, epigenetic reactions and transcription in regulation of meso-scale organization.

Thus, here we present a non-equilibrium thermodynamic continuum model of the meso-scale chromatin organization in the nucleus to bridge the gap in the understanding of the mechanistic relation between transcriptional and epigenetic regulation and the size-scaling of the meso-scale heterochromatin domains. Our model incorporates the energetics of chromatin-chromatin interactions which is constructed as a double-well function allowing the phase-separation of chromatin into compartments of distinct compactions. Along the nuclear periphery, the effect of chromatin-anchoring proteins such as LAP2β is captured via energetic chromatin-lamina interactions leading to the formation of LADs. Concomitant with the energetics, the chromatin organization is temporally driven by diffusion kinetics of nucleoplasm and the effective diffusion-like evolution of epigenetic marks. While the diffusion of nucleoplasm determines the level of chromatin compaction, such that higher local nucleoplasm content results in lesser chromatin compaction, diffusion of epigenetic marks results in accumulation of acetylated and methylated nucleosomes driving their segregation (Supplementary Fig. S4, Supplementary Sections S1.6, S1.7 of the Supplementary Information). Most importantly, we also account for the active reaction kinetics, which allow the interconversion of heterochromatin into euchromatin and vice-versa. The chromatin phase-interconversion can occur via the epigenetic regulation of chromatin in the nucleus via the acetylation or methylation of the histones. Finally, to capture the role of transcription mediated supercoiling-driven loop extrusion in the determination of the heterochromatin domain sizes, we incorporate a kinetic conversion of compacted chromatin into transcriptionally active euchromatin in presence of RNAPII (Fig. 7a, b).

Together the active transcriptional kinetics and epigenetic regulation determine the active interconversion of hetero-and euchromatin, thereby taking the chromatin organization in the nucleus to a dynamic steady-state configuration. Specifically, our theoretical analysis reveals that the active reaction kinetics aloneindependently of energetic interactionsoffers a significant control over the average extent of chromatin compaction in the nucleus, thereby breaking the detail balance of thermodynamic equilibrium. At the meso-scale, spanning individual heterochromatin domains, we theoretically observe that the distribution of epigenetic marks relevant to chromatin compaction exhibit a radial gradient which would drive an inward heterochromatin flux leading to ripening of the phase-separated domains. However, the presence of epigenetic and transcriptional regulation offers an opposition to the influx via -(a) acetylation of heterochromatin into euchromatin which is then pushed out via the diffusion of epigenetic marks, and (b) extrusion of loops of chromatin from the heterochromatin phase into euchromatin phase (refer Fig. 2e-g). The steady-state balance between the opposing fluxes leads to intrinsic emergence of a characteristic size-scaling of heterochromatin domains. Upon translating the theoretically and numerically obtained stable heterochromatin domain size-scale into physical dimensions (Supplementary Section S9, SI), we note that the predicted characteristic domain size is equivalent to that observed using multiple super-resolution imaging techniques 5,7,[67][68][69][70] . It is essential to note that without the active chromatin phase interconversionat thermodynamic equilibriumno inherent size-scale of heterochromatin domains would be observed.

Thus, our model predicts that transcriptional activity, synergistically with epigenetic regulation, controls the size and morphologies of

(a) Euchromatin Heterochromatin Cohesin loading Cohesin unloading Loop extrusion Transcription induced supercoiling STEP 1 STEP 2 Chromatin Loop Extrusion Euchromatin Heterochromatin Cohesin -loaded RNAPII Cohesin -unloaded NIPBL/MAU2 WAPL/PDS5 CTCF D N A lo o p s in e u c h r o m a t in DNA loops in heterochromatin (b) Control (c) HC domain size WAPLΔ ActD WAPLΔ + ActD Excessive cohesin loading HC domain size Transcription Loaded coehsin forms vermicelli structures. Domains become smaller Loss of chromatin loops, increased domain size. the heterochromatin domains. The key predictions of the model as summarized in Fig. 7c are: 1. Upon transcriptional inhibition, the characteristic sizes of heterochromatin domains increase due to loss of supercoiling mediated DNA loop extrusion (Fig. 7c, left panels). 2. The increased size of heterochromatin domains upon transcription abrogation are also observed in the vicinity of the nuclear lamina. 3. Transcriptional inhibition leads to reduction of DNA in the euchromatic phase. 4. Conversely, upon increased loop extrusion due to excessive cohesin loading, the size of the heterochromatin domains reduces in the interior as well as periphery of the nucleus (Fig. 7c, top panels). 5. Transcriptional inhibition in nuclei with excessive cohesin loaded, results in loss of loop extrusion resulting in increased domain sizes (Fig. 7c, bottom right panel).

Being founded on fundamental non-equilibrium thermodynamic principles, the predictions made by our model are cell-type agnostic. To validate cell-type independence, complementary techniques such as Chrom-STEM (for high-resolution chromatin conformation imaging) and PWS (for high-throughput nano-scale sensitive live-cell imaging) are carried out for BJ fibroblast cells and multiple epithelial cancer cell-lines -U2OS, HeLa, A549 and HCT116. We found that in vivo alterations in chromatin organization under transcriptional inhibition conditions are consistent with our model's predictions across all studied cell lines. A quantitative analysis of previously reported 19 superresolution STORM images of nuclei further gives a direct quantitative validation of the predicted effects of transcription abrogation on heterochromatin domain sizes. Of note, in all the reported cells, growth of condensed heterochromatin domains after ActD treatment is seen throughout the nucleus, including along the nuclear periphery where an increased LAD thickness is observed both in-silico and in cells. Lastly, our predictions on the changes in heterochromatin domain sizes upon over-extrusion of chromatin loops with and without transcription are quantitatively validated by domain size analysis of the previously reported 19 super-resolution STORM images of control and WAPLΔ nuclei after DMSO or ActD treatments.

In addition to imaging techniques using multiple modalities, previously reported observations of chromatin reorganization via chromatin conformation capture studies further confirm our model predictions 71,72 . While the loss of RNAPII only had a limited effect on the presence of chromatin loops, observed as off diagonal peaks on Hi-C contact maps [73][74][75] , investigations of contact maps at a finer resolution using Micro-C 72,76,77 revealed the existence of RNAPII-associated chromatin loops which are indeed disrupted upon transcription inhibition. Specifically, in agreement with our predictions, it was reported that depletion of RNAPII decreased genome-wide histone acetylation (specifically H3K27ac) levels, reduced local chromatin accessibility, and lead to loss of chromatin loops upon RNAPII depletion 72 . Further, Hi-C contact maps revealed an increase in chromatin loops in WAPL deficient nuclei 71 . However, it should be noted that not all the loops observed using Hi-C or Micro-C were RNAPII associated.

Beyond transcription induced supercoiling-driven loop extrusion, cohesin itself can play an active role in the formation, extrusion and maintenance of chromatin loops at different physiological length scales. Loops of chromatin, identified as peaks in chromatin contact mapping techniques like Hi-C and Micro-C [9][10][11]72 are formed at a length scale of topologically associated domains (TADs) or below. Such loops are observed within both active A 72 and inactive B 78 compartments. The mechanism of extrusion of such loops, could be distinct from RNAPII transcription induced supercoiling-driven loop extrusion, such as cohesin subunit SMC motor activity 20,79,80 or a passive cohesin diffusion along the chromatin polymer 81,82 . These multiple mechanisms of loop formation could be convergently cooperating forming chromatin loops at multiple physiological length scales.

However, here we are specifically interested in the meso-scale roles played by supercoiling-driven extrusion of chromatin loops from the silenced heterochromatin phase into transcriptionally active euchromatin region (Fig. 7a). These loops are specifically considered since the interconversion of heterochromatin to euchromatin will alter the sizes of the heterochromatin domains. Recent experimental evidence, as well as computational models, present strong evidence in favor of DNA loop extrusion mediated by the RNAPII driven transcription induced DNA supercoiling 12,14,19,20,57,58,79,83 . Negatively supercoiled DNA regions are particularly rich in transcription start sites (TSS) with a strong correlation seen between transcription and supercoiling 30 . Indeed, super-resolution images show high presence of RNAPII at the heterochromatin-euchromatin phase boundaries where loops would extrude from heterochromatin into euchromatin phase 22 .

Intriguingly, previous observations 19 also show that in HeLa nuclei WAPL deficiency introduces abnormalities in the peripheral distribution of lamin A/C. Since lamin A/C plays an integral role in the chromatin-lamina interactions via chromatin anchoring proteins such as LAP2β and emerin, it can be conjectured that WAPL treatment may affect the LAD organization. In the current study we have ignored such effects focusing purely on the role of supercoiling mediated chromatin loop extrusion. Our model can be easily modified to address the LAD alterations by introducing WAPL deficiency dependent modulations in the chromatin lamina interaction parameter V L in Eq. ( 2). Experiment guided modifications in the model will further strengthen our predictions of LAD formation. Further, the transcriptional machinery involves a highly complicated multi-stage process comprising recruitment of multiple transcription factors, RNAPII and gene regulatory elements, we have assumed the cohesin loading and RNAPII mediated supercoiling to be the rate defining steps which thereby govern the timescale for chromatin loop extrusion. A more refined kinetic model of transcription and loop extrusion could possibly be incorporated to predict the spatiotemporal chromatin arrangement in the nucleus. However, even without these inclusions, we believe that our model lays a fundamental computational framework to better understand the mechanistic role of transcription, and in general chemo-mechanical cell-signaling, on the meso-scale chromatin organization.

At the meso-scale, chromatin is organized into distinct transcriptionally dissimilar phases of euchromatin and heterochromatin as depicted schematically in Fig. 1a. We incorporate the energetic interactions between the nucleosomes depending on their epigenetic state, as discussed below. While on one hand entropic contributions push chromatin towards a homogenous organization, enthalpy arising from nucleosome-bridging via HP1 proteins 84 , via ionic interactions within chromatin phases 85,86 or local activity 87 oppose it. The entropicenthalpic competition comprising the chromatin-chromatin interactions drives the phase separation of chromatin domains. The emergent formation of domains occurring thermodynamically in our model is similar to the chromatin domains qualitatively postulated based on super-resolution imaging 5,68,69 . Near the nuclear periphery, heterochromatin can be further anchored to the nuclear lamina via anchoring proteins such as LAP2β [88][89][90] : This drives the formation of peripheral LADs. The interior and peripheral heterochromatin domain formation occurs spatiotemporally via free energy lowering diffusion of nucleoplasm, and diffusion-like evolution of acetylation or methylation marks on the histones. We incorporate the active interconversion between the eu-and heterochromatic phases in two ways: histone methylation or acetylation reactions can change the epigenetic distribution or transcription mediated supercoiling-driven chromatin loop extrusion.

As shown in Fig. 1a, supercoiling-driven DNA loop extrusion, can occur broadly in two regions where RNAPII is present 23,68,69,91 : within the euchromatin domains (red dashed circle in Fig. 1a) or at the interface of heterochromatin and euchromatin phases (black circle in Fig. 1a). Since the chromatin extrusion in the euchromatin phase maintains its transcriptionally active status and does not lead to any significant mesoscale changes in the epigenetic distribution, we focus on the domain interface. The chromatin extrusion at the interface is instrumental in the regulation of size of heterochromatic domains at the periphery to form euchromatin.

At any point x in the nucleus, at a time t, we consider three nuclear constituents, namely the nucleoplasm and the two phases of chromatin, euchromatin and heterochromatin with their volume fractions (refer Supplementary Information SI, Supplementary Section S1.1 for detailed definition) ϕ n x,t ð Þ, ϕ e x,t ð Þand ϕ h x,t ð Þ. We assume that these three constituents are space filling, and their volume fractions add up to unity, i.e., ϕ e + ϕ h + ϕ n = 1 (derived in SI, Supplementary Section S1.1). Hence, if the volume fractions of two of the constituents is known, the volume fraction of the third is determined by this constraint. The composition of the constituents can thus be defined in terms of two independent variables (refer to the methods for details) -(i) ϕ n x,t ð Þ volume fraction of the nucleoplasm, and (ii)

Þ which is the difference of the volume fractions of heterochromatin and euchromatin. Note that ϕ d < > ð Þ0 for the euchromatin (heterochromatin) rich phase and is therefore analogous to an order parameter. In terms of the chromatin composition variables, ϕ n and ϕ d , the free energy density at any point x can be expressed non-dimensionally as (refer Supplementary Section S1.5 for details on non-dimensionalization),

The construction of the free energy density function is discussed in more detail in the Supplementary Section S1.2. The first term, which is a Flory-Huggins type free energy density for chromatin, defines the competition between the enthalpy of the chromatin-chromatin interactions and entropic contributions of chromatin configuration. We discuss the choice of the form of chromatin-chromatin energetic interactions, and its similarity to the Flory-Huggins form of free energy density in the Supplementary Section S1.3. This term gives rise to the double-well potential describing the energy landscape of the possible chromatin distribution. The potential surface is visualized in Fig. 1b as a contour plot with well locations as ϕ h = 0 (euchromatin phase) and ϕ h = ϕ max h (heterochromatin phase). The well towards the bottom in Fig. 1b corresponds to the heterochromatin phase with a low water content and a higher chromatin compaction.

The methylated histone tails in heterochromatin phase can mediate inter-chromatin interactions via chromatin cross-linkers such as HP1α [92][93][94] : Such chromatin crosslinking lowers the enthalpy resulting in a heterochromatin phase well with a densely packed chromatin. On the other hand, the euchromatin well, corresponding to the energy minimum with a higher water content is marked with a more acetylated histone tails with a loosely packed chromatin conformation corresponding to a higher entropy.

The second term captures the interactions between the chromatin and the lamina via chromatin anchoring proteins (LAP2β, emerin, MAN1, etc.) [88][89][90] with parameter e V L denoting the rescaled strength of these anchoring interactions. Notably, these interactions are most robust at the nuclear periphery (distance from lamina d = 0) and vanish exponentially over a length scale d 0 . Since the chromatin domains preferentially associating with the nuclear lamina are linked to transcriptional repression and an increased histone methylation 89,[95][96][97] , the chromatin-lamina interactions are captured specifically towards heterochromatin phase. Lastly, the negative sign permits an energetic preference for the peripheral association of heterochromatin. Analogous discrete descriptions of chromatin-lamina interactions, via formation of strong bonds when the chromatin is within a characteristic distance from the lamina have been previously implemented 50,51 in polymer models of chromatin, although without the epigenetic or transcriptional kinetics.

The last term accounts for the interfacial energy which is not accounted in a Flory-Huggins model and penalizes the formation of sharp interfaces between the dissimilar phases (refer Supplementary Sections S1.2 and S1.5). The interfacial penalty competes with the energy of chromatin-chromatin interactions forming smooth interfaces of non-dimensional width δ (Supplementary Section S1.5).

Thus, the energetic considerations dictate that an initial chromatin configuration (light blue circle in Fig. 1b) spontaneously phaseseparates into the two energy wells to minimize the total free energy of the system. The driving force pushing the chromatin composition towards the energy wells is a measure of the gradients of the energy landscape and is called the chemical potential. Thus, the chemical potentials are obtained at each point in space by considering changes in energy density for small changes in the local volume fractions (labeled n or d): e μ nðdÞ x,t ð Þ= δ e W δϕ nðdÞ , as derived in Supplementary Eq. (S12). Here, the operator δ denotes the functional derivative, or the change in free energy density with respect to the volume fraction. Spatial gradients of chemical potential drive the diffusive flow of nucleoplasm to reduce the overall free energy of the system giving rise to nucleoplasm kinetics via Supplementary Eq. (S6) (Supplementary Section S1.4). By rescaling the evolution equation, via the methodology described in Supplementary Section S1.5, we obtain the nondimensional nucleoplasm kinetics as (Supplementary Eq. (S13),

Note that nucleoplasm diffusion kinetics in Eq. ( 7a) is conservative in nature, i.e., the net amount of water in the nucleus is conserved over time as long as no water enters or exits the nucleus.

The kinetics of epigenetic marks on the histonesacetylation or methylationcan have two contributions. Primarily, the epigenetic regulation via enzymes such as histone deacetylase (HDAC), histone methyltransferase (HMT), histone acetyltransferase (HAT) and histone demethylase (HDM) can result in interconversion of heterochromatin and euchromatin phases via acetylation and methylation reactions as shown in Supplementary Fig. S3. The reaction kinetics, inherently nonconservative (discussed in Supplementary Section S1.4) are captured via the second term in Eq. (7b').

The reaction kinetics however should also incorporate the contribution of chromatin-chromatin interactions, which determine how favorable the euchromatin-heterochromatin interconversion is depending on the epigenetic marks on the other nucleosomes in vicinity. In the SI, we qualitatively (Supplementary Section S1.6) and theoretically (Supplementary Section S1.7) describe how neighborhood dependent reaction-kinetics is effectively equivalent to diffusion-like evolution of epigenetic marks, which we incorporate in our model as 'diffusion of epigenetic marks' (first term in Eq. 7b'). Thus, the reaction kinetics also give rise to an effectively conservative contribution which allows for evolution of epigenetic marks without changing the overall amounts of heterochromatin and euchromatin in the nucleus.

Lastly, the kinetics of transcription-mediated supercoiling-driven chromatin extrusion localized at the heterochromatin domain boundaries is incorporated into the dynamics of epigenetic marks, giving rise to a second evolution equation of the form,

Note that Eq. (7b') is non-dimensionalized by rescaling all time variables with respect to the rate of histone tail acetylation (i.e., e t = Γ ac t) and all spatial variables with respect to the characteristic reaction diffusion length defined in Supplementary Section S1.5 (i.e., e x = x=l RD ). The second term in Eq. (7b') incorporates active first-order reaction kinetics of histone tail acetylation (see Eqs. (S7) and (S8) before rescaling) and that of histone methylation e Γ me leading to interconversion between hetero-and eu-chromatin.

The last term in Eq. (7b') accounts for the supercoiling-driven chromatin extrusion kinetics and the chromatin state changes resulting from it (Fig. 1c). Being transcription mediated, the kinetic rate of supercoiling-driven extrusion e Γ a e x À Á must be spatially dependent on local availability of RNAPII, which is prominently present at the boundaries of the compacted heterochromatin phase 23,68,69 . Although supercoiling-driven loop extrusion may also occur within the euchromatin phase, it does not contribute to interconversion of chromatin phases as euchromatin is already transcriptionally active. In contrast, at the interface of heterochromatin and euchromatin, supercoilingdriven loop extrusion can result in activation of otherwise inactive genes. Considering this spatial localization to the heterochromatin domain boundaries, we rewrite Eq. (7b') as,

Note that ϕ max h =2 is the volume fraction of heterochromatin at the domain boundary. A deviation of Δϕ from this value defines the width of the domain boundary, and the supercoiling-driven loop extrusion is spatially restricted to a narrow region at the boundary of heterochromatin domains. The last two terms of Eq. (7b) are responsible for the non-conservative dynamics and can alter the global heterochromatin to euchromatin ratio of the system. More detailed derivation of the chemical potential, contribution of passive diffusion kinetics, epigenetic and active loop extrusion can be found in the extended methods section in the SI (Supplementary Sections S1.1-S1.7).

Having developed the model to capture the spatiotemporal organization of chromatin in the nucleus, we numerically solve Eqs. (7a) and (7b) along with the equation defining the chemical potential (Supplementary Eq. ( S3)). As a boundary condition we ensure no exchange of water and chromatin between the nucleus and the surroundings. This condition can be suitably adjusted to allow flow of water from or into the nucleus. The parameters used in the model along with the initial and boundary conditions are described in detail in the SI (Supplementary Section S8) and listed in Supplementary Table S2. Note that the epigenetic rates e Γ me and the strength of chromatin-lamina affinity e V L are not modified throughout any of the simulations carried out, unless explicitly stated. This is to ensure that any predicted changes in chromatin organization occur specifically due to changes in supercoiling-driven chromatin loop extrusion.

To investigate the organization of chromatin in the nucleus, we develop a mathematical model for the phase separation of heterochromatin and euchromatin considering chromatin-chromatin interactions, chromatin-lamina interactions, epigenetic regulation of chromatin via histone acetylation or methylation, and the role of transcriptional regulators. We consider three nuclear constituents -nucleoplasm and chromatin in either heterochromatin or euchromatin phases.

At any point 𝒙 in the nucleus, at a time 𝑡, consider an infinitesimal observation window of volume 𝑉(𝒙, 𝑡). Let the volume of nucleoplasm, euchromatin and heterochromatin within this observation window be 𝑉 𝑛 (𝒙, 𝑡), 𝑉 𝑒 (𝒙, 𝑡) and 𝑉 ℎ (𝒙, 𝑡). Note that 𝑉 = 𝑉 𝑛 + 𝑉 𝑒 + 𝑉 ℎ . Thus,

The ratio of volume of each component to the total volume of the infinitesimal observation window is defined at the volume fraction of the component 𝜙 𝑖 (𝒙, 𝑡), for 𝑖 = 𝑛, 𝑒, ℎ, and determines the content of nuclear constituents -nucleoplasm, euchromatin and heterochromatin at each point 𝒙 and time 𝑡 in the nucleus. Thus,

Thus, the physical state of the nucleus at any point can be defined by volume fractions of any two nuclear constituents, with the third constrained via the above equation. Equivalently, the physical state of the nucleus at any point can be completely defined via two independent variables -(i) volume fraction of nucleoplasm 𝜙 𝑛 (𝒙, 𝑡), and (ii) difference between the volume fractions of heterochromatin and euchromatin ϕ 𝑑 (𝒙, 𝑡) = 𝜙 ℎ (𝒙, 𝑡) -ϕ 𝑒 (𝒙, 𝑡). 𝜙 𝑑 can be considered an order parameter which when negative implies a euchromatin rich phase and when positive implies a more condensed heterochromatin rich phase. While the change of the variables is entirely equivalent mathematically, physiologically such a description permits a natural definition of the movement of the two mobile species in the nucleus -nucleoplasm or water, and the epigenetic marks of acetylation or methylation.

In terms of the independent variables 𝜙 𝑑 (𝒙, 𝑡) and 𝜙 𝑛 (𝒙, 𝑡), the free energy density at any point 𝒙 can be expressed as 𝑊(𝜙 𝑛 , 𝜙 𝑑 , ∇𝜙 𝑛 , ∇𝜙 𝑑 ) where we have also incorporated the energetic considerations associated with the phase interfaces via the spatial gradients of the volume fractions. Specific form of the free energy can be invoked by considering the various energetic contributions in the nucleus such as,

• The first term in Eq S1 arises from the competition between entropy and enthalpy of mixing heterochromatin and euchromatin phases. It is equivalent to the Flory-Huggins free energy description, as discussed in subsection S1.3. This term gives rise to the double-well form of the free energy landscape, as shown in the contour plot in Figure 1b. The two wells, shown as red and blue dots, are the energy minima corresponding to the two stable phases of chromatin -a water-rich, loosely packed euchromatin phase (𝜙 ℎ = 0) or a compacted waterdevoid heterochromatin phase (𝜙 ℎ = ϕ h max , 𝜙 𝑛 ~0). Here ϕ h max denotes the extent of compaction in the heterochromatin phase. Any initial chromatin configuration will spontaneously phase separate into heterochromatin and euchromatin domains (red arrows). • The second term captures the interactions between the chromatin and the lamina via chromatin anchoring proteins (HDAC3, LAP2𝛽, emerin, etc. [1-3]) with parameter 𝑉 𝐿 denoting the strength of these anchoring interactions. Since these interactions are mediated by proteins anchored on the lamina, we take an interaction strength that is most robust at the nuclear periphery (distance from lamina 𝑑 = 0) and decays exponentially away from the lamina over a length scale 𝑑 0 (schematically shown in Figure S1). Note that the exact decay characteristics do not really affect the model qualitatively.

We choose exponential form as a generalized minimal-parameter decay. The negative sign permits an energetic attraction for the heterochromatin phase along the nuclear periphery. We chose the heterochromatin phase specifically to interact with the lamina since the chromatin domains preferentially associating with the nuclear lamina are linked to transcriptional repression and an increased histone methylation [2, 4-6].

• The last term in Eq S1 denotes the energy penalty associated with forming phase boundaries between the euchromatin and heterochromatin phases as they separate. The term 𝜂 𝑛,𝑑 is the increase in the energy due to formation of a unit width of the interface. Note that the term |∇𝜙 𝑛,𝑑 | is the magnitude of the slope of the interface. Notably, as 𝜂 𝑛,𝑑 increases, there is a greater penalty on formation of sharp interfaces, resulting in more smooth interfaces which are wider. Thus 𝜂 𝑛 and 𝜂 𝑑 directly control the width and the energy of the phase boundaries.

In our simulations, we choose 𝜂 𝑛 = 𝜂 𝑑 = 𝜂.

The total free energy of the nucleus of volume Ω can be written, after incorporating the work done in exchanging water between the nucleus and the cytoplasm, as,

Here, 𝜕Ω 𝑝 denotes the surface of the nucleus where pores and channels allow exchange of water between nucleus and cytoplasm by maintaining an 'external' chemical potential 𝜇̅ 𝑛 . 𝐼 𝑛 is the volume of water entering per unit surface area into the nucleus, as a result. The double-well energy landscape described by Eq 1 and Figure 1b, drive the time evolution (red arrows) of chromatin from an initial configuration (say corresponding to the red dot in Figure 1b) into the two energy minimal wells corresponding to the two chromatin phases. The gradients of the free energy in the 𝜙 𝑛 -𝜙 𝑑 variable space (shown by the contour plot in Figure 1b) provide the driving force for the time-evolution of chromatin organization towards the steady-state. This driving force is called the chemical potential and is written using variational principles as, 𝜇 𝑛 (𝒙, 𝑡) = 𝛿𝑊 𝛿𝜙 𝑛 = 𝜕𝑊 𝜕ϕ 𝑛 -𝛁 • ( 𝜕𝑊 𝜕∇ϕ 𝑛 ) 𝜇 𝑑 (𝒙, 𝑡) = 𝛿𝑊 𝛿𝜙 𝑑 = 𝜕𝑊 𝜕ϕ 𝑑 -𝛁 • ( 𝜕𝑊 𝜕∇ϕ 𝑑 ) (𝑆3)

where, 𝜇 𝑛 is the chemical potentials of nucleoplasm driving its kinetics as described in Section S1. 4. And 𝜇 𝑑 is the chemical potential for the order parameter 𝜙 𝑑 evolving the epigenetic marks in a conserved manner as discussed in the sub-section S1.5.

Before discussing the role of chemical potentials 𝜇 𝑑 and 𝜇 𝑛 in spatiotemporal evolution of chromatin organization, we describe the specific form of energetic contributions 𝑊 𝐶𝐶𝐼 (𝜙 𝑑 , 𝜙 𝑛 ) arising from the chromatin-chromatin interactions in Eq S1. The chosen form of energetic contribution must account for a ternary mixture of the three nuclear constituents, with volume fractions -𝜙 𝑒 , 𝜙 ℎ and 𝜙 𝑛 . An energy landscape can suitably be constructed to obtain two coexisting phases:

(i) Water rich euchromatin phase with 𝜙 𝑛 > 𝜙 𝑒 , 𝜙 ℎ , and 𝜙 ℎ → 0.

(ii) Water poor heterochromatin phase with 𝜙 ℎ > 𝜙 𝑒 ,𝜙 𝑛 and 𝜙 ℎ → 𝜙 ℎ 𝑚𝑎𝑥 . Also, 𝜙 𝑛 is very small.

The Flory Huggins model for the energetics of a mixture of polymer in a solvent can be used as a description of the energy landscape. It incorporates the competition between the change in entropy due to the mixing process and the change in enthalpy due to interactions such as bond formation between the polymer and the solvent. The change in free energy density of a ternary mixture separating into two phases, accounting for entropic and enthalpic contributions is given as [7, 8],

𝐶𝐶𝐼 (𝜙 𝑑 , 𝜙 𝑛 ) = 𝑘 𝐵 𝑇 Ω [𝜙 𝑒 2 + 𝜙 ℎ 𝑁 ln 𝜙 ℎ + 𝜙 𝑛 ln(𝜙 𝑛 ) ⏟ Entropic Contribution + 𝜒𝜙 ℎ 𝜙 𝑛 ⏟ Enthalpic Contribution ] (𝑆4)

where, 𝑁 is the degree of polymerization. For a particularly large polymer such as chromatin in the nucleoplasm, 𝑁 → ∞ thereby lowering the entropic contribution to the free energy. On a (𝜙 𝑑 , 𝜙 𝑛 ) phase space, the contour plot of the energy density 𝑊(𝜙 𝑑 , 𝜙 𝑛 ) is shown in Figure S2, left. Also note that, 𝑘 𝐵 is the Boltzmann constant, T is the temperature and Ω is the volume of individual chromatin particles. The coefficient 𝑘 𝐵 𝑇/Ω acts as an energy scaling factor.

The Flory-Huggins form of free energy is appropriate when considering only the enthalpic and entropic contributions to the free energy. It has been shown that locally driven energy consuming/producing activity, such as motor activity or transcription, or ATP consumption, can drive an 'activity-induced' segregation of the chromatin phases [9]. To account for such, more general, mechanisms of phase-separation it becomes more suitable to adopt a generalized simplified biquadratic form of energy landscape. Several polynomial descriptions of double-well energy functions emulating the Flory-Huggins form have been proposed and utilized previously [10-13]. To capture the biphasic (hetero-and euchromatin) phase-separation of a ternary mixture (𝜙 𝑒 , 𝜙 ℎ , 𝜙 𝑛 ) we adopt a simplified biquadratic double-well energy as, The contour plot of the adopted free energy form (Eq S5) on a (𝜙 𝑑 , 𝜙 𝑛 ) phase space (Figure S2, right), shows a similar location of two energy wells as the Flory-Huggins description. However, the simplification additionally allows (a) easier numerical implementation as the log terms may become undefined due to numerical errors, and (b) a better control on the location of wells, since the well locations are directly defined by the value of the parameter 𝜙 ℎ 𝑚𝑎𝑥 . Due to these advantages, we use the free energy description given by Eq S5.

The dynamic evolution of the nucleus towards the steady-state is governed in a time-dependent fashion by a combination of diffusion and reaction kinetics (Figure 1b, bottom panel). The local conservation of reactively inert nucleoplasm content relates the time evolution of local nucleoplasm volume fraction with its nominal volumetric flux 𝓙 𝑛 as, ϕ ̇𝑛 = -𝛁. 𝓙 𝑛 . Fick's first law gives the volumetric flux of nucleoplasm as 𝓙 𝑛 = -𝑀 𝑛 ∇𝜇 𝑛 in terms of the gradient of the chemical potential of nucleoplasm where 𝑀 𝑛 denotes the mobility of nucleoplasm in the nucleus, which is related to the dissipation that occurs when water flows through the porous nuclear microstructure. Thus,

The epigenetic reaction kinetics controls the methylation or acetylation levels of the histone tail marks. Heterochromatin, which is rich in methylation marks on the histone tails, can be converted into euchromatin phase, where histones are marked by an increased acetylation level. This process encompasses the removal of methylation marks on the histone tails called demethylation via proteins classified as histone demethylase (HDM), followed by acetylation of the histone tails, via histone acetyltransferase (HAT) as shown in Figure S3. We classify both these processes together as the acetylation of the histone -occurring at a rate Γ 𝑎𝑐 -which converts heterochromatin into euchromatin. Conversely euchromatin is converted into heterochromatin by first the deacetylation (via histone deacetylase, HDAC) followed by methylation (via histone methyltransferase, HMT) at a cumulative rate Γ 𝑚𝑒 as shown in Figure S3. Distinct from the conservative evolution in Eq S7, the interconversion of chromatin phases via epigenetic reactions changes the relative content of heterochromatin and euchromatin such that,

In the reaction kinetics via Eq S7, we have considered that any heterochromatineuchromatin interconversion occurs at a constant rate. However, this does not account for the chromatin-chromatin interaction energetics due to which heterochromatin-heterochromatin or euchromatin-euchromatin neighbors are more stable over euchromatin-heterochromatin neighbors. Therefore, if acetylation results in more euchromatin-euchromatin neighbors it will be more favorable than if it results in unlike-marked neighbors. Therefore, the acetylation and methylation reactions are additionally influenced by the specific location of the reaction site within the genome. The role of energetics in driving the reaction kinetics is discussed in more detail with specific schematic example in Section S1.6. Using qualitative examples given in Section S1.6, and with precise theoretical derivation in Section S1.7, we realize that such neighborhood dependent reaction kinetics effectively emulate a diffusion-like evolution of epigenetic marks. In other words, chromatin evolves in time as if the epigenetic marks of acetylation and methylation are diffusing spatially dependent on the gradients of their chemical potential 𝜇 𝑑 . Such evolution is written as,

Here, 𝑀 𝑑 is the mobility of epigenetic marks in the nucleus. Note that mobility 𝑀 is related to the diffusivity as 𝐷 = 𝑀𝑘 𝐵 𝑇 Ω .

Lastly, transcription by RNAPII activity requiring ATP-based energy expenditure, is known to supercoil chromatin fiber thereby resulting in active extrusion of DNA through cohesin rings. Since extruded chromatin loops are transcriptionally active, such extrusion within the euchromatin phase does not alter the gene expression. However, extrusion of transcriptionally silent heterochromatin into chromatin loops switches the transcriptional state of chromatin pulling out the silenced genes near the heterochromatin-euchromatin boundary into the euchromatin phase (Figure 1c). Chromatin extrusion at the phase boundaries happens in two steps (Figure 1c, 6a):

1. Cohesin rings entrap a portion of DNA fiber along the interface due to a balance between its loading (with reaction rate Γ 𝑙 ) via NIPBL/MAU2 and unloading (with reaction rate Γ 𝑢𝑙 ) via WAPL/PDS5. The overall rate of cohesin loading can be written as Γ 𝑐𝑜ℎ = Γ 𝑙 -Γ 𝑢𝑙 . 2. This is followed by the active extrusion of supercoiled loops of chromatin through the cohesin rings via the transcription due to RNAPII, at a rate denoted by Γ 𝑡𝑟 .

Altogether, Γ 𝑎 = Γ 𝑡𝑟 × Γ 𝑐𝑜ℎ denotes the overall the rate at which the chromatin extrusion converts heterochromatin into euchromatin via the two-step process (Figure 1c). Thus, we can write the transcriptionally-dependent conversion of the chromatin phases as,

Where, the exponential factor ensures that the transcription-based chromatin extrusion is spatially restricted to a narrow region where the volume fraction of heterochromatin is

⁄ + Δ𝜙, while peaking at the interface (𝜙 ℎ = ϕ h max /2). Combining Eq S7-S9, the timeevolution of the order parameter 𝜙 𝑑 can be written as,

Equations S3, S6 and S10 together form the time-dependent mathematical set of governing equations describing the spatiotemporal evolution of chromatin organization in the nucleus. As a ready reference, we have here listed the mathematical symbols used in our model, their physical meaning, and units in Table S1.

Table S1: List of mathematical symbols Symbol Physical Interpretation SI Unit Remarks 𝒙 Spatial variable, two-dimensional vector m 𝑡 Time variable s Volume Fractions 𝜙 ℎ Volume fraction of heterochromatin 1 Defined in Section S1.1 𝜙 𝑒 Volume fraction of euchromatin 1 𝜙 𝑛 Volume fraction of nucleoplasm 1 𝜙 𝑑 = 𝜙 ℎ -𝜙 𝑒 , difference between heterochromatin and euchromatin volume fractions, order parameter 1 Energetic contributions 𝑊 Total free energy density of chromatin organization 𝐽/𝑚 3 Defined via Eq S1 𝑊 𝐶𝐶𝐼 Free energy density due to chromatin-chromatin interactions 𝐽/𝑚 3 Defined via Eq S5 𝑉 𝐿 Strength of chromatin-lamina interactions 𝐽/𝑚 3 Defined in Section S1.2 𝑑 Distance of a point from the nuclear periphery m 𝑑 0 Length-scale of chromatin-lamina interactions -the distance over which effect of chromatin-anchoring proteins vanishes m 𝜂 Energy penalty on formation of sharp interface 𝐽/𝑚 Diffusion kinetics 𝜇 𝑛,𝑑 Chemical potential = 𝛿𝑊/𝛿𝜙 𝑛,𝑑 -driving force for energy-reducing passive diffusion kinetics 𝐽/𝑚 3 Defined via Eq S3 𝓙 𝑛,𝑑 Nominal volumetric flux of mobile species -defined as volume of species flowing per unit local area per unit time 𝑚 3 /𝑚 2 𝑠 Defined in Section S1.4 𝑀 𝑛,𝑑 Mobility of nucleoplasm or epigenetic marks in the nucleus 𝑚 2 /𝑠 𝐽/𝑚 3 𝐷 𝑛,𝑑 Diffusivity of nucleoplasm of epigenetic marks in the nucleus 𝑚 2 /𝑠 = 𝑀𝑘 𝐵 𝑇/Ω Reaction kinetics Γ 𝑚𝑒 Reaction rate of methylation, defined as deacetylation + methytransferase activity 𝑠 -1 Defined in Section S1.4 Γ 𝑎𝑐 Reaction rate of acetylation, defined as demethylation + acetyltransferase activity 𝑠 -1 Γ 𝑎 Rate of chromatin extrusion driven by transcriptionmediated supercoiling 𝑠 -1

The governing equations derived in the previous sections reveal intrinsic length and time scales within the model, which we next use to obtain the rescaled, non-dimensional set of governing equations.

The reaction-diffusion kinetics from Eq S10 results in a characteristic length scale determined by the reaction-diffusion kinetics i.e ℓ 𝑅𝐷 = √𝐷/Γ 𝑎𝑐 = √𝑀𝑘 𝐵 𝑇/ΩΓ 𝑎𝑐 . Interestingly, it can also be seen that another intrinsic length scale emerges from the competition between the interfacial and bulk mixing energies from Eq S1, i.e. the width of the interface ℓ 𝑖𝑛𝑡 = √𝜂Ω/𝑘 𝐵 𝑇. In our simulations, and from the theory discussed later in Section S4, we find that the reaction-diffusion length plays a significant role in determining the heterochromatin domain sizes and their spacing. Therefore, we choose to rescale all lengths with respect to ℓ 𝑅𝐷 , such that 𝒙 ̃= 𝒙/ℓ 𝑅𝐷 . Further, the reaction rates offer an intrinsic time scale for the system of equations, such that all times are rescaled as 𝑡 ̃= 𝑡Γ 𝑎𝑐 . Lastly, 𝑘 𝐵 𝑇/Ω which is the coefficient of energy of chromatin phase interactions in Eq S5 provides the energy scaling such that all energy densities are written as 𝑊 ̃= 𝑊Ω/𝑘 𝐵 𝑇.

Rescaling Eq S1,

̃= [𝜙 𝑒 2 + 𝜙 ℎ 2 (ϕ h max -𝜙 ℎ ) 2 ] ⏟ chromatin-chromatin interactions -𝑉 ̃𝐿𝜙 ℎ 𝑒 -𝑑 𝑑 0 ⏟ chromatin-lamina interactions + 1 2 𝜂Ω 𝑘 𝐵 𝑇 ΩΓ 𝑎𝑐 𝑀𝑘 𝐵 𝑇 |∇𝜙 𝑛 | 2 + 1 2 𝜂Ω 𝑘 𝐵 𝑇 ΩΓ 𝑎𝑐 𝑀𝑘 𝐵 𝑇 |∇𝜙 𝑑 | 2 ⏟ Interfacial energy

Here, 𝑉 ̃𝐿 = 𝑉 𝐿 Ω/𝑘 𝐵 𝑇 is the rescaled strength of chromatin-lamina anchoring interactions. Note that the coefficient of the interfacial energy terms can be rewritten in terms of the ratio of the length scales ℓ 𝑖𝑛𝑡 ℓ 𝑅𝐷 = 𝛿. Physically, the parameter 𝛿 is a rescaled measure of the width of the interface. Thus,

The rescaled chemical potentials from Eq S3 are,

Expanding these,

Lastly, we rescale the kinetics equations such that,

Note that in the second part of Eq S13, all reaction rates have also been rescaled with respect to the time scale such that Γ ̃𝑚𝑒 =

. To solve the rescaled set of governing equations, Eq S12 and Eq S13 are solved together subjected to boundary conditions of no flux of nucleoplasm or epigenetic marks across all boundaries. For the purposes of numerical implementation, these equations are converted into weak-form suitable for a suitable for finiteelement solver implementation. We then used COMSOL Multiphysics with a 'Weak Form PDE' module for the solution of the equations. The boundary conditions, and the list of non-dimensional parameters used in our model are discussed in detail in Section S8.

As discussed in the previous section, we have incorporated the kinetics of both diffusive and reactive nature -the former being conservative, i.e. it does not change the net amount of heterochromatin and euchromatin in the nucleus, while the latter non-conservative since it allows interconversion of the two phases (Figure S3). Here we explain in detail the kinetics underlying the non-conservative and conservative dynamics of nucleosomes in the chromatin polymer.

Firstly, we consider the diffusion of nucleoplasm via the conservative kinetics given by Eq S6. If there is no flux of water occurring across the nuclear lamina, the net amount of nucleoplasm in the nucleus 𝜙 ̅ 𝑛 remains a constant. Locally the effect of nucleoplasm movement allows methylated histones to come together resulting in their compaction (Figure S4a). Thus, the conservative diffusion of nucleoplasm allows coarsening of heterochromatin domains keeping the euchromatin phase water rich.

Before discussing the effective diffusion of epigenetic marks via the conservative kinetics (Eq S7), we focus on the epigenetic reaction kinetics described via Eq S8. The role of epigenetic reaction kinetics is non-conservative since it allows an interconversion of euchromatin to heterochromatin via a reaction rate Γ 𝑚𝑒 and heterochromatin to euchromatin via a reaction rate Γ 𝑎𝑐 (Figure S4b). Note that these interconversions require multiple steps, such as demethylation followed by acetylation, or deacetylation followed by methyltransferase activity (Figure S3). While these reactions do not affect the total amount of water or DNA in the nucleus, they change the individual amounts of heterochromatin and euchromatin. For instance, increase in Γ 𝑎𝑐 (Γ 𝑚𝑒 ) will increase the overall euchromatin (heterochromatin) in the nucleus, without changing the total amount of chromatin. The presence of epigenetic reactions has an additional energetics influenced component. When considering the reaction kinetics due to epigenetic regulation we have to account for the fact that the nucleosomes that are in close proximity to each other can interact. These interactions might stem from histone-bridging proteins such as HP1 [14, 15], binding complexes such as SAGA [16], or even ionic interactions due to the presence of post-translational modifications. The interactions are captured in our model as chromatin-chromatin interactions and give rise to the energy landscape described in Figure 1b. Because of chromatin-chromatin interaction energetics, heterochromatin-heterochromatin or euchromatin-euchromatin neighbors are more stable than euchromatin-heterochromatin neighbors. Thus, the acetylation and methylation reactions effectively depend on the specific location of the nucleosome within the genome.

The role of energetics can be better understood via the schematic shown in Figure S4d. Say a portion of chromatin polymer at some point in time looks like the initial configuration shown in Figure S4d(i) (in the figure, heterochromatin is shown as red and euchromatin blue). Without loss of generality, let us say that the first step is conversion of heterochromatin into euchromatin (acetylation). There are various heterochromatin sites which can be converted into euchromatin (shown in the grey box), but the configuration shown in Figure S4d (ii) has the highest probability of occurring. This is because it maximizes like-nucleosome neighbors and thus reduces the energy of chromatin-chromatin interactions the most. The next step could either involve converting heterochromatin to euchromatin or vice versa. If a heterochromatin to euchromatin flip occurs, in our meso-scale model its effect is captured by the non-conservative reaction (Eq S8).

On the other hand, if a euchromatin to heterochromatin flipping occurs, the possible configurations attainable are shown in the grey box. The configuration with highest probability of occurrence is if all heterochromatic nucleosomes are segregated from the euchromatic ones, as shown in Figure S4d (iii). Thus, over the course of these two steps, we see that an effective conservative rearrangement of heterochromatin and euchromatic nucleosomes can occur which leads to coarsening of the two phases. This conservative evolution effectively is captured in our model as 'diffusion of epigenetic marks'.

Note that the sequence of events described here gives rise to effective diffusion of epigenetic marks along the chain of polymer -labelled cis-diffusion in Figure S4c. A similar reaction kinetics driven epigenetic diffusion can occur between nucleosomes that are close in the 3D space, although not neighbors on the chromatin polymer chain. Such an event, as shown in Figure S4c, can be called trans-diffusion, and can similarly occur.

To support this qualitative description of events, in the next section we show how the presence of reaction rates dependent on the energetic interactions of nucleosomes can give rise to effective 'diffusion of epigenetic marks'.

In this subsection, we will theoretically prove that when the epigenetic reaction kinetics becomes neighborhood dependent, as discussed in subsection S1.6, we recover the spatiotemporal evolution described by Eq S10, used in our continuum chromatin dynamics model. To do so, first we will explain the stochastic events describing the conservative diffusive and non-conservative reaction kinetics.

Simple diffusion: In a 1D setting (to simulate a polymer), consider the series of empty sites that a randomly walking particle can occupy, as shown in Figure S4.1. Let us focus on a site at location 𝑥 with the left and right neighbors located at 𝑥 -𝑑𝑥 and 𝑥 + 𝑑𝑥.

Let 𝑃(𝑥) be the probability that the site at location 𝑥 is occupied. For brevity, the probability that the left site is occupied, 𝑃(𝑥 -𝑑𝑥) is denoted as 𝑃 -, while the probability that the right site is occupied 𝑃(𝑥 -𝑑𝑥) is denoted as 𝑃 + .

There are four possible events that can occur (Figure S4.1) -if the site at position 𝑥 is occupied, the particle can either move to the left or right given that those sites are empty, or if the site at position 𝑥 is free particle from left or right neighboring site can move in, provided that those sites are occupied. Then the probability of particle occupying site 𝑥 evolves in time 𝛿𝑡 as,

Expanding the left and right probabilities via Taylor expansion to the third order,

Thus, the spatiotemporal evolution of randomly walking particles follows the diffusive kinetics.

Reaction kinetics: In a 1D setting, now consider that the nucleosomes can be of two flavors: heterochromatic (H) and euchromatic (E). 𝑃 𝐻 is the probability that the nucleosome is heterochromatic. The probability that nucleosome is euchromatic is 𝑃 𝐸 = 1 -𝑃 𝐻 . The two nucleosome species can be interconverted at a rate Γ 𝑎𝑐 or Γ 𝑚𝑒 (dimension of 𝑠 -1 ). Note that nucleosomes are not moving in space, only their flavor changes in time. For now, let us assume that these conversion rates are constant. Then the probability that nucleosome at site 𝑥 is heterochromatic 𝑃 𝐻 (𝑥, 𝑡) evolves in time as,

It should be noted that Eq S15 is analogous to the non-conservative reaction kinetics in our continuum model as given by Eq S7. At steady state ( 𝜕𝑃 𝐻 𝜕𝑡 = 0), the probability of seeing a heterochromatic nucleosome at a site is 𝑃 𝐻 = Γ 𝑚𝑒 Γ 𝑎𝑐 +Γ 𝑚𝑒 . This is again analogous to the theoretical average heterochromatin content in the nucleus, derived from our model (see Section S3).

Reaction kinetics with neighborhood dependent energetics: Now let us consider the case where the rate of acetylation and methylation are not constant, but dependent on whether the neighbors are similar or dissimilar. This can happen because the conversion requires breaking bonds between neighbors. Say if the nucleosome is H (heterochromatic), and needs to convert to E (euchromatic), there are four possible options:

Configuration Change in energy Acetylation rate 1 HHH Need to break two heterohetero bonds ΔE = 2E 𝐻𝐻 Γ 𝑎𝑐 2 HHE Need to break one heterohetero and one eu-hetero bond ΔE = E 𝐻𝐻 + 𝐸 𝐸𝐻 ∼ Γ 𝑎𝑐 + ΔΓ 3 EHH Figure S4.2: Reaction Kinetics. 4 EHE Need to break two eu-hetero bonds ΔE = 2𝐸 𝐸𝐻 ∼ Γ 𝑎𝑐 + 2ΔΓ

Since like-marked nucleosomes interact strongly as compared to unlike marked nucleosomes, 𝐸 𝐸𝐻 < 𝐸 𝐻𝐻 , 𝐸 𝐸𝐸 . Thus, acetylation rate is higher (ΔΓ > 0) when a bond between unlike nucleosomes is broken than if that between like-marked nucleosome is broken. Similarly, we can have four possible options for Γ 𝑚𝑒 . Considering these eight possible options after taking into account the effects of the neighbors, the probability that nucleosome at site 𝑥 is heterochromatic evolves in time as,

Here, 𝑃 𝐻(𝐸) +(-) is the probability that the right (left) neighbor is heterochromatic (euchromatic). For simplification, we begin with the first bracket term on the right side of Eq S16:

(Next, we use the fact that 𝑃 𝐻 + 𝑃 𝐸 = 1, 𝑃 𝐻 + + 𝑃 𝐸 + = 1 and 𝑃 𝐻 -+ 𝑃 𝐸 -= 1.)

Similarly, the second bracket term can be simplified, thereby reducing Eq S16 as,

Note that if the neighborhood effect is switched off (i.e., ΔΓ = 0), Eq S17 reduces to Eq S15 derived above for the case of constant reaction rates. We can simplify Eq S17 further as,

To see the diffusion-like similarity of the second, we can expand 𝑃 𝐻 +/-via Taylor series expansion to obtain,

Comparing the first term to the Eq S15, we see that it obeys non-conservative reaction-like kinetics, while comparing the second term to Eq S14 we see that it obeys a conservative diffusionlike kinetics. Thus, we immediately see that the evolution decouples naturally into a reaction-like term and a diffusion-like term. This happens purely because the reaction rates are asserted to be neighborhood dependent.

Comparing Eq S18 with the evolution equation used in our model (Eq S10 or Eq 2b) clarifies the mechanistic origin of 'diffusion of epigenetic marks', and proves how epigenetic reaction kinetics influenced by energetics of chromatin-chromatin interactions can give rise to effective diffusionlike kinetics of the epigenetic flavors.

We can use the physical dimensions estimated in Section S9 to get an approximate scale of the energetic parameter Δ𝐸 and rate parameter ΔΓ. We see later in Section S4, that the diffusion prominently plays a role in the euchromatin region, where in association with the methylation, a reaction-diffusion flux is formed which we predict (see Section S4 for details) will drive the sizescaling of the heterochromatin domains. We estimated the effective diffusivity (see Table S3) to be of the order ∼ 10 -3 𝜇𝑚 2 𝑠 -1 . Considering the nucleosome spacing Δ𝑥 ∼ 50 nm within the euchromatin region, we obtain ΔΓ ∼ 2 𝑠 -1 . Compare this to Γ 𝑎𝑐 ∼ 10 -2 𝑠 -1 , we see that the heteroto-euchromatin conversion rate increases if the heterochromatin nucleosome is surrounded by one or more euchromatic nucleosomes. In other words, formation of euchromatin would be more stable in the vicinity of other euchromatin. Based on the definition of ΔΓ above, we can approximate ΔΓ Γ 𝑎𝑐 ∼ 𝑒 (𝐸 𝐻𝐻 -𝐸 𝐸𝐻 )/𝑘 𝐵 𝑇 . Based on our estimates of ΔΓ and Γ 𝑎𝑐 , thus we can get 𝐸 𝐻𝐻 -𝐸 𝐸𝐻 ∼ 5.29𝑘 𝐵 𝑇.

MATLAB was used for the analysis of STORM images. For chromatin density qualification, Voronoi tessellation-based segmentation was implemented to construct Voronoi polygon of each localization [17]. The corresponding Voronoi polygon of each locus represents a specific region where any points within this region are closer to this locus, such that the size of Voronoi polygon is inversely proportional to the local Voronoi density. Voronoi polygons located at the edge were omitted due to nearly infinite Voronoi area. The Voronoi density map was further constructed through calculating the reciprocal of Voronoi area map. To differentiate between hetero-and euchromatic regions, a density-based threshold was applied to filter out low density euchromatic region. The threshold density value was chosen such that 60~100 percentile of the Voronoi density distribution were classified as heterochromatin in the control group. The same density threshold was applied to all other treatments in comparison with the control. The remaining heterochromatin points cloud was then clustered using Density-based spatial clustering of applications with noise (DBSCAN) algorithm [18]. The heterochromatic point clouds were clustered into separated subdomains, such that each subdomain includes all the neighboring and connected Voronoi polygons which are above threshold. The hyperparameters of DBSCAN were set such that the minimal scan points are greater or equal to the dimension of dataset plus 1, namely 3 in our case [18].

Chromatin clusters obtained after the DBSCAN subdomain classification were categorized as LADs and non-LADs depending on proximity to the boundary of the nucleus image. The nucleus shape was detected using boundary detecting algorithm. A characteristic radius of nucleus (R) was then calculated. The minimal distance between each heterochromatin subdomain and nucleus boundary was calculated, such that any subdomain having a distance smaller than 0.15𝑅 is classified as LADs domain [19].

To quantify the size of non-LAD heterochromatin subdomains, the area of each non-LAD domain was calculated through detecting the boundary of point clouds which gives a polygon enveloping it. An approximating domain radius was calculated by assuming the subdomain to be a circular shape (𝑟 ≈ √Domain Area/𝜋).

The local LADs thickness was measured by sampling the LAD boundary along the nucleus periphery.

BJ cell lines (ATCC Manassas, VA) were cultured in Minimum Essential Media (ThermoFisher Scientific, Waltham, MA,#11095080) at physiological conditions (5% CO 2 and 37 °C). Cells were seeded on 35-mm glass-bottom Petri dishes (MatTek Corp.) until approximately 40-50% confluent and were given at least 24 hours to adhere to the dish before fixation. For ChromSTEM sample preparation, the previously published protocol was adapted [20]. Before fixation, cells were thoroughly rinsed three times in Hank's balanced salt solution without calcium and magnesium (EMS). Cells were fixed using 2.5% EM grade glutaraldehyde, 2% paraformaldehyde, 2 mM CaCl2 in 0.1 M sodium cacodylate buffer, pH = 7.4 at room temperature for 5 minutes and then replaced with fresh fixative and fixed on ice for an hour. The cells were then washed with 0.1 M sodium cacodylate buffer 5 times on the ice. The cells were then incubated in a blocking buffer (10 mM glycine, 10 mM potassium cyanide in 0.1 M sodium cacodylate buffer, pH = 7.4) for 15 minutes, followed by staining with 10 µM DRAQ5™ (Thermo Fisher) and 0.1% saponin solution in 0.1 M sodium cacodylate buffer, pH = 7.4 for 10 minutes. After washing with the blocking buffer twice, the cells were incubated in the blocking buffer on ice before photo-bleaching. During photobleaching on a cold stage using continuous epi-fluorescence illumination (150 W Xenon Lamp) with Cy5 red tilter with a 100x objective for 7 minutes, the cells were incubated in 2.5 mM of 3-5′-diaminobenzidine (DAB) solution (Sigma Aldrich) in 0.1 M sodium cacodylate buffer, pH = 7.4. The cells were washed with 0.1 M sodium cacodylate buffer five times and then stained with reduced osmium solution (EMS) containing 2% osmium tetroxide, 1.5% potassium ferrocyanide, 2 mM CaCl2 in 0.15 M sodium cacodylate buffer, pH = 7.4 for 30 minutes on ice. Then the cells were washed with double distilled water five times on ice. Serial ethanol dehydration (30% ethanol, 50%, 70%, 85%, 95%, An ultramicrotome (UC7, Leica) was used to section 100 nm thick slices that were deposited onto a copper slot grid with carbon/Formvar film. Then, 10 nm colloidal gold fiducial markers were carefully deposited on both sides of the sample. A 200 kV cFEG STEM (HD2300, HITACHI) with HAADF mode was used and while keeping the field of view constant, the sample was tilted from - 60° to 60° with 2° increments on two roughly perpendicular axes. The fiducial markers were used to align the tilt series in IMOD [21] and reconstructed using TomoPy [22] with a penalized maximum likelihood for 40 iterations independently.

The centers for individual chromatin domains were estimated from local maxima obtained from ChromSTEM projection with enhanced contrast in FIJI as previously described [23]. Domains occupying less than 50% volume in the z-plane were excluded from the analysis as they could be incomplete parts of other neighboring domains. Mass scaling analysis and radial density analysis were then performed originating from the identified individual domain centers. Average mass scaling originating from the individual domain centers was estimated using the area (mass) weighted by the grayscale ChromSTEM intensity within concentric circles with increasing distances from the domain centers. Similarly, radial chromatin density was estimated as the grayscale ChromSTEM intensity within concentric circles with increasing distances from the domain centers (Figure S5). We have shown that the radial mass density originating from the domain center decreases with increasing distance and approaching the domain boundary and then increases as the boundary of neighboring domains starts interacting. Both the mass scaling behavior and the radial chromatin density profile for each domain were then utilized to obtain the boundary or the approximate radius of the domain. The mass scaling approximately follows power-law scaling up to a given length scale and can be represented by a given slope or scaling exponent, D based on the linear regression of the mass scaling curve in the log-log scale. Beyond the domain regime, the slope gradually increases to reach the supra-domain regime. The radius of chromatin packing domains was estimated as the smallest length scale where the mass scaling curve deviates from the initial power law calculated from small length scales by 5% and the radial chromatin density starts to increase after a gradual decrease. The distributions for packing domain radius and density are shown as mean ± S.D. using violin super plots [24].

Prior to imaging, cells were cultured in 35-mm glass-bottom Petri dishes until approximately 70% confluent. All cells were given at least 24 hours to re-adhere before treatment (for treated cells) and imaging. HCT116, A549, HeLa, and U2OS cells were treated with Actinomycin D (Gibco, Cat: 11805017) for 1 hour at a final concentration of 5 µg/mL.

For live-cell measurements, cells were imaged and maintained under physiological conditions (5% CO2 and 37°C) using a stage-top incubator (In Vivo Scientific, Salem, SC; Stage Top Systems).

The PWS optical instrument was built on a commercial inverted microscope (Leica, Buffalo Grove, IL, DMIRB) supplemented with a Hamamatsu Image-EM CCD camera C9100-13. This camera was coupled to an LCTF (CRi Woburn, MA) for hyperspectral imaging. Spectrally resolved images of live cells were collected between 500 and 700 nm with a 2-nm step size. Broadband illumination was provided by an Xcite-120 light-emitting diode lamp (Excelitas, Waltham, MA). PWS is a highthroughput, label-free approach that measures the spectral standard deviation (Σ) of internal optical scattering originating from nuclear chromatin. The variations in the refractive index distribution Σ, are characterized by a mass density autocorrelation function (ACF) to calculate chromatin packing, scaling D.

Based on ChromSTEM [23], we have previously reported that chromatin packs into domains, wherein each domain exhibits a polymeric fractal-like behavior and can be described by an average packing scaling exponent. This implies that within the fractal regime, the genomic size of chromatin scales with its physical size following a power law relationship. Therefore, we estimated the upper bound of the power law regime as a measure of domain size. Thus, a power-law ACF which incorporates a lower and upper length scale limit of the power law regime was utilized for the subsequent approximations, where Γ(𝑥,𝑎) is the upper incomplete gamma function, and 𝑙min and 𝑙max characterize fractality's lower and upper length scales, respectively; 𝐵𝜌(𝑟=0) is 𝜎𝜌 2 , the variance of chromatin mass density; 𝐷𝐵 describes the shape of 𝐵𝜌 and is related to 𝐷, and r is the spatial separation Utilizing this previously described methodology [25], we evaluated 𝑙max, the upper length scale of chromatin mass density scaling to estimate the relative size of domains upon ActD treatment.

Our model predicts that the heterochromatin domains obtained at the steady state display a characteristic size. The length scale of the stable domains is regulated in tandem by the epigenetic reactions -acetylation as well as methylation -and the transcriptionally active extrusion of chromatin loops. As the levels of histone acetylation are increased (or conversely methylation is decreased), as shown in Figure S6a, we see that the sizes of the heterochromatin domains decrease. This trend is displayed not just by the interior heterochromatin domains, but also by the LADs near the nuclear periphery. Also note that the morphology of the heterochromatin domains is regulated by the balance of acetylation and methylation reaction kinetics. Under the conditions where acetylation predominates over methylation, the domains formed are more circular. On the other hand, when methylation is increased while acetylation is reduced, the domains become larger and predominantly lamellar.

In addition to methylation and acetylation, transcriptional activity also regulates the sizes of the heterochromatin domains (Figure S6b). As the rate of supercoiling-driven chromatin extrusion reduces, the heterochromatin domains become larger in size. This is discussed expansively in Sections 3.2 and 3.3 of the main manuscript (also see Figure 2d-g).

Here we show that the total amount of chromatin which falls into the individual phases, i.e., either euchromatin or heterochromatin, is determined solely by the epigenetic and chromatin extrusion reaction kinetics. The spatiotemporal evolution of order parameter 𝜙 𝑑 due to the presence of diffusion of epigenetic marks and reaction driven interconversion of eu-and heterochromatin phases, as described by Eq 2b (or equivalently, Eq S10) as,

Once steady state is reached, the evolution halts and 𝜕𝜙 𝑑 𝜕𝑡 = 0. At this stage, the relation between the average euchromatin content 𝜙 ̅ 𝑒 and the average heterochromatin 𝜙 ̅ ℎ content in the nucleus can be obtained by averaging Eq S19 (i.e., integrating over the entire nucleus or region of interest and divide by the area of the nucleus) as,

Note that the term ∇ 2 𝜇 𝑑 on the left-hand side is non-zero near the domain boundaries and after averaging over the entire volume of the nucleus can be approximated to zero. Thus,

The last term on the right-hand side can be resolved via integration by parts as,

The last term (in blue) when integrated over the domain measures the length of the interfaces between the heterochromatin and euchromatin domains. The term in red is non-zero only along the interface. For a narrow interface width Δ𝜙 → 0, this can be approximated as multiplying the integrand by a factor depending on the length of the interface giving a parameter ℓ 𝑖𝑛𝑡 . Thus,

Thus, the equation can be written as,

𝜅 being a non-trivial function dependent on ϕ h max , volume fraction change across the interface Δ𝜙, and the length of the interface between the two chromatin phases. By definition of the volume fractions, 𝜙 𝑒 = 1 -𝜙 𝑛 -𝜙 ℎ . Thus,

Or,

In the absence of transcription (Γ 𝑎 = 0), we can instead write the average heterochromatin and euchromatin contents as,

The average eu-and heterochromatin contents obtained via Eq S20 and S21 can be rescaled to obtain:

The steps involved in theoretical derivation of heterochromatin domain size have been enumerated in the main text (Section 3.2). Here we show the complete derivation of domain size determination in the interior of the nucleus, away from the periphery. To analyze the steady state size of the heterochromatin domain, we first examine the volume fraction fields within and around the droplet. As discussed in Section S3, the acetylation, methylation and chromatin extrusion together determine the mean heterochromatin (and euchromatin) volume fraction in the nucleus, given by Eq S20. In absence of any energetic considerations, this would give rise to a homogeneous mean chromatin composition (𝜙 ̅ ℎ , 𝜙 ̅ 𝑒 ). However, this composition (Figure 2b; light blue circle) lies in neither of the energy wells and is thus energetically unfavorable. As the system marches towards a steady state, its free energy must reduce, requiring the phase separation to initiate via nucleation of heterochromatin droplets (Figure S7). Under a dilute limit, i.e., when there is a lot more euchromatin than heterochromatin, we can assume that the droplet size is much smaller than the length scale of the interdomain spacing such that neighboring droplets are far enough to not interact with each other. Under such assumption the heterochromatin distribution would be spherically symmetric. Using a polar coordinate system with origin at the droplet center, we can determine radial distribution of the heterochromatin volume fraction 𝜙 ℎ (𝑅) as, 𝜙 ℎ (𝑅) = { 𝜙 ℎ 𝑚𝑎𝑥 𝑅 < 𝑅 𝑑 (driven by phase-separation) 𝜙 ℎ + = 𝜂 𝑅 𝑅 = 𝑅 𝑑 + (driven by interfacial surface tension) 𝜙 ̅ ℎ 𝑅 = ∞ (region far from the droplet stays undisturbed)

where 𝑅 𝑑 is the radius of the droplet at the current instance, while 𝜂 measures the interfacial energy. Figure S7 shows the distribution of heterochromatin volume fraction 𝜙 ℎ (𝑅) around a single spherical domain of heterochromatin (red) of radius 𝑅 𝑑 as it grows surrounded by euchromatin phase (in blue). The slope of the heterochromatin profile outside the droplet will drive an inward flux of heterochromatin into the droplet. The volume fraction field outside the droplet at steady state must follow the equation,

2 𝜙 ℎ -(Γ 𝑎𝑐 + Γ 𝑎 )𝜙 ℎ + Γ 𝑚𝑒 𝜙 𝑒 with boundary conditions 𝜙 ℎ | 𝑅 𝑑 + = 𝜂/𝑅 and 𝜙 ℎ | ∞ = 𝜙 ̅ ℎ , and thus must have the form, 𝜙 ℎ 𝑜𝑢𝑡 (𝑅) = 𝜙 ̅ ℎ + (𝜙 ℎ + -𝜙 ̅ ℎ ) 𝑅 𝑑 𝑅 𝑒 𝑅 𝑑 -𝑅 ℓ 𝑅𝐷 where ℓ 𝑅𝐷 is the characteristic reaction-diffusion length scale given under a dilute limit as ℓ 𝑅𝐷 = √ 𝐷 ℎ Γ 𝑎𝑐 . Note that the transcription driven active extrusion occurs only at the periphery and thus does not play a role in the reaction-diffusion length scale. Thus, 𝜕𝜙 ℎ 𝑜𝑢𝑡 𝜕𝑅 = (𝜙 ̅ ℎ -𝜙 ℎ + ) 𝑅 𝑑 (ℓ 𝑅𝐷 + 𝑅) ℓ 𝑅𝐷 𝑅 2 𝑒 𝑅 𝑑 -𝑅 ℓ 𝑅𝐷 Thus, 𝜕𝜙 ℎ 𝑜𝑢𝑡 𝜕𝑅 | 𝑅→𝑅 𝑑 = 𝜙 ̅ ℎ 𝑅 𝑑 -𝜂 𝑅 𝑑 2 (𝑤ℎ𝑒𝑛 𝑅 𝑑 ℓ 𝑅𝐷 ≪ 1) (𝑆24)

At the periphery of the droplet, the diffusive influx of heterochromatin into the droplet due to the reaction-diffusion phenomena outside is,

The inward diffusion is opposed by the outward flux of heterochromatin into euchromatin phase which occurs due to both acetylation of histones inside as well as supercoiling-driven chromatin extrusion along the domain boundary. Thus, the rate of change of the volume of the droplet 𝑉 𝑑 can be written as,

where ℓ 𝑖𝑛𝑡 is the width of the interface between the chromatin phases. Simplifying this we obtain,

neighboring heterochromatin domains do not interact with each other. Such assumption therefore requires that the domains be separated by a distance which scales with the reaction-diffusion length scale such that the spacing between the domains can be written as,

Thus, a quantitative image analysis of super-resolution images of heterochromatin foci in in-vitro nucleus can be used to quantitatively estimate the parameters Γ 𝑎𝑐 and Γ 𝑚𝑒 using Eq S28. As discussed in Section S6, we will use these relations to motivate the parameter choice for our numerical simulations.

For consistency with our non-dimensional model, we can rescale Eq S27-29 to obtain a nondimensional dependence of heterochromatin domain size on the epigenetic and transcriptional kinetics. All lengths are rescaled with respect to ℓ 𝑅𝐷 , while all times with respect to 1/Γ 𝑎𝑐 . Thus, Eq S27 becomes,

While Eq S28 and S29 become,

Note that although the qualitative effect of increasing acetylation (or decreasing methylation) is the same as increasing rate of supercoiling-drive chromatin extrusion, the quantitative ways in which their effects are felt are different (different scales in Eq S25, S26). This is immediately visible from Eq S30 and S31. With Γ ̃𝑚𝑒 = Γ 𝑚𝑒 /Γ 𝑎𝑐 , the domain radius scales as, 𝑅 ̃𝑑 𝑠𝑠 ∼ √ Γ ̃𝑚𝑒 1+Γ ̃𝑚𝑒 . On the other hand, the scaling of domain radius with Γ ̃𝑎 is more complex, but if the effect of acetylation is ignored, from Eq S26, we see that 𝑅 ̃𝑑 𝑠𝑠 ∼ 1 Γ ̃𝑎. This difference in scales can be exemplified by graphically seeing the change in 𝑅 ̃𝑑 𝑠𝑠 as these rates change, as shown in Figure S9. We see that the effect of acetylation is more pronounced than the effect of chromatin extrusion. However, this result is derived theoretically for a single domain. Since a single domain has only one continuous boundary where transcription occurs, this scaling holds. In a numerical simulation, there are many boundaries where transcription will drive chromatin extrusion, and thus change the scaling derived theoretically. The numerical effect of changing epigenetic reaction rate, and extrusion rate was discussed in Section S2, Figure S5. Theoretically this can be captured by increasing the value of 𝜅 in Eq S30. The parameter 𝜅, as discussed in Section S3, includes, amongst other effects, the role of multiple domain boundaries.

As discussed in Section 3.1, the steady state organization of chromatin the nucleus (Figure 2) comprises of many disconnected domains of heterochromatin phase (𝜙 ℎ = 𝜙 ℎ 𝑚𝑎𝑥 ) with a characteristic size. We also determined that the size of these domains is determined by the reaction kinetics -the rates of acetylation Γ 𝑎𝑐 , methylation Γ 𝑚𝑒 and active chromatin extrusion Γ 𝑎 (Eq 5).

To investigate the role of reactions numerically, we allow the phase-separation to occur from the same initial state as in Figure 2a, but without any reaction kinetics. The initial, intermediate and steady state chromatin organization thus obtained is shown in Figure S10a. While intermediate steps show the nucleation of multiple domains, much like in the presence of reactions (Figure 2a, center panel). However, as the organization evolves these domains merge. At the steady state, a single domain of heterochromatin remains. The growth rate of the domains given by Eq 4 is plotted in Figure S10b. In the absence of reactions, 𝑑𝑅 𝑑 /𝑑𝑡 never goes to zero, except at the critical radius. The critical radius only ensures that the domains above this size grow, while the rest shrink. Since no stable radius is predicted, the growing domains continue to grow ( 𝑑𝑅 𝑑 𝑑𝑡 > 0), until all heterochromatin merge into a single domain. This leads to full Ostwald ripening.

We have seen via the derivation in Section S4 and Eq S25 that a steady state heterochromatin domain (obtained by setting 𝑑𝑅/𝑑𝑡 = 0) is regulated via reaction kinetics, and apparently by the interfacial energy penalty 𝜂. Here we show that the contribution of 𝜂 in determining 𝑅 𝑑 𝑠𝑠 is very small.

As in the previous section, we plot the growth rate of the heterochromatin droplets with respect to their current radius, as obtained from Eq S25 (Figure S11). The stable radius of the domains 𝑅 𝑑 𝑠𝑠 is the point where the curve intersects with the x-axis such that the growth rate of the domains becomes zero, i.e. 𝑑𝑅 𝑑 /𝑑𝑡 = 0. As the contribution of the interfacial energy is changed by changing the energetic penalty 𝜂 on the formation of interfaces, we note that the change in the domain sizes is very small. Even when 𝜂 = 0, we see that the domain size does not change appreciably.

These results highlight a key role played by the epigenetic reactions as well as transcriptional regulation of chromatin extrusion. Interestingly, these results also highlight the key difference between our 'non-equilibrium' thermodynamic phase-separation model for chromatin organization as opposed to a more traditional energy-minimizing phase-separation. The traditional phase-field models result in domain formation by a competition between the energy reducing phase separation and energy increasing interface formation. However, in our model, the role of the interface formation is overshadowed by the role of reactions in a kinetics rather than energetics driven formation and maintenance of heterochromatin domains. This 'competition' between the energetic phase-separation and kinetics of interconverting reactions results in formation of heterochromatin domains of characteristic sizes.

Like the determination of the radius of the heterochromatin droplets in the interior of the nucleus, at the nuclear periphery, the epigenetic, transcriptional and diffusion kinetics balance regulates the thickness of the LADs.

We begin by examining the volume fraction fields within and around the LAD. As discussed in Section S4, the acetylation, methylation and chromatin extrusion together determine the mean heterochromatin (and euchromatin) volume fraction in the nucleus, given by Eq S20. However, a homogenous mean chromatin composition (𝜙 ̅ ℎ , 𝜙 ̅ 𝑒 ) lies in neither of the energy wells as shown in Figure 2b (light blue circle) and is thus energetically unfavorable. Nucleation of heterochromatin domains occurs due to the reduction of free energy as the system evolves. The interaction of heterochromatin with the nuclear lamina results in a preferential nucleation of heterochromatin domains along the lamina i.e., lamina associated domains, LADs. We assume that the LADs are formed uniformly along the lamina, and can grow normal to the lamina i.e., increase in thickness, as shown in Figure S12. We also assume that the LADs are far away from the neighboring interior heterochromatin domains, and do not interact with them.

Using a cartesian coordinate system with origin at the nuclear lamina directed normal to it, we determine spatial distribution of the heterochromatin volume fraction 𝜙 ℎ (𝑥). Note that this can be easily done by extending the derivations for the interior heterochromatin domains by setting 𝑅 → ∞, to obtain a linear continuous LAD.

where ℓ 𝑅𝐷 is the characteristic reaction-diffusion length scale given under a dilute limit as ℓ 𝑅𝐷 = √ 𝐷 ℎ Γ 𝑎𝑐

. The slope of the heterochromatin profile outside the droplet will drive an inward flux of heterochromatin into the droplet. As for the interior heterochromatin domains, the inward diffusion is opposed by the outward flux of heterochromatin into euchromatin phase which occurs due to both acetylation of histones inside as well as DNA loop extrusion along the domain periphery. Thus, the rate of change of the volume of the droplet 𝑉 𝑑 can be written as,

𝐿𝐴𝐷 𝑑𝑡 = 𝑆 𝑑𝑥 𝑡 𝑑𝑡 = 𝐽 𝑖𝑛 -Γ 𝑎𝑐 𝑆𝑥 𝑡 𝜙 ℎ 𝑚𝑎𝑥 -𝑆ℓ 𝑖𝑛𝑡 Γ 𝑎 𝜙 ℎ 𝑚𝑎𝑥 2 (𝑆33) Energetic parameters: The non-dimensional energy density (Eq 1 or Eq S11) involves a single parameter -the chromatin-lamina interaction strength 𝑉 ̃𝐿. We have previously [26] shown that, 𝑉 ̃𝐿 plays a role in deciding the LAD thickness and morphology, i.e. whether the LAD would be more droplet shaped or lamellar. The LAD morphology observed in in-vitro nuclei with transcriptional abrogation is used to estimate the value of 𝑉 ̃𝐿. For the previously estimated values of epigenetic reaction rates, we parametrically vary 𝑉 ̃𝐿 (Figure S13) to obtain a close match with the LAD distribution in ActD treated nucleus (Figure 3d, right panel) as reported in Table S2. All the simulations reported here use the same value for 𝑉 ̃𝐿, hence validating its choice with the cases where transcription is active.

We consider an initial spatially homogenous distribution of chromatin and nucleoplasm in the nucleus. The nucleoplasm content of the nucleus is estimated based on experimental images as 𝜙 𝑛 initial = 0.5, which is maintained a constant in the simulations as there is no exchange exchange of water across the boundary (equivalent to a boundary condition of no outward flux of nucleoplasm (𝛁𝜇 𝑛 . 𝒏 ̂|boundary = 0)). As chromatin is confined to the nucleus, a no flux boundary condition of the order parameter (𝛁𝜇 𝑑 . 𝒏 ̂|boundary = 0) ensures the conservation of epigenetic marks.

To mimic the spatial heterogeneities of the acetylation and methylation reactions, we add a Gaussian noise with a mean as the parameter values listed in Table SI and a 20% relative standard deviation. This gives us a chromatin domain distribution in agreement with the distribution of domain sizes observed in the experimental images. A similar Gaussian noise is also added to the strength of chromatin-lamina interactions to capture the heterogeneities in the anchoring of chromatin to the lamina. Lastly, we add a random uniform perturbation to the initial chromatin configuration to represent noise due to intrinsic heterogeneities present in the nucleus.

Table S2: Values of the parameters used in simulation. Parameter Description Value Remarks Initial Conditions 𝜙 𝑛 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 Initial water (nucleoplasm) content in the nucleus 0.5 Estimated from experimental H2B density 𝜙 𝑑 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 Initial value of order parameter 0.1 The choice does not affect the simulations in presence of reactions 𝜎 𝜙 𝑛𝑜𝑖𝑠𝑒 Range of uniform perturbation due to heterogeneities in initial conditions 0.01 29 Epigenetic Kinetics Γ ̃𝑚𝑒 Reaction rate of histone methylation (non-dimensionalised) 1.2 Estimated from H2B density in nuclei with transcription inhibited 𝜎 Γ 𝑛𝑜𝑖𝑠𝑒 Variance in spatial distribution of reaction rates 0.2 Value large enough to obtain a spatial variation in domain sizes Chromatinlamina Interactions 𝑉 ̃𝐿 Chromatin-lamina interaction strength (non-dimensionalised) 0.012 Parametric study performed to best approximate the LAD distribution in absence of transcription 𝜎 𝑉 𝐿 𝑛𝑜𝑖𝑠𝑒 Variance in spatial distribution of 𝑉 ̃𝐿 = 𝜎 Γ 𝑛𝑜𝑖𝑠𝑒 Extrusion Kinetics Γ ̃a Rate of active chromatin extrusion Control 0.03 Calibrated with experiments ActD treatment 0 Total inhibition of loop extrusion WAPLΔ 0.0075 Calibrated with experiments

We have used the rescaled model (see Section S1.5) to make non-dimensional quantitative predictions on the change in chromatin organization upon transcriptional perturbations via supercoiling driven loop extrusion. We can use rescaling parameters to obtain the physical interpretation of the timescales, length scales and the energy scales involved in our model. Table S3 below lists the values of the rescaling parameters, obtained from the published literature, and thus obtained values of time, length and energy scales.

Table S3: Experimental parameters are used to translate the non-dimensional model parameters into physical units Parameter Meaning Order of magnitude Reference Remarks Reaction timescale Γ 𝑎𝑐 Rate of histone acetylation reaction ∼ 10 -2 𝑠 -foot_3 [27] Calculated as reciprocal of timescale of reaction Γ 𝑚𝑒 Rate of histone methylation reaction ∼ 10 -3 𝑠 -1 [28] Length-scale 𝐷 Diffusivity of nucleosomes ∼ 10 -3 𝜇𝑚 2 /𝑠 [29] Diffusivity of nucleosomes calculated as slope of MSD curve ℓ 𝑅𝐷 Reaction-diffusion length scale ∼ 300 𝑛𝑚 Calculated as ℓ 𝑅𝐷 ∼ √𝐷/Γ 𝑎𝑐 𝑅 𝑑 𝑠𝑠 Characteristic steady state size of heterochromatin domains ∼ 100 𝑛𝑚 [29-34] 1 Calculated as ∼ √ 3𝐷 𝛤 𝑎𝑐 𝛤 𝑚𝑒 𝛤 𝑚𝑒 +𝛤 𝑎𝑐 (Eq S28) Diffusion timescale Time for diffusion across inter-domain space ∼ 5 𝑚𝑖𝑛 Calculated as ∼ (2ℓ 𝑅𝐷 ) 2 /𝐷 (Note that this is similar to timescale of reactions) Energy scale 𝑘 𝐵 𝑇 Ω Energy scale of chromatinchromatin interactions ~10 -24 𝐽/𝑛𝑚 3 Calculated. Assuming nucleosomes interact with each other over a length scale of 10 nm ℓ 𝑖𝑛𝑡 Width of the smooth boundary of heterochromatic domains ∼ 50 𝑛𝑚 Observed from ChromSTEM imaging (Figure 2c) 𝜂 Penalty associated with formation of interfaces ∼ 10 -21 𝐽/𝑛𝑚 Calculated as ℓ 𝑖𝑛𝑡 = √ 𝜂Ω 𝑘 𝐵 𝑇

, (Section S1.5)

Our model predicts that chromatin exists in two phases which are not only transcriptionally distinct but are also differentially constituted. For instance, in our model, the euchromatin phase primarily consists of acetylated chromatin with a volume fraction 𝜙 𝑒 𝐸𝐶 ∼ 0.2 -0.3 and is water rich with volume fraction of nucleoplasm 𝜙 𝑛 𝐸𝐶 ∼ 0.7 -0.8. However, the euchromatin phase has no methylated heterochromatin content, i.e., 𝜙 ℎ 𝐸𝐶 ∼ 0. On the other hand, the heterochromatin phase is rich in heterochromatic content 𝜙 ℎ 𝐻𝐶 = 𝜙 ℎ 𝑚𝑎𝑥 , and has a very little water content, with 𝜙 𝑛 𝐻𝐶 ∼ 0.02.

Note that the heterochromatin phase is only water-poor, and not water-free, i.e., 𝜙 𝑛 𝐻𝐶 ≠ 0. Being water-poor, the heterochromatin phase is highly compacted, but not so compacted as to exclude water. This distinction becomes important because even the highly compacted heterochromatin phase includes multiple chromatin associated proteins such as HP1 [15] and certain histone methyltransferases [31].

The extent of chromatin compaction in the heterochromatin phase i.e., the maximum heterochromatin volume fraction 𝜙 ℎ 𝐻𝐶 = 𝜙 ℎ 𝑚𝑎𝑥 is chosen suitably to allow the presence of some euchromatin and nucleoplasm within the heterochromatic phase. We observe a water content of 𝜙 ℎ 𝐻𝐶 ∼ 0.04. However, the exact values of these parameters, or the exact density (in Mbp/𝜇m 3 ) of chromatin within the heterochromatin phase, could be obtained by molecular dynamic simulations at a resolution of nucleosomes, calibrated to single nucleosome live-cell super-resolution imaging (such as [29, 32]). Single molecule localization microscopy (SMLM), which allows visualization of individual nucleosome clusters combined with super-resolution fluctuation-assisted bindingactivated localization microscopy (fBALM) has recently emerged as a promising avenue towards extracting the nucleus-wide 3D spatial DNA density (in Mbp/𝜇m 3 ) [35]. Further, a combined computational and imaging framework involving polymer-based modeling at nucleosome level and super-resolution imaging, called Modeling immuno-OligoSTORM (MiOS) has previously been developed by Neguembor, et al. [36], which could potentially be used for observing such finescale features of DNA organization.

As an example, we modify the overall water content of the nucleus (𝜙 ̅ 𝑛 ) to modulate the extent of chromatin compaction and water content within the heterochromatin phase (𝜙 𝑛 𝐻𝐶 ). The extent of chromatin compaction can also be modified in other ways, such as varying 𝜙 ℎ 𝑚𝑎𝑥 , but here we present one example to show that our model predictions are agnostic to the exact choice of these parameters. The chromatin organization predicted by the model is shown in Figure S15. Note that in Figure S15a the distribution of water in the nucleus is shown (as opposed to that of heterochromatin). However, the red region (with low water content) is still heterochromatin, and blue region (with high water content) is euchromatin. Further, we have annotated the exact value of local volume fraction of nucleoplasm within a heterochromatin domain, i.e., 𝜙 𝑛 𝐻𝐶 . It can be seen that as the nucleus becomes water rich, the heterochromatin phase compaction varies. This is because changing the water content in the nucleus also changes the location of the energy well on the energy landscape.

For a specific case of 𝜙 ̅ 𝑛 = 0.6, we see that as the transcription increases the heterochromatin domain size scale in the interior and along the periphery reduces -in qualitative agreement with predictions reported in the main manuscript (Figure S15b). Thus, while choosing the exact level of chromatin compaction does quantitatively improve our prediction accuracy, it does not modify our model results, nor does it change the underlying physics of the model.

The model proposed in this paper incorporates three nuclear constituents -nucleoplasm, heterochromatin and euchromatin. These constituents are mixed with each other and form two stable phases -heterochromatin which is compacted and prominently methylated, and euchromatin which is loosely packed and prominently acetylated.

However, euchromatin as well as heterochromatin can have different subtypes depending on location (e.g., lysine site) and extent (e.g., mono-, di-or tri-) of methylation. These posttranslational epigenetic variations can induce different functional properties to chromatin. As an example, H3K9me3 is expected to form a core of constitutive heterochromatin, which remains compacted at all stages of development in the cell, and all cell types [37, 38]. On the other hand, H3K27 in its trimethylated form H3K27me3 is a hallmark of Polycomb facultative heterochromatin, which can be reversibly switched between expressive (H3K27me1/2) or repressive (H3K27me3) forms [37, 38]. The formation of different classes of heterochromatin (constitutive and facultative) can involve different classes of epigenetic enzymes such as methyltransferase SUV39H1 and SUV39H2 for H3K9me3 or EZH2 for H3K27me3 [37].

As a simplified model to capture multiple states of chromatin, let us consider three nuclear constituents -euchromatin (which as a simplification is nucleoplasm rich by definition), constitutive heterochromatin and Polycomb-marked facultative heterochromatin. We consider that the three constituents can mix with different volume fractions to stably form three phaseseach phase rich in one of each constituent. The free-energy landscape of this triphasic system can be defined as (analogous to Eq S1) [39], Euchromatin well F. Heterochromatin well C. Heterochromatin well ( olume fraction of euchromatin) ENERGETI S REA TION INETI S ONSTITUENTS 1 olume Fraction of euchromatin 2 olume Fraction of facultative heterochromatin 3 olume Fraction of constitutive heterochromatin 1 2 3 23 32 Figure S16: The key ingredients of a ternary phase field model to capture three stable phases of chromatin -euchromatin, facultative and constitutive heterochromatin.