Introduction {#Sec1}
============

Through remarkable Hi-C experiments^[@CR1]--[@CR6]^, based on the Chromosome Conformation Capture (3C) technique^[@CR7]^, indirect glimpses of how the genome in a number of species is organized is starting to emerge. Because chromosome lengths are extremely large, ranging from tens of million base pairs in yeast to billion base pairs in human cells, they have to fold into highly compact structures in order to be accommodated in the cell nucleus. This requires that loci that are well separated along the one-dimensional genome sequence be close in three-dimensional (3D) space, which is made possible by forming a large number of loops. The high-throughput Hi-C technique and its variants are used to infer the probability of genome-wide contact formation between loci. In order to determine the contact probabilities between various loci in a genome, Hi-C experiments are performed in an ensemble of millions of cells. The readouts of the Hi-C experiment are contact frequencies between a large number of loci from instantaneous snapshots of each cell, which are then used to construct the contact maps (Hi-C maps). The contact map is a matrix (2D representation) in which the elements represent the probability of contact between two loci that are separated by a specified genomic distance. A high contact count between two loci means that they interact with each other more frequently compared with ones with low contact count.

A complementary and potentially a more direct way to determine genome organization is to measure spatial distances between loci using a low-throughput fluorescence in situ hybridization (FISH) technique^[@CR8],[@CR9]^. In addition to providing 3D distances in fixed cells, recently developed CRISPR--dCas9 FISH can be used to assay the dynamic behavior of loci in real time^[@CR10]--[@CR12]^. However, due to the current limitation of the number of distinct color probes, this method provides distance distribution information for only a small number of loci.

FISH and Hi-C, which are entirely different experimental techniques, provide data on different aspects of genome organization. As noted in recent reviews^[@CR13],[@CR14]^, there are problems associated with each method. It is difficult to reconcile Hi-C and FISH data for the following reasons. In interpreting the Hi-C contact map, one makes the intuitive assumption that loci with high probability contact must also be spatially close. However, it has been demonstrated using Hi-C and FISH data on the same chromosome that high contact frequency does not always imply proximity in space^[@CR13],[@CR15]--[@CR17]^. It should be noted that in most cases, the Hi-C and FISH measurements agree very well^[@CR8],[@CR9],[@CR18],[@CR19]^. However, from a purely theoretical perspective, even a single contradiction is intriguing if the experimental errors can be ruled out. An outcome of our theory is that the discordance between FISH and Hi-C data arises because of extensive heterogeneity, which is embodied by the presence of a variety of conformations adopted by chromosomes in each cell. There are a variety of reasons, including differing fixation conditions and presence of two or more subpopulation of cells in which the chromosomes are present in distinct conformations, which could give rise to the discordance between FISH and Hi-C data, as lucidly described recently^[@CR13],[@CR14]^. Contact between two loci could be a rare event, not present in all cells, which is captured in a Hi-C experiment by performing an ensemble average. We show using a precisely solvable model that due to the absence of a contact between two specific loci in a number of cells, those with higher contact frequency could be spatially farther on an average than two others with lower contact frequency. In contrast, the probability of contact formation using the FISH method can only be obtained if the tail (small distance) of the distance distribution between locus *i* and *j* can be accurately measured. For a variety of reasons, including the size of the probe and the signal strength, this not altogether straightforward using the FISH technique. Thus, in order to combine the data from the two powerful techniques, it is crucial to establish a theoretical basis with potential a practical link, between the contact probability and average spatial distance.

Setting aside the conditions under which FISH and Hi-C are performed (see recommendations for comparing the results from the two techniques with minimum bias which are described elsewhere^[@CR13]^) insights into the discordance between the two methods, when they occur, can be obtained using polymer physics concepts. Recently, Fudenberg and Imakaev^[@CR15]^ performed polymer simulations using a strong attractive energy between two labeled loci and a tenfold weaker interaction between two other loci that are separated by a similar genomic distance. In addition, they also reported simulations based on the loop extrusion model. Both these types of simulations showed there could be discordance between FISH and Hi-C, which we refer to as the FISH--Hi-C paradox. However, they did not provide any solution to the paradox, which is the principal goal of this work.

In addition, recent single-cell Hi-C^[@CR20]--[@CR22]^ and FISH experiments^[@CR8],[@CR9],[@CR18],[@CR19]^ have revealed that there are substantial cell-to-cell variations on genome organization. However, how to utilize the data reported in these experiments to enhance our understanding of 3D genome structural heterogeneity has not been unexplored. One approach is to create an appropriate polymer model based on Hi-C and imaging data, which would readily allow us to probe the structural variability using simulations^[@CR23]--[@CR26]^. Indeed, it has been shown, using Hi-C and FISH data as well as simulations^[@CR26]^, that if the conformation of the chromatin fiber is taken to be homogeneous then trends observed in the FISH data could not be predicted. However, using simulations and including two levels of chromatin organization (open and compact) qualitative trends observed in the FISH data could be recovered^[@CR26]^.

Here, we first establish a relationship between the contact probability and the mean spatial distance using an analytically solvable Generalized Rouse Chromosome Model (GRMC), which incorporates the presence of CTCF/cohein-mediated loops. The GRMC may be thought of as an ideal chromosome model, very much in the spirit of the Rouse model for polymers, in which conceptual issues such as the origin of the FISH--Hi-C paradox can be rigorously established. We first consider the solvable homogeneous limit, in which contacts are present in all the cells. In this case, precise numerical and analytical results show that there is a simple relation between the contact probability, *P*, and the ensemble mean 3D distance 〈*R*〉. However, the unavoidable heterogeneity in the cell populations in Hi-C experiments results in contacts between loci only in a fraction of cells. We first show that a direct consequence of the heterogeneity in both GRMC and chromosomes is that two loci (*m* and *n*) that have higher probability (*P*~*mn*~) of being in contact relative to another two loci (*k* and *l*) does not imply a direct spatial correlation, a finding that has already been qualitatively established in previous studies^[@CR13],[@CR15]^. In other words, the average spatial distance between *m* and *n* (〈*R*~*mn*~〉) could be larger than 〈*R*~*kl*~〉, the distance between loci *k* and *l*, even if *P*~*mn*~ \> *P*~*kl*~. These results provide a basis for understanding the origin of the FISH--Hi-C paradox.

We develop a fully theoretical approach, which allows us to provide quantitative insights into the extent of heterogeneity in chromosome organization. From our theory, it follows that the resolution of the FISH--Hi-C paradox requires invoking the notion of heterogeneity, which implies multiple populations of chromosomes coexist. By using the concepts that emerge from the study of the GRMC, we demonstrate that the information of cell subpopulations can be extracted by fitting the experimental FISH data using our theory, thus allowing us to calculate the Hi-C contact probabilities from the theoretically calculated cumulative distribution function of spatial distance (CDF)---a quantity that can be measured using FISH and super-resolution imaging methods. Our approach provides a theoretically based method to combine the available FISH and Hi-C data to produce a more refined characterization of the heterogeneous chromosome organization than is possible by using data from just one of the techniques. In other words, sparse data from both the experimental methods can be simultaneously harnessed to predict the 3D organization of chromosomes.

Results {#Sec2}
=======

Relating contact probability to mean spatial distance {#Sec3}
-----------------------------------------------------

The exact relationship between *P*~*mn*~ (contact probability between *m*th and *n*th locus) and the corresponding mean spatial distance, 〈*R*~*mn*~〉 for GRMC (see the Methods section for details of the derivation) is,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$P_{mn} = {\mathrm{erf}}\left( {\frac{{2r_{\mathrm{c}}}}{{\sqrt \pi \langle R_{mn}\rangle }}} \right) - \frac{4}{\pi }\frac{{r_{\mathrm{c}}}}{{\langle R_{mn}\rangle }}e^{ - \frac{{4r_{\mathrm{c}}^2}}{{\pi \left\langle {R_{mn}} \right\rangle ^2}}} \equiv R_0(\langle R_{mn}\rangle ).$$\end{document}$$The inverse of *R*~0~(〈*R*~*mn*~〉), the solution to Eq. ([1](#Equ1){ref-type=""}), gives the mean spatial distance 〈*R*~*mn*~〉 as a function of the contact probability *P*~*mn*~. Note that *m* and *n* are arbitrary locations of any two loci, and thus Eq. ([1](#Equ1){ref-type=""}) is general for any pair of loci.

A couple of conclusions, relevant to the application to the chromosomes, follow from Eq. ([1](#Equ1){ref-type=""}). (i) Note that Eq. ([1](#Equ1){ref-type=""}) is an exact one-to-one relation between the mean distance 〈*R*~*mn*~〉 and the contact probability *P*~*mn*~ provided *r*~c~ is known, and if the contacts are present in all the cells, which is not the case in experiments. For small *P*~*mn*~, it is easy to show from Eq. ([1](#Equ1){ref-type=""}) that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle R_{mn}\rangle \approx r_{\mathrm{c}}P_{mn}^{ - 1/3}$$\end{document}$. For the ideal GRMC, this implies that for any *m*, *n*, *k*, *l*, if *P*~*mn*~ \< *P*~*kl*~ then 〈*R*~*mn*~〉 \> 〈*R*~*kl*~〉, a consequence anticipated on intuitive grounds. (ii) If the value of the contact probability *P* and the threshold distance *r*~c~ are known precisely, then the distribution of the spatial distance can be readily computed by solving Eq. ([1](#Equ1){ref-type=""}) numerically. In Fig. [1b](#Fig1){ref-type="fig"}, we show the comparison between theory (Eq. ([1](#Equ1){ref-type=""})) and simulations (see the Methods section for details). The simulated curves are computed as follows: first collect (*P*~*mn*~, 〈*R*~*mn*~〉) for every pair labeled (*m*, *n*), where *P*~*mn*~ and 〈*R*~*mn*~〉 are computed using Eqs. ([17](#Equ17){ref-type=""}) and ([18](#Equ18){ref-type=""}) in the Methods section. The total number of pairs is *N*(*N*−1)/2. We then bin the points over the values of *P*~*mn*~. Finally, the mean value of 〈*R*~*mn*~〉 for each bin, 〈*R*〉 = *E*\[〈*R*~*mn*~〉\], is computed where E\[⋯\] is the binned average, which is computed using $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$(1/N_i)\mathop {\sum}\limits_{j = 1}^{N_i} {\left\langle {R_{mn}} \right\rangle ^j}$$\end{document}$, where *N*~*i*~ is the number of points in the *i*th bin. The bin size, Δ, is centered at *P*~*mn*~, spanning *P*~*mn*~−Δ/2 ≤ *P*~*mn*~ ≤ *P*~*mn*~ + Δ/2. Using this procedure, we find (Fig. [1](#Fig1){ref-type="fig"}) that the theory and simulations are in perfect agreement, which validates the theoretical result.Fig. 1Simulations demonstrate the power-law relation between contact probability and mean spatial distance and the effect of *r*~c~ on the inferred spatial distances. **a** A sketch of the Generalized Rouse Model for Chromosome (GRMC). Each bead represents a locus with a given resolution. Dashed lines represent harmonic bonds between loop anchors. **b** Mean spatial distance 〈*R*〉 as a function of the contact probability *P*. The solid lines are obtained using Eq. ([1](#Equ1){ref-type=""}) for different values of *r*~c~ (shown in the figure), the threshold distance for contact formation. The dots are simulation results. The agreement between simulations and theory is excellent. Asymptotically 〈*R*〉 approaches *r*~c~*P*^−1/3^ (dashed lines). The threshold for contact is expressed in terms of *a*, which is the equilibrium bond length in Eq. ([15](#Equ15){ref-type=""}). **c** Illustration of the sensitivity of *r*~c~ in determining the mean spatial distance 〈*R*〉. Blue and yellow curves are computed by solving 〈*R*〉 (Eq. ([1](#Equ1){ref-type=""})) for a given contact probability *P*~*mn*~ = 10^−3^, and *r*~c~. The calculated 〈*R*~*mn*~〉 is used in Eq. ([10](#Equ10){ref-type=""}) to obtain the distribution of the spatial distance *P*(*R*~*mn*~). Blue and yellow curves are for the same value of *P*, but different *r*~c~ values

Contact distance *r*~c~ affects the inferred value of the spatial distance {#Sec4}
--------------------------------------------------------------------------

However, in practice, the elements *P*~*mn*~ are measured with (unknown) statistical errors, and the value of the contact threshold *r*~c~ is only estimated. In the Hi-C experiments, contact probabilities and *r*~c~ by implication, are determined by a series of steps that start with cross-linking spatially adjacent loci using formaldehyde, chopping the chromatin into fragments using restriction enzymes, ligating the fragments with biotin, followed by sequence matching using deep-sequencing methods^[@CR5]^. Because of the inherent stochasticity associated with the overall Hi-C scheme, as well as the unavoidable heterogeneity (only a fraction of cells has a specific contact and the contact could be dynamic) in the cell population, the relationship *P*~*mn*~ and 〈*R*~*mn*~〉 is not straightforward.

To illustrate how the uncertainty in *r*~c~ affects the determination of the spatial distance in GRMC even when the population is homogeneous (all cells have a specific contact), we plot the distributions of distance for *r*~c~ = 0.02, 0.03 μm in Fig. [1c](#Fig1){ref-type="fig"}. A small change in *r*~c~ (from 0.02 to 0.03 μm) completely alters the distance distribution *P*(*R*), and hence the mean spatial distance (from ≈0.2 to ≈0.3 μm). For the exactly solvable GRMC, this can be explained by noting that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle R_{mn}\rangle \approx r_{\mathrm{c}}P_{mn}^{ - 1/3}$$\end{document}$ for small *P*~*mn*~. Because *P*~*mn*~ appears in the denominator, any uncertainty in *r*~c~ is amplified by *P*~*mn*~, especially when *P*~*mn*~ is small.

Heterogeneity causes paradox between FISH and Hi-C {#Sec5}
--------------------------------------------------

The expectation that the contact probability should decrease as the mean distance between the loci increases, which is the case in the exactly solvable ideal GRMC (*P*~*mn*~ ≈*r*~c~〈*R*~*mn*~〉^−3^), is sometimes violated when the experimental data^[@CR6]^ is analyzed^[@CR13],[@CR15]^. The paradox is a consequence of heterogeneity due to the existence of more than one population of cells, which implies that in some fraction of cells, contact between two loci exists while in others it is absent. Each distinct population has its own statistics. For instance, the probability distribution of the spatial distance between the *m*th and the *n*th loci, ~*Pi*,*mn*~(*r*), for one population of cells could be different from another population of cells *P*~*j*,*mn*~(*r*), where *i* and *j* are the indices for the two different populations (Fig. [2a](#Fig2){ref-type="fig"}). The Hi-C experiments yield only an average value of the contact probability. Let us illustrate the consequence of the inevitable heterogeneous mixture of cell populations by considering the simplest case in which only two distinct populations, one with probability *η* and the other 1 − *η*, are present (a generalization is presented below). For instance, in one population of cells, there is a CTCF loop between *m* and *n*, and it is absent in the other population. The probability distribution of spatial distance between the *m*th and the *n*th loci is a superposition of distributions for each population. Using Eq. ([10](#Equ10){ref-type=""}), the mixed distribution can be written as,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$P(R_{mn} = r) = \sqrt {\frac{2}{\pi }} \left( {\eta \frac{{r^2}}{{\sigma _{1,mn}^3}}e^{ - \frac{{r^2}}{{2\sigma _{1,mn}^2}}} + (1 - \eta )\frac{{r^2}}{{\sigma _{2,mn}^3}}e^{ - \frac{{r^2}}{{2\sigma _{2,mn}^2}}}} \right),$$\end{document}$$where *σ*~1,*mn*~ and *σ*~2,*mn*~ are the parameters with different values characterizing the two populations. In the GRMC, *σ*~1,*mn*~ and *σ*~2,*mn*~ are related to the mean spatial distances in the two populations by $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle R_{1,mn}\rangle = 2\sqrt {2/\pi } \sigma _{1,mn}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle R_{2,mn}\rangle = 2\sqrt {2/\pi } \sigma _{2,mn}$$\end{document}$. The mean spatial distance is, 〈*R*~*mn*~〉 = *η*〈*R*~1,*mn*~〉 + (1 − *η*)〈*R*~2,*mn*~〉, and the contact probability is *P*~*mn*~ = *ηP*~1,*mn*~ + (1 − *η*)*P*~2,*mn*~, where *P*~1,*mn*~ and *P*~2,*mn*~ are the contact probabilities for each population, given by Eq. ([1](#Equ1){ref-type=""}), which depends on the values of 〈*R*~1,*mn*~〉 and 〈*R*~2,*mn*~〉 as well as *r*~c~.Fig. 2Illustrating the FISH--Hi-C (\[*P*~*mn*~, 〈*R*~*mn*~〉\]) paradox. **a** Schematic illustration of the populations of two cells. There are two pairs of loci, pair 1 and pair 2. Cells 1 and 2 belong to two distinct populations such that pair 1 and pair 2 have different distributions of distances in the two cells. Pair 1 is always in proximity (contact is formed) in cell 1, whereas it is spatially separated (mean distance \>*r*~c~) in cell 2. Pair 2, on the other hand, has similar distributions of spatial distance in cells 1 and 2. Cell with two different populations gives rise to paradoxical behavior, which is illustrated by choosing *η*~1~ = 0.4 and *η*~2~ = 1 − *η*~1~ = 0.6. These are the probabilities for a cell belonging to population 1 and 2, respectively. The pair 1 has parameters *σ*~1~ = 0.3 μm and *σ*~2~ = 0.8 μm. The pair 2 has parameters *σ*~1~ = 0.4 μm and *σ*~2~ = 0.5 μm. See Eq. ([2](#Equ2){ref-type=""}) for the definition of *σ*~1~ and *σ*~2~. **b** The distribution of distance for pair 1 (thick blue) and pair 2 (thick orange), respectively. The distributions for the two different populations are shown separately for pair 1 (dashed lines) and pair 2 (dotted lines). **c ** Cumulative distribution of the spatial distance. The horizontal dashed line indicates the median distance. **d** Mean distances for pair 1 is larger than for pair 2. **e** Pair 1 has larger contact probability than 2, which is paradoxical since the distance between the loci in pair 1 is larger than in 2. The threshold for determining contact is *r*~c~ = 20 nm

If the values of 〈*R*~1,*mn*~〉 and 〈*R*~2,*mn*~〉 are unknown (as is the case in Hi-C experiments), and only the value of the contact probability between the two loci is provided, one can not uniquely determine the values of the mean spatial distances. This is the origin of the Hi-C and FISH data paradox. In Fig. [2b--e](#Fig2){ref-type="fig"}, we show an example of the paradox for a particular set of parameters (*η*, *σ*~1,*mn*~, *σ*~2,*mn*~). Pair \#1 has a larger contact probability than pair \#2, while also exhibiting a larger mean spatial distance. The GRMC explains in simple terms the origin of the paradox.

To systematically explore the parameter space, we display 〈*R*~*mn*~〉 and *P*~*mn*~ as heatmaps showing 〈*R*〉~1,*mn*~ versus 〈*R*〉~2,*mn*~ for different values of *η* (Fig. [3](#Fig3){ref-type="fig"}). When there is a single homogenous population (*η* = 0.0), the mean spatial distance 〈*R*~*mn*~〉 and contact probability *P*~*mn*~ depend only on the value of 〈*R*~2,*mn*~〉 (upper panel in Fig. [3](#Fig3){ref-type="fig"}). In this case, there is a precise one-to-one mapping between 〈*R*~*mn*~〉 and *P*~*mn*~. However, if *η* ≠ 0 (*η* = 0.3, lower panel in Fig. [3](#Fig3){ref-type="fig"}) then the relation between *P*~*mn*~ and 〈*R*~*mn*~〉 is complicated. The contour lines for *P*~*mn*~ cross the contour lines of 〈*R*~*mn*~〉, which implies that for a given value of *P*~*mn*~, one cannot infer the value of 〈*R*~*mn*~〉 without knowing the value of *η*, 〈*R*~1,*mn*~〉, and 〈*R*~2,*mn*~〉. For instance, the triangle and circle shown for *η* = 0.3 in Fig. [3](#Fig3){ref-type="fig"} demonstrate an example of the paradox, in which $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle R(\blacktriangledown )\rangle ( = 57a) > \langle R( \bullet )\rangle ( = 40a)$$\end{document}$ whereas $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$P(\blacktriangledown )( \approx 7.7 \times 10^{ - 4}) > P( \bullet )( \approx 3.9 \times 10^{ - 4})$$\end{document}$.Fig. 3Plots of mean distance 〈*R*~*mn*~〉 and the contact probability *P*~*mn*~ as heatmaps computed using *r*~c~ = 2*a*. The colorbars on the right show the values of 〈*R*~*mn*~〉 and *P*~*mn*~. The results for *η* = 0(≠0) is shown on the top panel (bottom panel). Two specific pairs are marked as triangle and circle in the lower-left panel. These loci pairs illustrate the \[*P*~*mn*~, 〈*R*~*mn*~〉\] paradox

Extracting cell subpopulation information from FISH data {#Sec6}
--------------------------------------------------------

Can we extract the information about subpopulations from experimental data so that the result from two vastly different techniques can be reconciled? To answer this question, we first generalize our theory for the GRMC to real chromatins. The generalization of Eq. ([2](#Equ2){ref-type=""}) is,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$P(R_{mn} = r) = \eta P(r|\langle R_{1,mn}\rangle ) + (1 - \eta )P(r|\langle R_{2,mn}\rangle ),$$\end{document}$$where *P*(*r*\|〈*R*~1,*mn*~〉) and *P*(*r*\|〈*R*~2,*mn*~〉) are the Redner-des Cloizeaux distribution of distances for polymers^[@CR27],[@CR28]^ (Supplementary Note [1](#MOESM1){ref-type="media"} and Supplementary Fig. [1](#MOESM1){ref-type="media"}). The distribution *P*(*r*\|〈*R*~*mn*~〉) is rigorously known for self-avoiding homopolymer in a good solvent, generalized Rouse model (Eq. ([10](#Equ10){ref-type=""}) in the Methods section), and a semi-flexible polymer^[@CR29],[@CR30]^. However, a simple analytic expression for chromosomes is not known. By assuming that the Redner-des Cloizeaux form for *P*(*r*\|〈*R*~*mn*~〉) also holds for chromosomes (see Supplementary Eq. ([1)](#MOESM1){ref-type="media"} for details), we find that *g* = 1 and *δ* = 5/4 in Supplementary Eq. ([1)](#MOESM1){ref-type="media"}. These parameters were previously extracted using the experimental data^[@CR8]^, and the Chromosome Copolymer Model (CCM) for chromosomes^[@CR24]^. The value of *g* is inferred from the scaling relationship between mean spatial distance 〈*R*〉 and contact probability *P*, *P*  \~〈*R*〉^3+*g*^. The value of *δ* is computed as *δ* = 1/(1 − *ν*). *ν* is inferred from scaling 〈*R*(*s*)〉 \~*s*^*ν*^, where *s* is the genomic distance.

The integral of Eq. ([3](#Equ3){ref-type=""}) up to *R*, which is the cumulative distribution function CDF(*R*), can be used to fit the FISH data. Thus, the probability of contact formation can be computed as, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\int_0^{r_{\mathrm{c}}}} P (r|\langle R\rangle ){\mathrm{d}}r$$\end{document}$, where *r*~c~ is the contact threshold. Using the data in ref. ^[@CR6]^, the CDF(*R*) for two pairs of loci are shown in Fig. [4a](#Fig4){ref-type="fig"}. By fitting the two experimentally measured curves to the theoretical prediction (see Supplementary Note [2](#MOESM1){ref-type="media"}), we obtain *η* ≈0.42 for peak4-loop and *η* ≈0.97 for peak3-control. The parameters obtained can then be used to compute the contact probability. Since the Hi-C experiments measure the number of contact events instead of contact probability and the value of *r*~c~ is unknown, we compare the relative contact frequency, which is computed as ~*Pi*~/〈*P*〉, where ~*Pi*~ is the contact probability computed using the model or the contact number measured in Hi-C for the *i*th pair and 〈*P*〉 is the mean value for all the pairs considered. First, we fit all the eight CDF(*R*) curves in ref. ^[@CR6]^. the excellent agreement between theory and experiments is vividly illustrated in Supplementary Fig. [2](#MOESM1){ref-type="media"} and also manifested by the Kolmogorov--Smirnov statistics (Supplementary Note [5](#MOESM1){ref-type="media"} and Supplementary Table [1](#MOESM1){ref-type="media"}). Second, we calculate their corresponding relative contact frequency (Fig. [4b](#Fig4){ref-type="fig"}). Comparison of the theoretical calculations with Hi-C measurements shows excellent agreement (Fig. [4b](#Fig4){ref-type="fig"}) with the Pearson correlation coefficient being 0.87. The contact probability is computed using *r*~c~ = 10 nm. Note that any value of *r*~c~ ≤ 10 nm gives similar results (Supplementary Fig. [3](#MOESM1){ref-type="media"}). The goodness of fits using different sets of *g* and *δ* is summarized in Supplementary Table [2](#MOESM1){ref-type="media"}. The set of *g* = 0 and *δ* = 2 gives equivalent good fits as the set of *g* = 1 and *δ* = 5/4. It is also important to note that fitting the FISH data with the assumption that cell population is homogeneous leads to unphysical values of *g* and *δ* and the Kolmogorov--Smirnov statistics are inferior (see Supplementary Note [4](#MOESM1){ref-type="media"}, Supplementary Fig. [5](#MOESM1){ref-type="media"} and Supplementary Table [3](#MOESM1){ref-type="media"}).Fig. 4Extracting statistics of subpopulations from FISH data. **a** Cumulative distribution function of the spatial distance, CDF(*R*) for two pairs of loci, labeled peak3-control and peak4-loop in ref. ^[@CR6]^. The excellent agreement between theory and experiments shows the usefulness of the relationship between *P*~*mn*~ and *R*~*mn*~ obtained using GRMC. The solid curves are the experiment data^[@CR6]^. The dashed lines are the fits to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\int_0^R} P (r){\mathrm{d}}r$$\end{document}$ (the needed expressions are in Eq. ([3](#Equ3){ref-type=""}) and Supplementary Eq. ([1)](#MOESM1){ref-type="media"}). The best fit parameters are *η*~peak3−control~ ≈ 0.97, 〈*R*~1,peak3−control~〉 ≈ 0.67 μm, 〈*R*~2,peak3−control~〉 ≈ 4.08 μm, *η*~peak4−loop~ ≈ 0.42, 〈*R*~1,peak4−loop~〉 ≈ 0.30 μm and 〈*R*~2,peak4−loop~〉 ≈ 1.21 μm. **b** Relative contact frequency computed from the fits of CDF(*R*) for eight pairs of loci investigate experimentally^[@CR6]^ (orange bars). For each pair of loci, the contact probability is calculated as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$P_{mn} = {\int}_0^{r_{\mathrm{c}}} P (r){\mathrm{d}}r$$\end{document}$ (Eq. ([3](#Equ3){ref-type=""})) using the parameters obtained by fitting CDF(*R*) with *r*~c~ = 20 nm. Comparison of the CDF(*R*)s between theory and experiments for the eight pairs of loci are displayed in Supplementary Fig. [2](#MOESM1){ref-type="media"}. Blue bars are computed using the contact number from Hi-C measurements in ref. ^[@CR6]^. The relative contact frequency is calculated as ~*Pi*~/〈*P*〉, where ~*Pi*~ is the contact probability computed using the model or the contact number measured in Hi-C for *i*th pair, and 〈*P*〉 is the mean value for all the pairs considered. p1-loop/p1-control/... are the ones referred to peak1-loop/peak1-control/... in ref. ^[@CR6]^

Interestingly, the values of 〈*R*~1~〉 obtained from fitting the four CTCF/cohesin-mediated loops (peak(1, 2, 3, 4)-loop) are all about 0.25−0.35 μm (*R*~1,peak1--loop~ ≈0.24 μm, *R*~1,peak2--loop~ ≈0.33  μm, *R*~1,peak3--loop~ ≈0.35 μm, *R*~1,peak4--loop~ ≈0.30 μm) regardless of their genomic separation (see Supplementary Table [1](#MOESM1){ref-type="media"}), suggesting that the mechanism of looping between CTCF motifs are similar with a mean spatial distance ≈0.3 μm. The physically reasonable value of 〈*R*~*mn*~〉 ≈0.3 νm for all peak--loop pairs shows that these CTCF-mediated contacts describe molecular interactions between loci that are separated by a few hundred kilo base pairs. It has been shown that these contacts, referred to as peaks^[@CR6]^ are significantly closer in space than others that are separated by similar genomic distance. The peak--loop contacts correspond to chromatin loops with the loci in the peaks being the anchor points between a specific loop. In sharp contrast, the distances between peak*i* and control (*i* goes from 1 to 4), which are greater than the distances between peak loci, vary ranging from ≈0.47 to ≈0.67 μm (see Supplementary Table [1](#MOESM1){ref-type="media"}). It is likely that these contacts are more dynamic because they are not be anchored by CTCF-binding proteins.

Massive heterogeneity in chromosome organization {#Sec7}
------------------------------------------------

In a recent study^[@CR19]^, which combined Hi-C and high-throughput optical imaging to map contacts within single chromosomes in human fibroblasts, revealed massive heterogeneity. Such extensive existence of a large number of conformations, leading to multiple or nearly continuous distribution of subpopulations, was much greater than previously anticipated. Although, the results in Fig. [4](#Fig4){ref-type="fig"} quantitatively reveal heterogeneity associated with CTCF loops by considering only two dominant subpopulations, the most recent experiment requires a generalization of the theory. In principle, our theory also applies to interactions of any nature, not only the CTCF loops. In doing so, it may be more reasonable to assume a continuous distribution of subpopulations, *P*(〈*R*〉), (see Supplementary Notes [6](#MOESM1){ref-type="media"} and [7](#MOESM1){ref-type="media"}, and Supplementary Fig. [6](#MOESM1){ref-type="media"} for generalization) instead of two discrete subpopulations, 〈*R*~1~〉 and 〈*R*~2~〉, which of course is much simpler and may suffice in many cases as the results in Fig. [4](#Fig4){ref-type="fig"} illustrate. As a proof of concept of our theory, we solve *P*(〈*R*〉) for the eight pairs of contacts analyzed in the previous section. The results are shown in Supplementary Fig. [7](#MOESM1){ref-type="media"}. In all cases, *P*(〈*R*〉)s are found to be multimodal. For peak1/2/3/4-control and peak3-loop, *P*(〈*R*〉) yield peaks located at positions very close to 〈*R*~1~〉 and 〈*R*~2~〉 shown in Supplementary Table [1](#MOESM1){ref-type="media"}, justifying the effectiveness of the theory. To show that our theory has a broader range of applicability, we use the FISH data from the recent study^[@CR19]^, which reports spatial distance measurements for 212 pairs of loci. *P*(〈*R*〉) is solved for each of a total of 212 pairs of loci. To illustrate our results, we compare in Fig. [5](#Fig5){ref-type="fig"} the predicted CDF(*r*) and the experimentally measured CDF(*r*), as well as the *P*(〈*R*〉) obtained by fitting for six pairs of loci as examples in Fig. [5](#Fig5){ref-type="fig"}. The results show substantial variations in 〈*R*〉, manifested by the multiple peaks and wide spread variations in *P*(〈*R*〉). Remarkably, the calculated CDF(*r*) (without any adjustable parameters) and the measured CDF(*r*) are in excellent agreement for the six loci pairs, which were arbitrarily chosen for illustration purposes. The residual errors between the two, shown as insets in Fig. [5](#Fig5){ref-type="fig"}, are extremely small.Fig. 5Exampled fits of CDF(*r*) using Supplementary Eq. ([9)](#MOESM1){ref-type="media"} to the experimental data^[@CR19]^. The six exampled pairs of loci are indicated above each subfigure. Orange lines, showing the fits using our theory, are indistinguishable from the experiment (the differences between fitted and experimental curve are shown in the insets). The distribution *P*(〈*R*〉) given in the integral equation (Supplementary Eq. ([9)](#MOESM1){ref-type="media"}) is solved using nonnegative Tikhonov Regularization (Supplementary Note [7](#MOESM1){ref-type="media"}). As shown here, *P*(〈*R*〉) have multi-peaks and are widespread, which is a manifestation of heterogeneity. We set *g* = 1 and *δ* = 5/4

In Fig. [6a,](#Fig6){ref-type="fig"} we show the normalized distributions *P*(〈*R*〉/*μ*(〈*R*〉)) for each of the 212 pairs of loci (see Supplementary Fig. [8](#MOESM1){ref-type="media"} for each pair as a separate figure). We expect that *P*(〈*R*〉/*μ*(〈*R*〉)) should be narrowly distributed around value 1 if there is only one population. However, many *P*(〈*R*〉/*μ*(〈*R*〉)) show multiple peaks with large variations. To further quantify the extent of heterogeneity, we calculate the coefficient of variation, CV = *σ*(〈*R*〉)/*μ*(〈*R*〉), where *σ*(〈*R*〉) and *μ*(〈*R*〉) are the standard deviation and the mean of 〈*R*〉, respectively. If there is only one population associated with 〈*R*〉, CV should have a value of around zero. Figure [6b](#Fig6){ref-type="fig"} shows the histogram of CV for all 212 pairs of loci. The CV values are widely distributed, suggesting that 3D structural heterogeneity is common and is associated with many pairs of loci rather than a few. Thus, the analyses of experimental data are not possible without taking heterogeneity into account. The theory presented here is sufficiently general and simple that it can be used to calculate the measurable quantities readily.Fig. 6Chromosome conformations are extensively heterogeneous. **a** Normalized distribution *P*(〈*R*〉/*μ*(〈*R*〉)) (*μ*(〈*R*〉) is the mean of 〈*R*〉) for all the 212 pairs of loci reported in ref. ^[@CR19]^. For almost every pair of loci, the associated *P*(〈*R*〉/*μ*(〈*R*〉)) has multiple peaks, and is widespread. **b** Histogram of the coefficient of variations CV for all 212 pairs of loci probed in ref. ^[@CR19]^. The CV values are calculated for each pair of loci, using CV = *σ*(〈*R*〉)/*μ*(〈*R*〉), where *σ*(〈*R*〉) is the standard deviation of 〈*R*〉. For a large number of loci pairs, CV exceeds 0.5, which is a quantitative measure of the extensive heterogeneity noted in the experiment^[@CR19]^

The role of loop extrusion in chromosome heterogeneity {#Sec8}
------------------------------------------------------

What is the origin of heterogeneity in the individual cell populations? There are two possibilities. The first one is static heterogeneity: each subpopulation explores a distinct region of the genomic folding landscape (GFL) (Fig. [7a](#Fig7){ref-type="fig"}). The second is the dynamic heterogeneity. Each cell explores a local minimum of the GFL before transiting to another local minimum (Fig. [7b](#Fig7){ref-type="fig"}). The only assumption in the application of our theory to genome organization is that there must be more than one population of cells, which does not violate the observation that the Hi-C experiment report only the average contact probability over millions of cells. Dynamic looping would be an example of the dynamic heterogeneity where the CTCF/cohein-mediated loops are formed and broken dynamically on a fast time scale compared with the lifetime of a cell. Such a picture is supported by recent single-cell molecule experiments^[@CR31],[@CR32]^. The average residence time of CTCF/cohesin complex is shown to be in the range of a few to tens of minutes, which is much smaller compared with the time scale of the cell cycle (15--30 h). Loop extrusion model^[@CR33]--[@CR35]^ is another possible origin of dynamic heterogeneity. In the loop extrusion model, it is thought that cohesins extrude loops along the chromosome fiber, which could detach stochastically. At any given time, there would be many subpopulations, each characterized by a distinct set of loops in the chromosome. Indeed, our analyses of the most recent high-throughput optical imaging data lend credence to the notion that multiple subpopulations in chromosomes arise because of massive dynamic heterogeneity. Our theory also gives an indirect theoretical justification for the work in ref. ^[@CR15]^, in which the authors found the loop extrusion model could lead to the \[*P*~*mn*~, 〈*R*~*mn*~〉\] paradox.Fig. 7Schematic of the Genomic Folding Landscape (GFL). **a** Static heterogeneity: cell subpopulation occupies distinct local minima in the GFL, with each minimum representing a stable organization. The energy barrier is too large for transitions between different local minima on a biological time scale (one cell cycle). **b** Dynamical heterogeneity: the energy barrier between local minima on the landscape is small enough which allows the dynamic transition between different subpopulations. **c** Combination of two different types of heterogeneity. In all three scenarios, the \[*P*~*mn*~, 〈*R*~*mn*~〉\] paradox arises. The loci contacts are in orange. The polymer conformation sketches are not shown in this scenario due to insufficient space

Single-cell temporal information is necessary to determine whether the loops are static or dynamic or a combination of the two (Fig. [7c](#Fig7){ref-type="fig"}). Hence, the combination of the dynamic FISH technique such as CRISPR--dCas9 FISH and single-cell Hi-C would be crucial for us to fully understand the organization of genomes. Our theory provides a theoretically rigorous method based on polymer physics to connect the results from measurements using the two vastly different techniques.

Discussion {#Sec9}
==========

From polymer physics for single chains it follows that in a homogeneous system, the contact probability and mean 3D distances are linked, resulting in a power-law relation connecting the two quantities that can be measured using Hi-C and FISH techniques. However, the one-to-one mapping does not hold in Hi-C experiments because of the presence of a mixture of distinct cell subpopulations each characterized by its own statistics leads to heterogeneity, which in turn gives rise to the \[*P*~*mn*~, 〈*R*~*mn*~〉\] paradox. We show that the theory based on precisely solvable GRMC could be used to solve the paradox in practice. The theory can be readily used to analyze data from experiments, provided the FISH and Hi-C experiments are done under similar conditions^[@CR6]^. The central result of the theory in Eq. ([3](#Equ3){ref-type=""}) can be used to analyze the available sparse FISH data. We show that the fraction of cell subpopulations (*η* in Eq. ([3](#Equ3){ref-type=""})) and the generalization derived in Supplementary Note [6](#MOESM1){ref-type="media"} can be extracted by fitting the FISH data using our theory. From Eq. ([3](#Equ3){ref-type=""}), we calculate the Hi-C contact probabilities, thus establishing that the theory resolves the \[*P*~*mn*~, 〈*R*~*mn*~〉\] paradox.

In this work, we confine ourselves to two-point interactions, which allows us to consider one pair of loci at a time. However, recent experiments probing multipoint interactions have suggested that formations of loops are likely to be cooperative^[@CR9],[@CR36]^, such that the formation of one loop could facilitate the formation of a nearby loop. Such cooperative loop formation was previously shown in an entirely different context involving the folding of proteins directed by disulfide bond formation^[@CR37]^. It can be shown within our framework that the formation of one loop can certainly increase the probability of formation of another loop. The theoretical basis for this finding is given in the Supplementary Note [8](#MOESM1){ref-type="media"}.

The reconciliation of the FISH and Hi-C data using polymer physics concepts is the first key step in integrating the data from these experimental techniques to construct the 3D structures of chromosomes. The work described here provides a theoretical basis for accomplishing this important task. Finally, our results suggest that heterogeneity in contact formation is an intrinsic property of genome organization, and hence the acquisition of single-cell experimental data is crucial to further our understanding of both the dynamics and the heterogeneous structural organization of chromosomes.

Methods {#Sec10}
=======

Generalized rouse model for chromosome {#Sec11}
--------------------------------------

In order to derive an approximate relationship connecting contact probabilities between loci and the three-dimensional distances, we use a variant of the random loop model^[@CR38],[@CR39]^. We first consider a minimal cross-linked phantom chain model, which incorporates the presence of CTCF/cohein-mediated loops^[@CR6]^. The model, originally introduced for describing physical gels^[@CR38]^, and more recently used for chromosome dynamics in a number of insightful studies^[@CR23],[@CR39]^, could be viewed as a Generalized Rouse Chromosome Model (GRMC)^[@CR40],[@CR41]^. The cross-links modeling the CTCF/cohein-mediated loops here are not random. Their locations are predetermined by the Hi-C data^[@CR6]^.

The equation of motion for the GRMC is^[@CR42]^$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\xi \frac{{{\mathrm{d}}{\mathbf{R}}}}{{{\mathrm{d}}t}} = {\mathbf{AR}} + {\mathbf{F}},$$\end{document}$$where *ξ* is the friction coefficient, **R** = \[**r**~1~, **r**~2~, ..., **r**~*N*~\]^T^ with **r**~*i*~ being the position of the *i*th locus. The vector **F** = \[**f**~1~, **f**~2~, ..., **f**~*N*~\]^T^ (T is the transpose), where **f**~*i*~ is the Gaussian random force acting on the *i*th locus, characterized by 〈*f*~*n*~(*t*)〉 = 0 and 〈*f*~*nα*~(*t*)*f*~*mβ*~(*t*′)〉 = 2*ξk*~B~*Tδ*~*nm*~*δ*~*αβ*~*δ*(*t* − *t*′); **A** is the *N* × *N* connectivity matrix, embedding the information of chain connectivity and the location of the loops connecting two loci (Fig. [1a](#Fig1){ref-type="fig"})$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$A_{mn} = \left\{ {\begin{array}{*{20}{l}} { - 2\kappa - |\Sigma _m|\omega ,} \hfill & {{\mathrm{if}}\,m = n \ne 1\;{\mathrm{or}}\;N} \hfill \\ { - \kappa - |\Sigma _m|\omega ,} \hfill & {{\mathrm{if}}\,m = n = 1{\mathrm{ or }}N} \hfill \\ {\kappa ,} \hfill & {{\mathrm{if}}\,|m - n| = 1} \hfill \\ {\omega ,} \hfill & {{\mathrm{if}}\,|m - n| > 1,{\mathrm{and}}\,{\mathrm{connected}}\,{\mathrm{in}}\,\Sigma } \hfill \\ {0,} \hfill & {{\mathrm{if}}\,{\mathrm{otherwise}}} \hfill \end{array}} \right.$$\end{document}$$where Σ is the set of indices representing the loci pairs specifying the CTCF facilitated loop anchors, and \|Σ~*m*~\| is the number of loops connected to the *m*th locus. The spring constant *κ* enforces chain connectivity, and *ω* is the associated spring constant for a CTCF pair. Note that the GRMC model does not account for excluded volume interactions, which in the modeling of chromatin is often justified by noting that topoisomerases enable chain crossing. Our purpose is to use GRMC to first illustrate concretely the challenges in going from the measured average contact map to spatial organization, precisely. More importantly, using the insights from the study of the GRMC, we provide a solution to the FISH--Hi-C paradox.

Since **A** in Eq. ([5](#Equ5){ref-type=""}) is a real symmetric matrix, it can be diagonalized using the orthonormal matrix **V**$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathbf{VAV}}^T = {\mathbf{\Lambda}} = {\mathrm{diag}}(\lambda _0,\lambda _1,...,\lambda _{N - 1}),$$\end{document}$$where *λ*~0~, *λ*~1~, ..., *λ*~*N*−1~ are the eigenvalues of **A**. By defining **X** = **VR** and using **R** = **V**^*T*^**X** and **VV**^*T*^ = **I**, we obtain the equations of motion of the normal coordinates **X**$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\xi \frac{{\mathrm{d}}{\mathbf{X}}}{\mathrm{d}}t = {\mathbf{\Lambda}}{\mathbf{X}} + {\mathbf{f}}.$$\end{document}$$Because **Λ** is a diagonal matrix, the normal coordinates of the GRMC **X**~*p*~ are decoupled. Using the normal modes, **X**, the physical quantities associated with the polymer can be readily calculated. Therefore, for GRMC with a predetermined set of CTCF/cohein-mediated loops, we can solve for the eigenvalues of the connectivity matrix **A**, and the orthonormal matrix **V** numerically, and thus calculate the contact probability and spatial distance precisely.

Relation between contact probability and mean spatial distance {#Sec12}
--------------------------------------------------------------

The vector between the positions of the *m*th and the *n*th loci may be written as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathbf{R}}_m - {\mathbf{R}}_n = \mathop {\sum}\limits_{p = 0}^{N - 1} {(V_{pm} - V_{pn})} {\mathbf{X}}_p,$$\end{document}$$where *V*~*pm*~ and *V*~*pn*~ are the elements of orthonormal matrix **V**. The equilibrium solution of Eq. ([7](#Equ7){ref-type=""}) yields, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\mathop {{\lim }}\nolimits_{t \to \infty } X_{p,\alpha }(t)\sim {\cal{N}}(0, - \frac{{k_{\mathrm{B}}T}}{{\lambda _p}})$$\end{document}$, where *α* = *x*, *y*, *z*, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\cal{N}}$$\end{document}$ is Gaussian distribution. Therefore$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\mathop {{\lim }}\limits_{t \to \infty } R_{mn,\alpha }(t)\sim {\cal{N}}(0, - \mathop {\sum}\limits_{p = 0}^{N - 1} {(V_{pm} - V_{pn})^2} \frac{{k_{\mathrm{B}}T}}{{\lambda _p}}) \equiv {\cal{N}}(0,\sigma _{mn,\alpha }^2).$$\end{document}$$where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\sigma _{mn,\alpha } = - \mathop {\sum}\limits_{p = 0}^{N - 1} {(V_{pm} - V_{pn})^2} (k_{\mathrm{B}}T/\lambda _p)$$\end{document}$. Since the model is isotropic, it follows that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\sigma _{mn,x}^2 = \sigma _{mn,y}^2 = \sigma _{mn,z}^2 \equiv \sigma _{mn}^2$$\end{document}$. The mean distance 〈*R*~*mn*~〉 is related to *σ*~*mn*~ through $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle R_{mn}\rangle = 2\sqrt {2/\pi } \sigma _{mn}$$\end{document}$. The distribution of the distance between the *m*th and the *n*th loci, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\mathop {{\lim }}\nolimits_{t \to \infty } |{\mathbf{R}}_{mn}(t)| = \mathop {{\lim }}\nolimits_{t \to \infty } \sqrt {\mathop {\sum}\limits_\alpha {R_{mn,\alpha }^2} (t)}$$\end{document}$ is a non-central chi distribution (we will neglect the notation $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\mathop {{\lim }}\nolimits_{t \to \infty }$$\end{document}$ from now on)$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$P(R_{mn} = r) = \sqrt {\frac{2}{\pi }} \frac{1}{{\sigma _{mn}}}e^{ - r^2/(2\sigma _{mn}^2)}\frac{{r^2}}{{\sigma _{mn}^2}}.$$\end{document}$$The contact probability *P*~*mn*~, for a given threshold *r*~c~ (contact exists if *r* ≤ *r*~c~), computed using Eq. ([10](#Equ10){ref-type=""}) yields$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\begin{array}{c}P_{mn} = {\int_0^{r_{\mathrm{c}}}} {\mathrm{d}} r\sqrt {\frac{2}{\pi }} \frac{1}{{\sigma _{mn}}}e^{ - r^2/(2\sigma _{mn}^2)}\frac{{r^2}}{{\sigma _{mn}^2}}\\ = {\mathrm{Erf}}\left( {\frac{{r_{\mathrm{c}}}}{{\sqrt 2 \sigma _{mn}}}} \right) - \sqrt {\frac{2}{\pi }} e^{ - \frac{{r_{\mathrm{c}}^2}}{{2\sigma _{mn}^2}}}\frac{{r_{\mathrm{c}}}}{{\sigma _{mn}}}.\end{array}$$\end{document}$$The mean spatial distance 〈*R*~*mn*~〉 is given by$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle R_{mn}\rangle = {\int_0^\infty} {\mathrm{d}} rr\sqrt {\frac{2}{\pi }} \frac{1}{{\sigma _{mn}}}e^{ - r^2/(2\sigma _{mn}^2)}\frac{{r^2}}{{\sigma _{mn}^2}} = 2\sqrt {\frac{2}{\pi }} \sigma _{mn}.$$\end{document}$$Using Eqs. ([11](#Equ11){ref-type=""}) and ([12](#Equ12){ref-type=""}), the desired relation between *P*~*mn*~ and 〈*R*~*mn*~〉 becomes$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$P_{mn} = {\mathrm{erf}}\left( {\frac{{2r_{\mathrm{c}}}}{{\sqrt \pi \langle R_{mn}\rangle }}} \right) - \frac{4}{\pi }\frac{{r_{\mathrm{c}}}}{{\langle R_{mn}\rangle }}e^{ - \frac{{4r_{\mathrm{c}}^2}}{{\pi \left\langle {R_{mn}} \right\rangle ^2}}} \equiv R_0(\langle R_{mn}\rangle ).$$\end{document}$$Eq. ([13](#Equ13){ref-type=""}) is identical with Eq. ([1](#Equ1){ref-type=""}) in the main text.

Simulations {#Sec13}
-----------

The energy function for the GRMC is$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$U({\mathbf{r}}_1,...,{\mathbf{r}}_N) = \mathop {\sum}\limits_{i = 1}^{N - 1} {U_i^{\mathrm{S}}} + \mathop {\sum}\limits_{\{ p,q\} } {U_{\{ p,q\} }^{\mathrm{L}}} .$$\end{document}$$For the bonded stretch potential, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$U_i^{\mathrm{S}}$$\end{document}$, we use$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$U_i^{\mathrm{S}} = \frac{\kappa }{2}(|{\mathbf{r}}_{i + 1} - {\mathbf{r}}_i| - a)^2,$$\end{document}$$where *a* is the equilibrium bond length. The interaction between the loop anchors is also modeled using a harmonic potential$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$U_{\{ p,q\} }^{\mathrm{L}} = \frac{\omega }{2}(|{\mathbf{r}}_p - {\mathbf{r}}_q| - a)^2,$$\end{document}$$where the spring constant is associated with the CTCF facilitated loops, and {*p*, *q*} represents the indices of the loop anchors, which are taken from the Hi-C data^[@CR6]^ (Supplementary Note [3](#MOESM1){ref-type="media"}). We simulate the chromosome segment from 146 to 158 Mbps of Chromosome 5. Each monomer represents 1200 bps, resulting the total number of coarse-grained loci *N* = 10,000.

In order to accelerate conformational sampling, we perform Langevin Dynamics simulations at low friction^[@CR43]^. We simulate each trajectory for 10^8^ time steps, and save the snapshots every 10,000 time steps. We generate ten independent trajectories, which are sufficient to obtain reliable statistics (Supplementary Fig. [4](#MOESM1){ref-type="media"}).

Data analyses {#Sec14}
-------------

The contact probability between the *m*th and *n*th loci in the simulation is calculated using$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$P_{mn} = \frac{1}{{TM}}\mathop {\sum}\limits_{a = 1}^M {\mathop {\sum}\limits_{t = 1}^T \Theta } (r_{\mathrm{c}} - |{\mathbf{r}}_m^{(a)}(t) - {\mathbf{r}}_n^{(a)}(t)|),$$\end{document}$$where Θ(⋅) is the Heaviside step function, *r*~c~ is the threshold distance for determining the contacts, the summation is over the snapshots along the trajectory, and the total *M* number of independent trajectories, and *T* is the number of snapshots for a single trajectory. The mean spatial distance between the *i*th and the *j*th loci in the simulation is calculated using$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle R_{mn}\rangle = \frac{1}{{TM}}\mathop {\sum}\limits_{a = 1}^M {\mathop {\sum}\limits_{t = 1}^T {|{\mathbf{r}}_m^{(a)}(t) - {\mathbf{r}}_n^{(a)}(t)|} } .$$\end{document}$$The objective is to go from *P*~*mn*~ to 〈*R*~*mn*~〉, and to determine, if in doing so, we get reasonable results. Because these quantities can be computed precisely in the GRMC, the \[*P*~*mn*~, 〈*R*~*mn*~〉\] relationship can be tested, which allows us to obtain the needed cues to solve the FISH--Hi-C paradox.

Reporting summary {#Sec15}
-----------------

Further information on research design is available in the [Nature Research Reporting Summary](#MOESM3){ref-type="media"} linked to this article.

Supplementary information
=========================

 {#Sec16}

Supplementary Information Peer Review File Reporting Summary

**Peer review information**: *Nature Communications* thanks Bogdan Bintu and other anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information
=========================

**Supplementary Information** accompanies this paper at 10.1038/s41467-019-11897-0.

We are grateful to the National Science Foundation (CHE 19-00093) and the Collie-Welch Regents Chair (F-0019) for supporting this work.

G.S. and D.T. designed and performed the research, G.S. and D.T. analyzed the data, G.S. and D.T. wrote the paper.

All relevant data supporting the findings of this study are available within the article and its [Supplementary Information](#MOESM1){ref-type="media"} files or upon requests from the corresponding author. The Hi-C and FISH experimental data used in this study are publicly available from GEO database under accession number [GSE63525](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63525) and from 4DN portal at <https://data.4dnucleome.org/publications/80007b23-7748-4492-9e49-c38400acbe60/>. The processed data are available upon request from the authors.

The polymer simulations are performed using LAMMPS Molecular Dynamics Simulation software^[@CR44]^, which is an open-source code available at <http://lammps.sandia.gov>. The codes used to analyze data in the present study are deposited to Github repository <https://github.com/anyuzx/chromosome-heterogeneity-analysis>.

The authors declare no competing interests.
