Bifractal nature of chromosome contact maps
Simone Pigolotti, Mogens H. Jensen, Yinxiu Zhan, Guido Tiana
CChromosome contact maps are bifractal
S. Pigolotti , M. H. Jensen , Y. Zhan , , and G. Tiana Biological Complexity Unit, Okinawa Institute of Science and Technology Graduate University,Onna, Okinawa 904-0495, Japan. The Niels Bohr Institute, Blegdamsvej 16, 2100 Copenhagen Ø, Denmark. Friedrich Miescher Institute for Biomedical Research, Maulbeerstrasse 66, 4058 Basel, Switzer-land. University of Basel, Petersplatz 1, 4001 Basel, Switzerland Department of Physics and Center for Complexity and Biosystems, Universit`a degli Studi di Mi-lano and INFN, via Celoria 16, 20133 Milano, Italy.
Modern biological techniques such as Hi–C permit to measure probabilities that differentchromosomal regions are close in space. These probabilities can be visualised as matricescalled contact maps. Contact maps are characterized by self-similar blocks along the di-agonal, corresponding to a hierarchy of tightly packed chromosomal domains. Statisticalproperties of these maps, and thus of chromosomal interactions, are usually summarized bythe scaling of the contact probability as a function of the genomic distance. This approachis fraught with limitations, since such scaling is usually rather noisy. Here, we introduce amultifractal analysis of chromosomal contact maps. Our analysis reveals that Hi–C maps arebifractal, i.e. they are complex geometrical objects characterized by two distinct fractal di-mensions. To rationalize this observation, we introduce a model that describes chromosomes a r X i v : . [ q - b i o . B M ] M a r s a hierarchical set of domains nested in each other and we solve it exactly. The predictedmultifractal spectrum is in excellent quantitative agreement with experimental data. Ourtheory yields to a more robust estimation of the scaling exponent of the contact probabil-ity than existing methods. By applying this method to experimental data, we detect subtleconformational changes among chromosomes during differentiation of human stem cells. During cellular interphase, mammalian chromosomes assume a globular structure in the nu-cleus
1, 2 . Their conformational properties can be studied in vivo with a set of techniques calledchromosome conformation capture (3C) , most notably their genome-wide version called Hi–C .Results of Hi–C experiments can be represented as matrices whose elements are proportional tothe probability that two chromosome regions are in contact in space . Hi–C measurements havepaved the way for a mechanistic understanding of chromosome folding. In particular, they haverevealed that mammalian chromosomes are characterized by a hierarchy of nested, tightly con-nected structures . Structures below the megabase scale are usually called topological associatingdomains (TADs). They are particularly relevant for gene regulation
6, 8 , although they do not seemto be privileged over other levels in the hierarchy from a structural point of view . It has beensuggested that the hierarchical structure of Hi–C matrices can be studied by means of an iterativecoarse-graining procedure .Information contained in Hi–C matrices is often characterized by the decay of the contactprobability p ( (cid:96) ) of pairs of chromosomal regions i and j with respect to their genomic distance2 ≡ | i − j | . Such decay appears to follow a power law p ( (cid:96) ) ∼ (cid:96) − β . (1)The contact probability exponent β is often estimated to be slightly smaller than one
4, 11, 12, 26 . Suchlow values are incompatible with simple equilibrium homopolymeric models . In contrast, non-equilibrium models such as the crumpled globule
4, 13 are able to account for such low exponents.Other proposed mechanisms include mediation of polymer interaction by other molecules , activeloop-extrusion , and finite–size effects in heteropolymers . In any case, the apparent power-lawrange of the contact probability is usually limited to one decade or less, see e.g.
4, 11, 12 . Therefore,estimates of the exponent β are rather sensitive to the choice of the fitting range.In this paper, we robustly characterize scale-invariant properties of chromosomes using thetheory of multifractals. This theory has been developed to characterize heterogeneous systemscharacterized by scale invariance . This analysis reveals that chromosome contact maps are bifractal , i.e. multifractal objects characterized by two distinct fractal dimensions, a behaviorpreviously reported in studies of surface roughness , distribution of matter in the universe andturbulence .We introduce our theory using a Hi–C map of chromosome 1 in mouse embryonic stem cells,see Fig. 1 a . The contact probability seems to decay as a power law, at least for relatively shortgenomic distances, see Fig. 1 b . However, the local logarithmic slope of the contact probabilitydoes not present the clear plateau characteristic of a power law behavior, see Fig. 1 c .3o characterize scaling properties of chromosomes in a more rigorous way, we study Hi–Cmaps as multifractals. A multifractal is a system described in terms of a density, that in our case isthe density of number of counts in the Hi–C map. We study the scaling property of this density byconstructing two-dimensional histograms of the Hi–C map on grids of different linear resolution (cid:15) ,see Fig. 1 d . Geometrical structures of linear size smaller than (cid:15) are not resolved in these maps. Thesmallest possible value of (cid:15) is the resolution of the original Hi–C map, in our case (cid:15) = 4 · bp .We define the probability p ij ( (cid:15) ) in bin at coordinates i , j at resolution (cid:15) . We always work withnormalized maps, so that (cid:80) ij p ij ( (cid:15) ) = 1 for all choices of (cid:15) .We assume that the system is scale invariant, at least for (cid:15) small enough. This means that theprobability in a bin scales as a power law of the resolution: p ij ( (cid:15) ) ∼ (cid:15) α , (2)where ∼ denotes the leading behavior and α is the scaling exponent associated to the density. Sincethe map is not homogeneous, different bins can be in principle characterized by different values of α . The number of bins N ( α ) associated to a given value of α must also scales as a power of (cid:15)N ( α ) ∼ ρ ( α ) (cid:15) − D ( α ) . (3)The exponent D ( α ) characterizes the scaling of the number of bins of linear size (cid:15) necessary tocover the set with density exponent α and therefore can be interpreted as the fractal dimensionassociated to this set. The quantity ρ ( α ) is a prefactor independent of (cid:15) .Computing D ( α ) directly is often unpractical. A convenient related quantity is the partition4unction Z ( q, (cid:15) ) , defined by
19, 20 Z ( q, (cid:15) ) = (cid:88) ij [ p ij ( (cid:15) )] q . (4)The name “partition function” originates from an analogy with equilibrium statistical mechanics,where the exponent q plays the role of an inverse temperature. Indeed, in an inhomogeneoussystem, for q → (high temperature), all bins give similar contributions to the sum in Eq. (4),whereas for large q (low temperature) the sum is dominated by relatively few terms characterizedby largest values of the measure p ij . We link the partition function with the fractal dimensions bycollecting in Eq. (4) all terms with the same value of α : Z ( q, (cid:15) ) ∼ (cid:88) α ρ ( α ) (cid:15) qα − D ( α ) . (5)In the limit of small (cid:15) , we evaluate this sum with the saddle point method. This calculation showsthat, for a scale-invariant system, one also expects a power-law scaling of the partition function asa function of the resolution Z ( q, (cid:15) ) ∼ (cid:15) K ( q ) (cid:15) → . (6)The function K ( q ) is the multifractal spectrum and characterizes the scaling behavior with re-spect to the resolution of the partition function . The saddle point evaluation of Eq. (5) leads toconclude that the dimensions D ( α ) and the multifractal spectrum K ( q ) are related by a Legendretransform: K ( q ) = max α [ qα − D ( α )] . (7)A linear multifractal spectrum indicates that the system is homogeneous, i.e. all of its parts arecharacterized by the same fractal dimension α = D and thus by the same scaling behavior N ∼ − D and p ij ∼ (cid:15) D . Conversely, a non-linear multifractal spectrum signals a diversity of scalingexponents and associated dimensions.The partition function Z ( q, (cid:15) ) computed from Hi–C data presents a clear power law behavioras a function of (cid:15) , as predicted by Eq. (6), see Fig. 1 e . In this case, for (cid:15) slightly larger than itsminimum value, the local logarithmic slope of the partition function is essentially flat, see Fig. 1 f .This means that the power–law behaviour of Z ( q, (cid:15) ) is much more robust than that of p ( (cid:96) ) .The multifractal spectrum K ( q ) of chromosome 1 presents two different linear regimes, oneat low moments with the properties of a two-dimensional object and one at high moments charac-terized by a fractal dimension between one and two, see Fig. 1 g . A system with such propertiesis called a bifractal
21, 22, 22 . Since the exponent q is analogous to an inverse temperature, the sharpchange of slope in the spectrum of a bifractal system is akin to a phase transition in equilibriumstatistical physics .To rationalize these observations, we introduce a hierarchical domain model of Hi–C maps.We define the model in terms of an iterative transformation of a measure on the unit square [0 , × [0 , . At the first iteration, the measure is given by a × matrix with diagonal element a andoff–diagonal elements b , with a > b > . At each following iteration n we generate a n × n matrix. The off-diagonal element of the matrix at the previous iteration becomes a × block ofidentical values in the matrix at the n th iteration. The diagonal blocks are further multiplied bythe original matrix (see Fig. S1 in the Supp. Mat). We impose b = 1 / − a , so that the measureremains normalized at each iteration. This means that, effectively, the model is defined in terms6f a single free parameter a . Because of the normalization and the requirement that the measureshould be concentrated along the diagonal, such parameter falls in the range / ≤ a ≤ / .The parameter a represents the weight of domains compared to the rest of the Hi–C map.Matrices with larger a have a measure more concentrated along the diagonal, whereas matriceswith lower a are characterized by a more uniform measure, see Fig. 2 a . In the limiting case a = 1 / the matrix is uniform at all iterations, whereas for a = 1 / the matrix is characterized bya uniform measure on the diagonal, whereas all the off-diagonal elements are zero.The partition function of the hierarchical domain model can be expressed as Z ( q, n ) = n (cid:88) k =0 exp[ ξ ( k )] , (8)where we defined the exponent ξ ( k ) = [ n + max( n − k − , kq ln a ++ (1 − δ kn ) q ln b − max( n − k − , q ln 4 . (9)A saddle–point estimation of the partition function gives dξdk = − ln 2 + q ln 4 a, (10)which is positive for q > q c = (ln 2) / (ln 4 a ) . This means that, for q ≥ q c , the leading term is either ξ ( n ) or ξ ( n − . Since ξ ( n ) − ξ ( n −
1) = q ln( a/b ) , the maximum of the exponent is attainedat k = n . Thus, for large n , the partition function scales as Z ( q, n ) ∼ (cid:15) K ( q ) with the length scale (cid:15) = 2 − n and the multifractal spectrum K ( q ) = − q ln a ln 2 − . (11)7or a = 1 / , the matrix p ij ( (cid:15) ) is diagonal at each iteration. In this case, Eq. (11) predicts alinear spectrum with slope D = 1 , consistent with the fact that the geometric set is equivalent to aone dimensional line. For a = 1 / the distribution is uniform on the square, and Eq. (11) correctlyreturns D = 2 . Between these two limiting cases, Eq. (11) predicts a fractal distribution with adimension D = − ln a/ ln 2 between and .We now focus on the “high temperature phase” where − ln 2 + q ln 4 a < . In this case, theterm dominating the scaling is k = 0 , so that K ( q ) = 2( q − . (12)Since in the high temperature phase the scaling is determined by the terms far from the diagonal,the spectrum is that of a regular two-dimensional set.Summarizing, our calculation predicts a multifractal spectrum characterized by two linearregimes: one with slope − log( a ) / log(2) for q > q c = (ln 2) / (ln 4 a ) , Eq. (11), and one with slope for q < q c , Eq. (12). Such predictions are in excellent agreement with numerical simulations, asshown in Fig. 2 a , with very small discrepancies for high values of q and low values of a arisingfrom finite size effects. These results confirm the validity of our saddle point approximation.Strikingly, our theory predicted extremely well also the multifractal spectra of real chromo-somes , with a value of a ≈ . and very little variability among the first three chromosomes ofmouse embryonic stem cells (mean squared displacement MSD ≈ . in all three cases), seeFig. 2 c . To verify whether the multifractal spectrum in Fig.2 is truly a signature of a hierarchical8echanism, we computed numerically the spectra of two null models. In both null models, thecontact probability decays as a power law of the distance, to which we add domains at a singlelength scale M in the first null model and random noise in the second. Additional details on thenull models are presented as Supplementary Information. The solution of Eqs. (11) and (12) pro-vides a poor fit to this first null model (MSD in the range . ∼ . depending on parameters,see Supplementary Fig. S2). The second null model provides a better fit, with MSD ≈ . for small noise and and β ≈ , see Supplementary Fig. S2. This suggests that a pure power–lawbehavior with the addition of noise provides an alternative explanation for the observed scalingbehavior.Our theory accurately describes also Hi–C maps of Drosophila chromosomes (Supplemen-tary Fig. S3, average MSD ≈ . for the first four chromosomes) and of human chromosomes(Supplementary Fig. S4, average MSD ≈ . for all chromosomes except chromosome Y).These observations support that the bifractal spectrum predicted by the hierarchical domain modelis compatible with that observed in a broad class of higher organismsFrom the model, we also predict the scaling of the contact probability with the genomicdistance. In our framework, the contact probability P ( (cid:96) ; n ) at the n -th iteration can be expressedby summing the probabilities of blocks at a genomic distance (cid:96) from the diagonal: P ( (cid:96) ; n ) = (cid:88) ij p ij ( n ) δ | i − j | ,(cid:96) . (13)From this expression, we derive that the decay exponent of the contact probability with the genomic9istance is β = ln(4 a )ln(2) . (14)Derivation of Eq. (14) is presented in SI. For / ≤ a ≤ , Eq. (14) predicts contact probabilityexponents in the range β ∈ [0 , . This prediction is supported by numerical simulations, see Fig.S3 in the Supp. Mat. This means that the hierarchical structure of contact matrices automaticallyleads to exponents β < , as observed in chromosomes.We test more extensively our method using datasets from different Hi–C experiments usingmouse embryonic stem cells (mESC). We fit the multifractal spectrum of each dataset and obtainthe corresponding exponent β via Eq. (14). We first test the robustness of the multifractal analysisacross two Hi–C experiments in mESC from two different labs
17, 30 , both following the Hi–Cprotocol in solution , see Fig. 3 a . Both the multifractal and the power law methods predict that thevariability of the exponent β across chromosomes is significantly larger than the variability of theexponent across replicates. To figure out which method performs best, we use bootstrapping totest against the null hypothesis that the values of β of different chromosomes are paired at random.The multifractal method permits to exclude random pairing of chromosomes (p–value . · − )much better than the direct fit (p–value . , see Supplementary Table S1). Moreover, the meansquared error of β between the two replicates is smaller in the multi–fractal case compared to directfit ( . vs . ). The χ associated with the direct fit is affected by a strong systematic error,although remaining quite correlated (see Supplementary Fig. S6). This effect is much milder inthe multifractal approach. The multifractal and direct fit methods are similarly robust with respectto varying the resolution of Hi–C maps (see Supplementary Fig. S7).10e apply the two methods to attempt to distinguish among experiments on the same celllines, but following different experimental protocols. These different protocols mainly differ inthe ligation step of digested DNA fragments. In the original “in solution” protocol, the ligation isperformed in a diluted solution . In other protocols, the ligation step is either carried out in intactnuclei (“in situ” protocol ) or on solid substrates (“TCC” protocol ). The two latter protocolsare able to produce Hi-C maps with better signal-to-noise ratio compared to the in solution protocol
32, 36 . Both the multifractal and the direct fit methods show that the values of β obtained from in situand TCC experiments are markedly different from those obtained in solution (see SupplementaryFig. S6). The multifractal method estimates compatible scaling exponents for in situ and TCCprotocols (the p–value associated with the null hypothesis that chromosomes are paired at randomis . , see Supplementary Table S2), while a direct fit does not permit to draw this conclusion (p–value . ). This result suggests that the in situ and TCC protocols are able to produce statisticallycompatible Hi–C maps due to their high signal-to-noise ratio.We compare the two methods in detecting differences between wild–type and mutant celllines, in which either the Nipbl gene or the CTCF gene is knocked down (see Fig. 3 b ).Knock–down of these genes has been shown to disrupt chromosome folding. In particular, Nipblknock–down leads to loss of TAD structures and global changes in scaling properties . We quan-tify average differences in exponents between the wild–type and the mutant in terms of χ of pairsof chromosomes, weighted by their mean squared errors (see Supplementary Table S2 and Sup-plementary Fig. S7). Multifractal analysis detects more marked differences ( χ = 214 ) comparedto the direct fit ( χ = 4 . ), although both methods highlight a statistically significant difference11etween the two sets (p–values . · − for the direct fit and < − for the multifractal). Knock–down of CTCF also causes a loss of TAD structure, but without a clear effect on genome–wide scal-ing . Nonetheless, the values of β of different chromosomes obtained by multifractal analysisreveal a large ( χ = 86 ) and significant (p–value < − ) difference between the scaling prop-erties of wild–type and CTCF–deficient cells. Also in this case, the direct fit detects less markeddifferences ( χ = 25 , see Supplementary Table S2).We apply the multifractal analysis to elucidate how chromosome structure changes duringcellular differentiation. To this aim, we analyze Hi–C data obtained from different human cell linesat different stages of early differentiation . In temporal order, the cell lines are: trophectoderm,embyonic stem cells (ESC), mesoderm, neuronal precursor cells (NPC) and mesenchymal cells.Mesoderm and NPC can be regarded as in similar development stages. The value of β obtained bymultifractal analysis (see Fig. 3 c ) tends to increase upon differentiation(p–value of Kendall’s tautest applied to all chromosome pairs < − ). A similarly clear conclusion cannot be obtainedfrom the values of β estimated by direct fit (see Fig. 3 d , p–value . ). The increase of β pointsto the idea that more differentiated cell types require more insulated domains to achieve morespecialised function.Multifractal analysis of Hi–C data is a powerful statistical tool to characterize scaling prop-erties of chromosomes. We have shown that contact maps of chromosomes are bifractal, i.e. theyare characterized by two distinct fractal dimensions. To understand this observation, we proposeda hierarchical domain model, whose analytical solution is in striking quantitative agreement with12bserved multifractal spectra. Our theory implicates a power-law scaling of the contact probabilitywith the genomic distance with exponents β lower than one, also in agreement with experimentalobservations. The multifractal method is sensitive enough to discard a null model with power-law decay of the contact probability, but domains on a single length scale only. The predictedform of the multifractal spectrum provides a more stringent prediction than the contact probabilityexponent alone. The analysis proposed here constitutes an excellent benchmark to select amongdifferent polymer models that provide similar values of β .Our results indicate that scaling properties of chromosomes are a direct consequence of thehierarchical structure of chromosome domains
9, 27 . Recent work has suggested that such domainstructure is generated by a ”hierarchical folding” mechanism, mediated by different proteins suchas cohesin and CTCF . The precise mechanism driving the folding of chromosomes has beensubject of debate . It will be interesting to test whether the activity of these factors can produceself-similar structures compatible with our observations. The shape of the multifractal spectrumis controlled by a single parameter, that also controls the contact probability exponent β via Eq.(14). The determination of β from the fit of the multifractal spectrum is a much more robust wayof characterizing the Hi–C map than the direct fit of β . Multifractal analysis is importantly alsomore sensitive in highlighting subtle structural differences as shown in the case of Nipbl and CTCFknock–down. The analysis of the multifractal spectrum is simple and robust enough to become aroutinely tool to analyze contact maps, both from polymer models and experiments, and capturesubtle differences among them. 13. L. Giorgetti, R. Galupa, E. P. Nora, T. Piolot, F. Lam et al. , Predictive polymer modelingreveals coupled fluctuations in chromosome conformation and transcription.
Cell , 950(2014).2. G. Tiana, A. Amitai, T. Pollex, T. Piolot, D. Holcman et al. Structural Fluctuations of theChromatin Fiber within Topologically Associating Domains.
Biophys J. , 1234 (2016)3. J. Dekker, K. Rippe, M. Dekker and N. Kleckner,
Capturing chromosome conformation.
Sci-ence, , 1306 (2002).4. E. Lieberman-Aiden, N. L. van Berkum, L. Williams, M. Imakaev, T. Ragoczy, A. Telling et al. , Comprehensive mapping of long-range interactions reveals folding principles of thehuman genome.
Science , 289 (2009).5. J. Redolfi, Y. Zhan, C. Valdes, M. Kryzhanovska, I. Misteli Guerreiro et al. , DamC revealsprinciples of chromatin folding in vivo without crosslinking and ligation.
Nature Str. Mol.Biol. , 471 (2019)6. Q. Szabo, et al Principles of genome folding into topologically associating domains , ScienceAdvances, , (2019)7. J. Nuebler, G. Fudenberg, M. Imakaev, N Abdennur and L. A. Mirny, Chromatin Organizationby an Interplay of Loop Extrusion and Compartmental Segregation.
Proc. Natl. Acad. Sci.USA , E6697 (2018). 14. D. G. Lupi´a˜nez, K. Kraft, V. Heinrich, P. Krawitz, F. Brancati, E. Klopocki et al. Disruptionsof topological chromatin domains cause pathogenic rewiring of gene-enhancer interactions.
Cell , 1012 (2015).9. Y. Zhan, L. Mariani, I Barozzi, E. G. Schulz, N. Bl¨uthgen et al. , Reciprocal insulation analysisof Hi-C data shows that TADs represent a functionally but not structurally privileged scale inthe hierarchical folding of chromosomes.
Genome Research , 479 (2017).10. A. M. Chiariello, S. Bianco, C. Annunziatella, A. Esposito and M. Nicodemi, The scaling fea-tures of the 3D organization of chromosomes are highlighted by a transformation la Kadanoffof Hi-C data.
Europhys. Lett. et al. Chromatin extrusion explainskey features of loop and domain formation in wild-type and engineered genomes.
Proc. Natl.Acad. Sci. USA , E6456 (2017).12. B. Bonev, N. M. Cohen Q. Szabo, L. Fritsch, G. L. Papadopoulos et al. , Multiscale 3D GenomeRewiring during Mouse Neural Development.
Cell , 557.e1 (2017).13. L. Mirny,
The fractal globule as a model of chromatin architecture in the cell.
Chrom. Res. ,37 (2011).14. M. Barbieri, M. Chotalia, J. Fraser, L. M. Lavitas, J. Dostie, A. Pombo and M. Nicodemi, Complexity of chromatin folding is captured by the strings and binders switch model.
Proc.natl. Acad. Sci. USA , 16173 (2012). 155. G. Fudenberg, M. Imakaev, C. Lu, A. Goloborodko, N. Abdennur and L. Mirny,
Formation ofChromosomal Domains by Loop Extrusion.
Cell Rep. , 2038 (2016).16. Y. Zhan, L. Giorgetti and G. Tiana, Looping probability of random heteropolymers helps tounderstand the scaling properties of biopolymers.
Phys. Rev. E , 032402 (2016).17. L. Giorgetti, B. R. Lajoie, A. C. Carter, M. Attia, et al. , Structural organization of the inactiveX chromosome in the mouse.
Nature , 575 (2016).18. G. Parisi and U. Frisch in ”Turbulence and Predictability in Geophysical Fluid Dynamics”, p.8487, eds. M. Ghil et al (1985).19. R. Benzi, G. Paladin, G. Parisi, and A. Vulpiani,
On the multifractal nature of fully developedturbulence and chaotic systems.
Journal of Physics A , 3521 (1984).20. T.C. Halsey, M.H. Jensen, L.P. Kadanoff, I. Procaccia, and B. Shraiman, Fractal measures andtheir singularities: The characterization of strange sets.
Phys. Rev. A , 1141 (1986).21. Bhushan B, Majumdar A. Elastic-plastic contact model for bifractal surfaces.
Wear, (1),53-64. (1992).22. Balian, R., and R. Schaeffer.
Scale-invariant matter distribution in the universe. ii-bifractalbehaviour.
Astronomy and Astrophysics , 373-414 (1989).23. Iyer, K.P., Sreenivasan, K.R. and Yeung, P.K.,
Circulation in high Reynolds number isotropicturbulence is a bifractal. arXiv preprint arXiv:1902.07326 (2019).164. T. Bohr and M.H. Jensen,
Order parameter, symmetry breaking, and phase transitions in thedescription of multifractal sets.
Phys. Rev. A , 4904 (1987).25. G. Tiana and L. Giorgetti, Integrating experiment, theory and simulation to determine thestructure and dynamics of mammalian chromosomes.
Curr. Opin. Struct. Biol. , 11 (2018).26. Y. Zhan, L. Giorgetti and G. Tiana, Modelling genome-wide topological associating domainsin mouse embryonic stem cells.
Chromosome Res. et al , A 3D map of the human genome at kilobaseresolution reveals principles of chromatin looping.
Cell,
Organization and function of the 3D genome.
Nat. Rev. Gen. , 661(2016).29. M. J. Rowley and V. G. Corces, Organizational principles of 3D genome architecture.
Nat.Rev. Gen. , 789 (2018).30. E. P. Nora, A. Goloborodko, A. L. Valton, J. H. Gibcus, A. Uebersohn, N. Abdennur et al.Targeted Degradation of CTCF Decouples Local Insulation of Chromosome Domains fromGenomic Compartmentalization. Cell, , 930 (2017).31. W. Schwarzer, N. Abdennur, A. Goloborodko, A., Pekowska, G. Fudenberg, Y. Loe-Mie et alTwo independent modes of chromatin organization revealed by cohesin removal.
Nature, ,1270 (2017). 172. R. Kalhor, H. Tjong, N. Jayathilaka, F. Alber and L. Chen,
Genome architectures revealed bytethered chromosome conformation capture and population-based modeling.
Nature Biotech-nology , 90 (2012).33. J. Redolfi, Y. Zhan, C. Valdes-Quezada, et al DamC reveals principles of chromatin folding invivo without crosslinking and ligation. Nat. Struct. Mol. Biol, , 471-480, (2019).34. E. Lieberman-Aiden , et al Comprehensive mapping of long-range interactions reveals foldingprinciples of the human genome. Science, , 5950, (2019).35. J. R. Dixon, et al Chromatin architecture reorganization during stem cell differentiation.
Na-ture, , 331336, (2015).36. T. Nagano, et al Comparison of Hi-C results using in-solution versus in-nucleus ligation.
Genome Biology, , (2015). Acknowledgements
We thank Luca Giorgetti for help with the data.
Competing Interests
The authors declare that they have no competing financial interests.
Correspondence
Correspondence and requests for materials should be addressed to S.P. (email: [email protected]). a , Hi–C map of chromosome 1 of mouse embryonic stem cells (mESC) . Darker shadesof red correspond to higher values of the contact probability p ij b , Scaling of the average contactprobability p ( (cid:96) ) as a function of the distance (cid:96) = | i − j | . The continuous line is a best-fit of a powerlaw, Eq. (1), in the range of distances [10 , ] , yielding β ≈ . . c , Local logarithmic slope d ln P ( (cid:96) ) / d(ln (cid:96) ) of the contact probability. d , Coarse-graining of the Hi–C map at increasingvalues of the resolution (cid:15) = 10 , , bp . e , Scaling of Z ( q ) as a function of resolution, seeEqs. (4) and (6). Different lines correspond to different values of q , increasing from top to bottomfrom q = 0 to q = 5 at intervals of ∆ q = 0 . . f , Local logarithmic slopes of the first momentaof Z ( q ) . g , Corresponding multifractal spectrum K ( q ) . Here and throughout the paper, spectraare obtained by a fit in log-log scale of Z ( q, (cid:15) ) versus (cid:15) in the range (cid:15) ∈ [10 , ] unless specifiedotherwise. 19 n = n = a = 0.3 a n = a = 0.4 -2 0 2 4 6 8 b K ( q ) qHierarchical Domain Modela = 0.3a = 0.35a = 0.4a = 0.45-2 0 2 4 c K ( q ) qExperimental Datachr 1chr 2chr 3 Figure 2: a , Construction of the hierarchical domain model for two different values of the parame-ter a . b , Colored symbols are the multifractal spectrum calculated by the hierarchical domain modeat different values of a , the solid curve is the prediction of Eqs. (11) and (12). c , The multifractalspectra of the first three chromosomes of mESC (colored symbols) are barely distinguishable fromeach other. The solid curve is a fit of Eqs. (11) and (12), giving a ≈ . .20 .000.800.60 β multifractal direct fi t wt wt wt wtmut mut mut mut Nora et al 2017Nora et al 2017Schwarzer et al2017 Schwarzer et al 2017 dataset β multifractal direct fi t in solution in solution type of experiment r1r1 r1r2 r2 a in situ TCC in situTCC b trophectodermESC mesodermNPC mesenchymal0.600.700.800.901.00 β trophectodermESC mesodermNPC mesenchymal0.600.70 β c d multifractal direct fi t r1r2 r1r2 Figure 3: a Exponents β calculated from the multifractal spectrum using Eq. (14) and from adirect fit of the data (as in Fig. 1B) with the alternative Hi–C methods. One replicate of in situ , TCC and two independent replicates of in solution
17, 30
Hi–C are shown in situ
32, 33 and insolution
17, 30 . The black bar indicates the median over the chromosomes, the colored box indicatesthe first and third quartiles and black dots are outliers. b Value of β obtained for the wild–type (wt)system and for a mutant (mut) in which either CTCF or Nipbl is depleted in experiments. c Exponents β obtained with the multifractal method for cell lines at different stages of human stemcell development . dd