Conformational Heterogeneity in Human Interphase Chromosome Organization Reconciles the FISH and Hi-C Paradox

TL;DR

This paper resolves the apparent contradiction between FISH and Hi-C data on chromosome organization by introducing a heterogeneity-based theory that explains contact probability and spatial distance variations across cell populations.

Contribution

It presents a novel, parameter-free theoretical framework based on the Generalized Rouse Model to interpret FISH and Hi-C data considering cell heterogeneity.

Findings

Heterogeneity explains FISH-Hi-C paradox.

The model accurately reproduces Hi-C data from FISH distributions.

Genome organization varies significantly across cells.

Abstract

Hi-C experiments are used to infer the contact probabilities between loci separated by varying genome lengths. Contact probability should decrease as the spatial distance between two loci increases. However, studies comparing Hi-C and FISH data show that in some cases the distance between one pair of loci, with larger Hi-C readout, is paradoxically larger compared to another pair with a smaller value of the contact probability. Here, we show that the FISH-Hi-C paradox can be resolved using a theory based on a Generalized Rouse Model for Chromosomes (GRMC). The FISH-Hi-C paradox arises because the cell population is highly heterogeneous, which means that a given contact is present in only a fraction of cells. Insights from the GRMC is used to construct a theory, without any adjustable parameters, to extract the distribution of subpopulations from the FISH data, which quantitatively reproduces the Hi-C data. Our results show that heterogeneity is pervasive in genome organization at all length scales, reflecting large cell-to-cell variations.

Equations

$$P_{mn}=\mathrm{erf}\bigg{(}\frac{2r_{\mathrm{c}}}{\sqrt{\pi}\langle R_{mn}\rangle}\bigg{)}-\frac{4}{\pi}\frac{r_{\mathrm{c}}}{\langle R_{mn}\rangle}e^{-\frac{4r_{\mathrm{c}}^{2}}{\pi\langle R_{mn}\rangle^{2}}}\equiv R_{0}(\langle R_{mn}\rangle).$$

$$P(R_{mn}=r)=\sqrt{\frac{2}{\pi}}\bigg{(}\eta\frac{r^{2}}{\sigma^{3}_{1,mn}}e^{-\frac{r^{2}}{2\sigma^{2}_{1,mn}}}+(1-\eta)\frac{r^{2}}{\sigma^{3}_{2,mn}}e^{-\frac{r^{2}}{2\sigma^{2}_{2,mn}}}\bigg{)}$$

$$P(R_{mn}=r)=\eta P(r|⟨R_{1,mn}⟩)+(1-\eta )P(r|⟨R_{2,mn}⟩)$$

$$\xi \frac{\mathrm{d}\mathrm{R}}{\mathrm{d}t}=\mathrm{A}\mathrm{R}+\mathrm{F}$$

$$A_{mn}=\{\begin{array}{ll}-2\kappa -|{\Sigma }_m|\omega ,& \text{if }m=n\ne 1\text{ or }N\\ -\kappa -|{\Sigma }_m|\omega ,& \text{if }m=n=1\text{ or }N\\ \kappa ,& \text{if }|m-n|=1\\ \omega ,& \text{if }|m-n|>1,\text{and connected in }\Sigma \\ 0,& \text{if otherwise}\end{array}$$

$$\mathrm{V}\mathrm{A}{\mathrm{V}}^T=\mathrm{\Lambda }=\mathrm{diag}({\lambda }_0,{\lambda }_1,...,{\lambda }_{N-1})$$

$$\xi \frac{\mathrm{d}\mathrm{X}}{\mathrm{d}t}=\mathrm{\Lambda }\mathrm{X}+\mathrm{f}.$$

$${\mathrm{R}}_m-{\mathrm{R}}_n=\sum _{p=0}^{N-1}(V_{pm}-V_{pn}){\mathrm{X}}_p$$

$$\underset{t\to \infty }{\lim }{R}_{mn,\alpha }(t)\sim \mathcal{N}(0,-\sum _{p=0}^{N-1}(V_{pm}-V_{pn}{)}^2\frac{k_{\mathrm{B}}T}{{\lambda }_p})\equiv \mathcal{N}(0,{\sigma }_{mn,\alpha }^2).$$

$$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}.$$

$$\begin{split}P_{mn}&=\int_{0}^{r_{\mathrm{c}}}\mathrm{d}r\sqrt{\frac{2}{\pi}}\frac{1}{\sigma_{mn}}e^{-r^{2}/(2\sigma^{2}_{mn})}\frac{r^{2}}{\sigma^{2}_{mn}}\\ &=\mathrm{Erf}\bigg{(}\frac{r_{\mathrm{c}}}{\sqrt{2}\sigma_{mn}}\bigg{)}-\sqrt{\frac{2}{\pi}}e^{-\frac{r_{\mathrm{c}}^{2}}{2\sigma_{mn}^{2}}}\frac{r_{\mathrm{c}}}{\sigma_{mn}}.\end{split}$$

$$⟨R_{mn}⟩={\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}.$$

$$P_{mn}=\mathrm{erf}\bigg{(}\frac{2r_{\mathrm{c}}}{\sqrt{\pi}\langle R_{mn}\rangle}\bigg{)}-\frac{4}{\pi}\frac{r_{\mathrm{c}}}{\langle R_{mn}\rangle}e^{-\frac{4r_{\mathrm{c}}^{2}}{\pi\langle R_{mn}\rangle^{2}}}\equiv R_{0}(\langle R_{mn}\rangle).$$

$$U({\mathrm{r}}_1,...,{\mathrm{r}}_N)=\sum _{i=1}^{N-1}{U}_i^{\mathrm{S}}+\sum _{\{p,q\}}{U}_{\{p,q\}}^{\mathrm{L}}.$$

$$U_i^{\mathrm{S}}=\frac{\kappa }{2}(|{\mathrm{r}}_{i+1}-{\mathrm{r}}_i|-a{)}^2,$$

$$U_{\{p,q\}}^{\mathrm{L}}=\frac{\omega }{2}(|{\mathrm{r}}_p-{\mathrm{r}}_q|-a{)}^2$$

$$P_{mn}=\frac{1}{TM}\sum _{a=1}^M\sum _{t=1}^T\Theta (r_{\mathrm{c}}-|{\mathrm{r}}_m^{(a)}(t)-{\mathrm{r}}_n^{(a)}(t)|),$$

$$⟨R_{mn}⟩=\frac{1}{TM}\sum _{a=1}^M\sum _{t=1}^T|{\mathrm{r}}_m^{(a)}(t)-{\mathrm{r}}_n^{(a)}(t)|.$$

Full text

Introduction

Through remarkable Hi-C experiments LiebermanAiden2009 ; Dixon2012 ; Sexton2012 ; Jin2013 ; Dekker2013 ; rao20143d , based on the Chromosome Conformation Capture (3C) technique dekker2002capturing , 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 to 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 wang2016spatial ; Bintu2018 . 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 Chen2013 ; Ma2015 ; Ma2016 . 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 giorgetti2016closing ; Fraser15MicrbiolMolbiolReview , 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 giorgetti2016closing ; fudenberg2017fish ; Bickmore2013 ; williamson2014spatial . It should be noted that in most cases, the Hi-C and FISH measurements agree very well wang2016spatial ; Bintu2018 ; szabo2018tads ; finn2019extensive . 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 giorgetti2016closing ; Fraser15MicrbiolMolbiolReview . Contact between two loci could be a rare event, not present in all cells, which is captured in 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 FISH technique. Thus, in order to combine the data from the two powerful techniques, it is crucial to establish a theoretical basis with potential 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 giorgetti2016closing ) insights into the discordance between the two methods, when they occur, can be obtained using polymer physics concepts. Recently, Fudenberg and Imakaev fudenberg2017fish performed polymer simulations using a strong attractive energy between two labelled loci and a ten fold 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 principle goal of this work.

In addition, recent single-cell Hi-C Stevens2017 ; Flyamer2017 ; Tan2018 and FISH experiments wang2016spatial ; szabo2018tads ; Bintu2018 ; finn2019extensive 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 shukron2017transient ; shi2018interphase ; liu2018chain ; Buckle2018 . Indeed, it has been shown, using Hi-C and FISH data as well simulations Buckle2018 , 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 Buckle2018 .

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 $\langle R\rangle$. 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 giorgetti2016closing ; fudenberg2017fish . In other words, the average spatial distance between $m$ and $n$ ($\langle R_{mn}\rangle$) could be larger than $\langle R_{kl}\rangle$, the distance between loci $k$ and $l$, even if $P_{mn}>P_{kl}$. These results provide a basis for understanding the origin of 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

Relating contact probability to mean spatial distance. The exact relationship between $P_{mn}$ (contact probability between $m^{th}$ and $n^{th}$ locus) and the corresponding mean spatial distance, $\langle R_{mn}\rangle$ for GRMC (see Methods section for details of the derivation) is,

[TABLE]

The inverse of $R_{0}(\langle R_{mn}\rangle)$, the solution to Equation (1), gives the mean spatial distance $\langle R_{mn}\rangle$ as a function of the contact probability $P_{mn}$. Note that $m$ and $n$ are arbitrary locations of any two loci, and thus Equation (1) is general for any pair of loci.

A couple of conclusions, relevant to the application to the chromosomes, follow from Equation (1). (i) Note that Equation (1) is an exact one-to-one relation between the mean distance $\langle R_{mn}\rangle$ and the contact probability $P_{mn}$ provided $r_{\mathrm{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 Equation (1) that $\langle R_{mn}\rangle\approx r_{\mathrm{c}}P_{mn}^{-1/3}$. For the ideal GRMC, this implies that for any $m,n,k,l$, if $P_{mn}<P_{kl}$ then $\langle R_{mn}\rangle>\langle R_{kl}\rangle$, a consequence anticipated on intuitive grounds. (ii) If the value of the contact probability $P$ and the threshold distance $r_{\mathrm{c}}$ are known precisely, then the distribution of the spatial distance can be readily computed by solving Equation (1) numerically. In Fig.1(b), we show the comparison between theory (Equation (1)) and simulations (see Methods section for details). The simulated curves are computed as follows: first collect $(P_{mn},\langle R_{mn}\rangle)$ for every pair labeled $(m,n)$ where $P_{mn}$ and $\langle R_{mn}\rangle$ are computed using Equations (17-18) 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 $\langle R_{mn}\rangle$ for each bin, $\langle R\rangle=E[\langle R_{mn}\rangle]$, is computed where $\mathrm{E}[\cdots]$ is the binned average, which is computed using $(1/N_{i})\sum_{j=1}^{N_{i}}\langle R_{mn}\rangle^{j}$ where $N_{i}$ is the number of points in the $i^{th}$ bin. The bin size, $\Delta$, is centered at $P_{mn}$, spanning $P_{mn}-\Delta/2\leq P_{mn}\leq P_{mn}+\Delta/2$. Using this procedure, we find (Fig.1) that the theory and simulations are in perfect agreement, which validates the theoretical result.

Contact distance $r_{\mathrm{c}}$ affects the inferred value of the spatial distance. However, in practice, the elements $P_{mn}$ are measured with (unknown) statistical errors, and the value of the contact threshold $r_{\mathrm{c}}$ is only estimated. In the Hi-C experiments, contact probabilities and $r_{\mathrm{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 Dekker2013 . 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 $\langle R_{mn}\rangle$ is not straightforward.

To illustrate how the uncertainty in $r_{\mathrm{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_{\mathrm{c}}=0.02,0.03\mathrm{\ }\mu\mathrm{m}$ in Fig.1(c). A small change in $r_{\mathrm{c}}$ (from 0.02 $\mu\mathrm{m}$ to 0.03 $\mu\mathrm{m}$) completely alters the distance distribution $P(R)$, and hence the mean spatial distance (from $\approx 0.2\mu m$ to $\approx 0.3\mu m$). For the exactly solvable GRMC, this can be explained by noting that $\langle R_{mn}\rangle\approx r_{\mathrm{c}}P_{mn}^{-1/3}$ for small $P_{mn}$. Because $P_{mn}$ appears in the denominator, any uncertainty in $r_{\mathrm{c}}$ is amplified by $P_{mn}$, especially when $P_{mn}$ is small.

Heterogeneity causes paradox between FISH and Hi-C. 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}\approx r_{\mathrm{c}}\langle R_{mn}\rangle^{-3}$), is sometimes violated when the experimental data rao20143d is analyzed giorgetti2016closing ; fudenberg2017fish . 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, $P_{i,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.2(a)). 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 $\eta$ and the other $1-\eta$, 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 Equation 10, the mixed distribution can be written as,

[TABLE]

where $\sigma_{1,mn}$ and $\sigma_{2,mn}$ are the parameters with different values characterizing the two populations. In the GRMC, $\sigma_{1,mn}$ and $\sigma_{2,mn}$ are related to the mean spatial distances in the two populations by $\langle R_{1,mn}\rangle=2\sqrt{2/\pi}\sigma_{1,mn}$ and $\langle R_{2,mn}\rangle=2\sqrt{2/\pi}\sigma_{2,mn}$. The mean spatial distance is, $\langle R_{mn}\rangle=\eta\langle R_{1,mn}\rangle+(1-\eta)\langle R_{2,mn}\rangle$, and the contact probability is $P_{mn}=\eta P_{1,mn}+(1-\eta)P_{2,mn}$ where $P_{1,mn}$ and $P_{2,mn}$ are the contact probabilities for each population, given by Equation (1), which depends on the values of $\langle R_{1,mn}\rangle$ and $\langle R_{2,mn}\rangle$ as well as $r_{\mathrm{c}}$.

If the values of $\langle R_{1,mn}\rangle$ and $\langle R_{2,mn}\rangle$ 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 Figs.2(b)-2(e) we show an example of the paradox for a particular set of parameters ($\eta,\sigma_{1,mn},\sigma_{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 $\langle R_{mn}\rangle$ and $P_{mn}$ as heat maps showing $\langle R\rangle_{1,mn}$ versus $\langle R\rangle_{2,mn}$ for different values of $\eta$ (Fig.3). When there is a single homogenous population ($\eta=0.0$), the mean spatial distance $\langle R_{mn}\rangle$ and contact probability $P_{mn}$ depend only on the value of $\langle R_{2,mn}\rangle$ (upper panel in Fig.3). In this case, there is a precise one-to-one mapping between $\langle R_{mn}\rangle$ and $P_{mn}$. However, if $\eta\neq 0$ ($\eta=0.3$, lower panel in Fig.3) then the relation between $P_{mn}$ and $\langle R_{mn}\rangle$ is complicated. The contour lines for $P_{mn}$ cross the contour lines of $\langle R_{mn}\rangle$, which implies that for a given value of $P_{mn}$, one cannot infer the value of $\langle R_{mn}\rangle$ without knowing the value of $\eta$, $\langle R_{1,mn}\rangle$ and $\langle R_{2,mn}\rangle$. For instance, the triangle and circle shown for $\eta=0.3$ in Fig.3 demonstrate an example of the paradox in which $\langle R(\blacktriangledown)\rangle(=57a)>\langle R(\bullet)\rangle(=40a)$ whereas $P(\blacktriangledown)(\approx 7.7\times 10^{-4})>P(\bullet)(\approx 3.9\times 10^{-4})$.

Extracting cell subpopulation information from FISH data. 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 derived from GRMC to real chromatins. The generalization of Equation (2) is,

[TABLE]

where $P(r|\langle R_{1,mn}\rangle)$ and $P(r|\langle R_{2,mn}\rangle)$ are the Redner-des Cloizeaux distribution of distances for polymers jannink1990polymers ; redner1980distribution (Supplementary Note 1 and Supplementary Figure 1). The distribution $P(r|\langle R_{mn}\rangle)$ is rigorously known for self-avoiding homopolymer in a good solvent, generalized Rouse model (Equation 10 in the Methods section), and a semi-flexible polymer Hyeon2006 ; Wilhelm1996 . However, a simple analytic expression for chromosomes is not known. By assuming that the Redner-des Cloizeaux form for $P(r|\langle R_{mn}\rangle)$ also holds for chromosomes (see Supplementary Equation 1 for details), we find that $g=1$ and $\delta=5/4$ in Supplementary Equation 1. These parameters were previously extracted using experimental data wang2016spatial , and the Chromosome Copolymer Model (CCM) for chromosomes shi2018interphase . The value of $g$ is inferred from the scaling relationship between mean spatial distance $\langle R\rangle$ and contact probability $P$, $P\sim\langle R\rangle^{3+g}$. The value of $\delta$ is computed as $\delta=1/(1-\nu)$. $\nu$ is inferred from scaling $\langle R(s)\rangle\sim s^{\nu}$ where $s$ is the genomic distance.

The integral of Equation (3) up to $R$, which is the cumulative distribution function $\mathrm{CDF}(R)$, can be used to fit the FISH data. Thus, the probability of contact formation can be computed as, $\int_{0}^{r_{\mathrm{c}}}P(r|\langle R\rangle)\mathrm{d}r$ where $r_{\mathrm{c}}$ is the contact threshold. Using the data in rao20143d , the $\mathrm{CDF}(R)$ for two pairs of loci are shown in Fig.4(a). By fitting the two experimentally measured curves to the theoretical prediction (see Supplementary Note 2), we obtain $\eta\approx 0.42$ for peak4-loop and $\eta\approx 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_{\mathrm{c}}$ is unknown, we compare the relative contact frequency, which is computed as $P_{i}/\langle P\rangle$ where $P_{i}$ is the contact probability computed using the model or the contact number measured in Hi-C for the $i^{th}$ pair and $\langle P\rangle$ is the mean value for all the pairs considered. First, we fit all the eight $\mathrm{CDF}(R)$ curves in rao20143d . the excellent agreement between theory and experiments is vividly illustrated in Supplementary Figure 2 and also manifested by the Kolmogorov-Smirnov statistics (Supplementary Note 5 and Supplementary Table 1). Second, we calculate their corresponding relative contact frequency (Fig.4(b)). Comparison of the theoretical calculations with Hi-C measurements shows excellent agreement (Fig.4(b)) with the Pearson correlation coefficient being 0.87. The contact probability is computed using $r_{\mathrm{c}}=10\mathrm{\ nm}$. Note that any value of $r_{\mathrm{c}}\leq 10\mathrm{\ nm}$ gives similar results (Supplementary Figure 3). The goodness of fits using different sets of $g$ and $\delta$ is summarized in Supplementary Table 2. The set of $g=0$ and $\delta=2$ gives equivalent good fits as the set of $g=1$ and $\delta=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 $\delta$ and the Kolmogorov-Smirnov statistics are inferior (see Supplementary Note 4, Supplementary Figure 5 and Supplementary Table 3).

Interestingly, the values of $\langle R_{1}\rangle$ obtained from fitting the four CTCF/cohesin mediated loops (peak(1,2,3,4)-loop) are all about $0.25-0.35\mu\mathrm{m}$ ($R_{1,\mathrm{peak1-loop}}\approx 0.24\mu m$, $R_{1,\mathrm{peak2-loop}}\approx 0.33\mu m$, $R_{1,\mathrm{peak3-loop}}\approx 0.35\mu m$,$R_{1,\mathrm{peak4-loop}}\approx 0.30\mu m$) regardless of their genomic separation (see Supplementary Table 1), suggesting that the mechanism of looping between CTCF motifs are similar with a mean spatial distance $\approx 0.3\mathrm{\ \mu m}$. The physically reasonable value of $\langle R_{mn}\rangle\approx 0.3\mathrm{\nu 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 rao20143d 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$-control ($i$ goes from 1 to 4), which are greater than the distances between peak loci, vary ranging from $\approx 0.47\mathrm{\ \mu m}$ to $\approx 0.67\mathrm{\ \mu m}$ (see Supplementary Table 1). It is likely that these contacts are more dynamic because they are not be anchored by CTCF binding proteins.

Massive heterogeneity in chromosome organization. In a recent study finn2019extensive , 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 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(\langle R\rangle)$, (see Supplementary Note 6, Supplementary Note 7, and Supplementary Figure 6 for generalization) instead of two discrete subpopulations, $\langle R_{1}\rangle$ and $\langle R_{2}\rangle$, which of course is much simpler and may suffice in many cases as the results in Fig. 4 illustrate. As a proof of concept of our theory, we solve $P(\langle R\rangle)$ for the eight pairs of contacts analyzed in the previous section. The results are shown in Supplementary Figure 7. In all cases, $P(\langle R\rangle)$s are found to be multi-modal. For peak1/2/3/4-control and peak3-loop, $P(\langle R\rangle)$ yield peaks located at positions very close to $\langle R_{1}\rangle$ and $\langle R_{2}\rangle$ shown in Supplementary Table 1, 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 finn2019extensive , which reports spatial distance measurements for 212 pairs of loci. $P(\langle R\rangle)$ is solved for each of a total of 212 pairs of loci. To illustrate our results, we compare in Fig.5 the predicted $\mathrm{CDF}(r)$ and the experimentally measured $\mathrm{CDF}(r)$, as well as the $P(\langle R\rangle)$ obtained by fitting for six pairs of loci as examples in Fig.5. The results show substantial variations in $\langle R\rangle$, manifested by the multiple peaks and wide spreaded variations in $P(\langle R\rangle)$. Remarkably, the calculated $\mathrm{CDF}(r)$ (without any adjustable parameters) and the measured $\mathrm{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, are extremely small.

In Fig. 6a we show the normalized distributions $P(\langle R\rangle/\mu(\langle R\rangle))$ for each of the 212 pairs of loci (see Supplementary Figure 8 for each pair as a separate figure). We expect that $P(\langle R\rangle/\mu(\langle R\rangle))$ should be narrowly distributed around value 1 if there is only one population. However, many $P(\langle R\rangle/\mu(\langle R\rangle))$ show multiple peaks with large variations. To further quantify the extent of heterogeneity, we calculate the coefficient of variation, $\mathrm{CV}=\sigma(\langle R\rangle)/\mu(\langle R\rangle)$ where $\sigma(\langle R\rangle)$ and $\mu(\langle R\rangle)$ are the standard deviation and the mean of $\langle R\rangle$, respectively. If there is only one population associated with $\langle R\rangle$, $\mathrm{CV}$ should have a value of around zero. Fig. 6b shows the histogram of $\mathrm{CV}$ for all 212 pairs of loci. The $\mathrm{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.

The role of loop extrusion in chromosome heterogeneity. 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). The second is the dynamic heterogeneity. Each cell explores a local minimum of the GFL before transiting to another local minimum (Fig.7b). 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 to the life time of a cell. Such a picture is supported by recent single-cell molecule experiments Hansen2017elife ; Hansen2017nucleus . 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 to the time scale of the cell cycle (15-30 hours). Loop extrusion model alipour2012self ; sanborn2015chromatin ; fudenberg2016formation is another possible origin of dynamic heterogeneity. In the loop extrusion model, it is thought that cohesins extrude loops along the chromsome 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 mutiple subpopulations in chromosomes arise because of massive dynamic heterogeneity. Our theory also gives an indirect theoretical justification for the work in fudenberg2017fish in which the authors found the loop extrusion model could lead to the $[P_{mn},\langle R_{mn}\rangle]$ paradox.

Single-cell temporal information is necessary to determine whether the loops are static or dynamic or a combination of the two (Fig.7c). 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

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}$,$\langle R_{mn}\rangle$] 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 rao20143d . The central result of the theory in Equation (3) can be used to analyze the available sparse FISH data. We show that the fraction of cell subpopulations ($\eta$ in Equation (3)) and the generalization derived in Supplementary Note 6 can be extracted by fitting the FISH data using our theory. From Equation (3) we calculate the Hi-C contact probabilities, thus establishing that the theory resolves the [$P_{mn}$,$\langle R_{mn}\rangle$] 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 multi-point interactions have suggested that formations of loops are likely to be cooperative quinodoz2018higher ; Bintu2018 , 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 Camacho95PNAS . 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.

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

Generalized Rouse Model for Chromosome. 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 bryngelson1996internal ; shukron2017statistics . We first consider a minimal cross-linked phantom chain model, which incorporates the presence of CTCF/cohein mediated loops rao20143d . The model, originally introduced for describing physical gels bryngelson1996internal , and more recently used for chromosome dynamics in a number of insightful studies shukron2017statistics ; shukron2017transient , could be viewed as a Generalized Rouse Chromosome Model (GRMC) hyeon2008force ; o2009accurate . The cross-links modeling the CTCF/cohein mediated loops here are not random. Their locations are predetermined by the Hi-C data rao20143d .

The equations of motion for the GRMC is doi1988theory ,

[TABLE]

where $\xi$ is the friction coefficient, $\boldsymbol{\mathrm{R}}=[\boldsymbol{\mathrm{r}}_{1},\boldsymbol{\mathrm{r}}_{2},...,\boldsymbol{\mathrm{r}}_{N}]^{\mathrm{T}}$ with $\boldsymbol{\mathrm{r}}_{i}$ being the position of the $i^{th}$ locus. The vector $\boldsymbol{\mathrm{F}}=[\boldsymbol{\mathrm{f}}_{1},\boldsymbol{\mathrm{f}}_{2},...,\boldsymbol{\mathrm{f}}_{N}]^{\mathrm{T}}$ ($\mathrm{T}$ is the transpose), where $\boldsymbol{\mathrm{f}}_{i}$ is the Gaussian random force acting on the $i^{th}$ locus, characterized by $\langle f_{n}(t)\rangle$=0 and $\langle f_{n\alpha}(t)f_{m\beta}(t^{\prime})\rangle=2\xi k_{\mathrm{B}}T\delta_{nm}\delta_{\alpha\beta}\delta(t-t^{\prime})$; $\boldsymbol{\mathrm{A}}$ is the $N\times N$ connectivity matrix, embedding the information of chain connectivity and the location of the loops connecting two loci (Fig. 1(a));

[TABLE]

where $\Sigma$ is the set of indices representing the loci pairs specifying the CTCF facilitated loop anchors, and $|\Sigma_{m}|$ is the number of loops connected to the $m^{th}$ locus. The spring constant $\kappa$ enforces chain connectivity, and $\omega$ 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 $\boldsymbol{\mathrm{A}}$ in Equation (5) is a real symmetric matrix, it can be diagonalized using the orthonormal matrix $\boldsymbol{\mathrm{V}}$,

[TABLE]

where $\lambda_{0},\lambda_{1},...,\lambda_{N-1}$ are the eigenvalues of $\boldsymbol{\mathrm{A}}$. By defining $\boldsymbol{\mathrm{X}}=\boldsymbol{\mathrm{V}}\boldsymbol{\mathrm{R}}$ and using $\boldsymbol{\mathrm{R}}=\boldsymbol{\mathrm{V}}^{T}\boldsymbol{\mathrm{X}}$ and $\boldsymbol{\mathrm{V}}\boldsymbol{\mathrm{V}}^{T}=\boldsymbol{\mathrm{I}}$, we obtain the equations of motion of the normal coordinates $\boldsymbol{\mathrm{X}}$,

[TABLE]

Because $\boldsymbol{\mathrm{\Lambda}}$ is a diagonal matrix, the normal coordinates of the GRMC $\boldsymbol{\mathrm{X}}_{p}$ are decoupled. Using the normal modes, $\boldsymbol{\mathrm{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 $\boldsymbol{\mathrm{A}}$, and the orthonormal matrix $\boldsymbol{\mathrm{V}}$ numerically, and thus calculate the contact probability and spatial distance precisely.

Relation between contact probability and mean spatial distance. The vector between the positions of the $m^{th}$ and the $n^{th}$ loci may be written as,

[TABLE]

where $V_{pm}$ and $V_{pn}$ are the elements of orthonormal matrix $\boldsymbol{\mathrm{V}}$. The equilibrium solution of Equation (7) yields, $\lim_{t\to\infty}X_{p,\alpha}(t)\sim\mathcal{N}(0,-\frac{k_{\mathrm{B}}T}{\lambda_{p}})$, where $\alpha=x,y,z$, $\mathcal{N}$ is Gaussian distribution. Therefore,

[TABLE]

where $\sigma_{mn,\alpha}=-\sum_{p=0}^{N-1}(V_{pm}-V_{pn})^{2}(k_{\mathrm{B}}T/\lambda_{p})$. Since the model is isotropic, it follows that $\sigma^{2}_{mn,x}=\sigma^{2}_{mn,y}=\sigma^{2}_{mn,z}\equiv\sigma^{2}_{mn}$. The mean distance $\langle R_{mn}\rangle$ is related to $\sigma_{mn}$ through $\langle R_{mn}\rangle=2\sqrt{2/\pi}\sigma_{mn}$. The distribution of distance between the $m^{th}$ and the $n^{th}$ loci, $\lim_{t\to\infty}|\boldsymbol{\mathrm{R}}_{mn}(t)|=\lim_{t\to\infty}\sqrt{\sum_{\alpha}R^{2}_{mn,\alpha}(t)}$ is a non-central chi distribution (we will neglect the notation $\lim_{t\to\infty}$ from now on),

[TABLE]

The contact probability $P_{mn}$, for a given threshold $r_{\mathrm{c}}$ (contact exists if $r\leq r_{\mathrm{c}}$), computed using Equation 10 yields,

[TABLE]

The mean spatial distance $\langle R_{mn}\rangle$ is given by,

[TABLE]

Using Equations (11) and (12), the desired relation between $P_{mn}$ and $\langle R_{mn}\rangle$ becomes,

[TABLE]

Equation 13 is identical with Equation (1) in the main text.

Simulations. The energy function for the GRMC is,

[TABLE]

For the bonded stretch potential, $U_{i}^{\mathrm{S}}$, we use,

[TABLE]

where $a$ is the equilibrium bond length. The interaction between the loop anchors is also modeled using a harmonic potential,

[TABLE]

where the spring constant is associated with the CTCF facilitated loops, and $\{p,q\}$ represent the indices of the loop anchors, which are taken from the Hi-C data rao20143d (Supplementary Note 3). We simulate the chromosome segment from 146 Mbps to 158 Mbps of Chromosome 5. Each monomer represents 1200 bps, resulting total number of coarse-grained loci $N=10,000$.

In order to accelerate conformational sampling, we perform Langevin Dynamics simulations at low friction honeycutt1992nature . 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 Figure 4).

Data analyses. The contact probability between the $m^{th}$ and $n^{th}$ loci in the simulation is calculated using,

[TABLE]

where $\Theta(\cdot)$ is the Heaviside step function, $r_{\mathrm{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,

[TABLE]

The objective is to go from $P_{mn}$ to $\langle R_{mn}\rangle$, and to determine, if in doing so, we get reasonable results. Because these quantities can be computed precisely in the GRMC, the $[P_{mn},\langle R_{mn}\rangle]$ relationship can be tested, which allows us to obtain the needed cues to solve the FISH-Hi-C paradox.

Data availability. All relevant data supporting the findings of this study are available within the article and its Supplementary Information 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 is available upon request from the authors.

Code availability. The polymer simulations are performed using LAMMPS Molecular Dynamics Simulation software Plimpton1995 , 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.