Measurements of Contact Probability Using Genome-Wide Averages Are Inconsistent with an Ordinary Polymer at Equilibrium.

In our original Hi-C study (5), we characterized the behavior of chromatin at the megabase scale by analyzing the contact probability function, I(s), described above, based on Hi-C data, analytical estimates, and in silico studies. The data for human chromatin showed a power-law relationship of the form I(s) s⁻ between 500 kb and 7 Mb, with = 1.08. We showed that values of can be used to discriminate among distinct polymer states. Specifically, = 1.08 is inconsistent with the classic structure of a globular polymer at equilibrium (an “equilibrium globule,” which has = 1.5), but is consistent with a dense, scale-invariant, isotropic, long-lived polymer state known as the fractal globule (5). Because the fractal globule’s unknotted topology makes it easier to access individual genomic loci, it furnishes an appealing model for chromatin structure.

When we repeated the above analysis on our new, kilobase-resolution maps, we observed a scaling of = 1.27 between 300 kb and 3 Mb (SI Appendix, Fig. S2A). This slightly higher value is consistent with our previous conclusion that chromatin does not fold into an equilibrium globule, and falls within the range we predicted for a fractal globule (17).

Genome-Wide Measurements of Chromatin Folding Inside Individual Contact Domains Reveal a Polymer State Characterized by ≈ 0.75.

In our original study, we could not discern local folding features at scales smaller than ~1 Mb. With our new maps, which contain 200- to 1,000-fold more data, we had the opportunity to study folding within contact domains, which are contiguous genomic intervals in which there is an enhanced probability of contact among all loci (Fig. 2A). We found that folding measurements differ sharply within contact domains vs. across contact domains.

Open in a separate window

Fig. 2.

Contact domains exhibit a contact probability scaling with ≈ 0.75. (A, Left) Contact domains from a region on chromosome 4 of GM12878 lymphoblastoid cells. (A, Right) Number of contacts (Top) incident on a 50-kb window at the center of a domain (Bottom). (B) Contact probability vs. genomic distance for 473 individual domains, measured with respect to a 50-kb locus at the domain’s center. A power law (reference slope of −0.75, gray dashed line) is consistently observed inside domains, whose boundary is indicated by a vertical dashed line. A single black line shows contact probability for the domain from A. Domains are grouped by size; each group is vertically shifted by an order of magnitude for visual clarity. (Inset) Contact probability vs. genomic distance, excluding pairs of loci that lie in different contact domains. A power-law ( = 0.76) is seen over two orders of magnitude. (C) Histogram of values for contact domains across six human cell types. (Inset) Representative microscopy images (maximum Z projections) of four cell types, showing chromatin (blue, DAPI stain) and cytoplasm (red, CellTracker CMTPX dye [Thermo Fisher]). (Scale bar: 10 μm.)

We began by calculating I_same(s) using our genome-wide averaging technique, but only including pairs of loci that were in the same contact domain. We obtained a markedly lower value: = 0.76 (Fig. 2B).

Next, we measured the decline in contact probability with distance relative to a fixed DNA locus. For loci longer than 50 kb, the results were highly reproducible (SI Appendix, Fig. S2G). We focused on 1,057 distinct 50-kb loci, each situated at the midpoint of a contact domain. The resulting contact probability plots consistently exhibited two distinct regimes. Contact frequency within a domain exhibited a power-law decline whose was centered on 0.75 (SD = 0.05) and always smaller than 1 (Fig. 2B and SI Appendix, Fig. S2D). Outside the domain, the power-law regularity was replaced by a more heterogeneous monotonic decline (Fig. 2B and SI Appendix, Fig. S2E).

Our findings suggest that because the frequency of contact between two loci declines markedly when a contact domain boundary is crossed, I(s), which is calculated predominantly using pairs of loci separated by such a boundary, tends to overestimate for contact domains (SI Appendix, Fig. S2B).

To check whether depended on nuclear volume, we compared four human cell types using confocal microscopy and in situ Hi-C. Despite threefold variation in nuclear volume (from smallest to largest, GM12878: 237 ± 84 μm³, IMR90: 381 ± 157 μm³, NHEK: 440 ± 90 μm³, HMEC: 728 ± 307 μm³), the intradomain measurements were indistinguishable (Fig. 2C and SI Appendix, Table S1). The results were also similar in different nuclear compartments (A/B) (5) and subcompartments (A1/A2/B1/B2/B3) (8) (SI Appendix, Fig. S2C), in mouse (CH12-LX lymphoblasts), and in experiments with and without cross-linking (SI Appendix, Table S2).

Because site-directed recombination relies on the spatial proximity of DNA sites, we reexamined published experiments probing the relationship between flippase recombination frequency and genomic distance in humans (18). Interestingly, we noted a power-law scaling with = 0.75 (SI Appendix, Fig. S2F).

Taken together, these results suggest that chromatin in contact domains is characterized by ≈ 0.75. Next, we sought to understand the implications of this exponent.

New Mathematical Theorem Indicates That Chromatin Folding Inside Contact Domains Is Not Strictly Fractal.

A difficulty in interpreting is the uncertainty about which values are consistent with a fractal globule. Approximate methods and physical simulations suggested values ranging from 1 to 1.33 (5, 17). However, no rigorous bounds are known. We therefore sought to derive such bounds.

We proved mathematically that must lie between 1 and 2 for any fractal structure. To prove this result, we analyzed mathematical functions (denoted f) that continuously map (in other words, fold) the unit segment [0,1] into a higher dimensional space. We focused on fractal curves, which are generated by applying a simple folding rule to a simple initial state and repeating ad infinitum. When the folding rule is applied identically at all scales, the resulting fractal curves have no characteristic length scale. Because such curves continuously transform a 1D segment into a higher dimensional object, they have been of interest to mathematicians ever since the first space-filling curves, which map the unit segment onto the unit square, were discovered by Giuseppe Peano (the “Peano curve” in 1890) and David Hilbert (the “Hilbert curve” in 1891). If the repetition process is terminated after a finite number of steps, the resulting curve corresponds to a physically realizable polymer chain. For this reason, finite iterations of fractal curves are often used to model the fractal globule (5). By bounding the values of that can be obtained from fractal curves, we can test whether our experimental data are consistent with a strict fractal globule.

When characterizing a fractal curve, a commonly used measure is the Minkowski (or “box-counting”) dimension, denoted dim(X), which generalizes the common notion of dimension to non-integer values. Just as the number of line segments with width 1/N needed to cover the 1D unit segment scales as N¹ and the number of squares with width 1/N needed to cover the 2D unit square scales as N², dim(X) is defined so that the number of boxes with width 1/N needed to cover X scales as N^dim(X). For instance, the Minkowski dimension of a crumpled sheet of paper (≈2.51) provides a measure of its packing density (19). The Minkowski dimension of Great Britain’s coastline (≈1.25) is a measure of its roughness (20). The Minkowski dimension can also be less than 1: The set of points in [0,1] without an odd digit in their decimal expansion (i.e., 0.86, 0.22222) has dimension 0.699.

We proved that folding the 1D unit segment [0,1] into a d dimensional fractal curve scales the Minkowski dimension of all subsets of [0,1] by a factor of d; that is, any k-dimensional subset of the unit segment will fold into a k ^* d dimensional shape. Our results can be summarized in the following theorem and corollary, whose proofs appear in SI Appendix:

• Theorem: For any self-similar fractal curve f([0,1]), dimf(X) = d ⋅ dimX for any $X\subseteq \left[\mathrm{0,1}\right]$.^*

• Corollary: The contact probability of a fractal curve satisfies I(s) s⁻ with = 2 − (d_surf/d), where s is the linear distance along the curve, d_surf is the Minkowski dimension of the curve’s surface (i.e., its surface roughness), and d is the Minkowski dimension of the curve as a whole.

An illustration of the theorem is the 2D Dragon curve, which doubles the Minkowski dimension of all subsets in its domain (Fig. 3A). Mathematically, our result is a deterministic analog of Henry McKean’s well-known “dimension-doubling” theorem, which states that Brownian motion doubles the dimension of subsets (13).

Open in a separate window

Fig. 3.

New mathematical theorem demonstrates that chromatin folding inside contact domains is not strictly fractal. (A) Successively applying a simple folding rule transforms a 1D line segment (Left) into a 2D Dragon curve (Right). Our theorem reveals that just as the curve doubles the dimension of the line segment, it doubles the dimension of all subsets of the line segment. Thus, when we intersect the Dragon curve with a line to create a 1D feature, the corresponding points in the original segment must have a Minkowski dimension of 1/2. A corollary of this theorem makes it possible to calculate the contact probability scaling exponent for any fractal curve. (B) Contact probability vs. distance for various structures. Our corollary predicts that for fractal curves satisfies = 2 − (d_surf/d), where d_surf is the dimension of the curve’s surface and d is the dimension of its interior. (Left) Comparing simulations (solid lines) with theoretical predictions (dashed lines). Two-dimensional Hilbert curve (purple; d_surf = 1, d = 2, = 1.5), 3D Hilbert curve (blue; d_surf = 2, d = 3, = 1.33), inside-out Hilbert curve, rank 3 (teal; d_surf = 1.5, d = 2, = 1.25), and fractal globule (green). As the rank of an inside-out Hilbert curve increases, its boundary becomes nearly 2D and its draws asymptotically close to 1. (Right, Top to Bottom) Three contact domains: chromosome 12: 46.2–46.4 Mb, chromosome 4: 21.8–22.4 Mb, and chromosome 5: 2.1–4.8 Mb.

The corollary may be illustrated by measuring for classic fractal curves, such as the 2D Hilbert curve (d_surf = 1, d = 2, = 3/2; Fig. 3B, purple), the 3D Hilbert curve (d_surf = 2, d = 3, = 4/3; Fig. 3B, blue), and many others (SI Appendix, Fig. S3 and Table S3). The corollary also implies that for curves with extremely rough surfaces (d_surf close to d), can draw arbitrarily close to unity. Because no such curves are known, we generalized the Hilbert curve, constructing a class of “inside-out” Hilbert curves (Fig. 3B, teal) whose boundaries are arbitrarily rough and whose values come arbitrarily close to 1 (SI Appendix, Figs. S6 and S7).

Because 0 ≤ d_surf/d < 1, the corollary proves that for a fractal curve must lie between 1 and 2. Thus, our measurements of ≈ 0.75 inside contact domains (Fig. 3B, red) are inconsistent with the hypothesis that contact domains tend to form fractal globules.

Physical Simulations Suggest That = 0.75 Is Consistent with an Unknotted, Nonequilibrium State That Is Anisotropic Rather than Fractal.

Another way of exploring the significance of a particular value of is by computationally modeling chromatin as a polymer comprising many identical monomers, each of which represents a fixed number of bases. By simulating the dynamics of a condensing polymer chain and the surrounding mixture under various physical assumptions, it is possible to test whether a particular set of conditions leads to a realistic value.

In our original models, we simulated a condensation process in which the collapse of the polymer was driven by external forces (i.e., the crowding of a stretch of chromatin by other components of the nucleoplasm). Through an excluded volume interaction, these components crush the polymer chain. Such forces can be modeled using a potential function that attracts all monomers equally toward a single point (5). Because this potential has no characteristic length scale, the resulting dynamics are scale-invariant, and the polymers collapse isotropically into a fractal globule.

Notably, our earlier models did not examine the effects of interactions among the monomers on polymer condensation. Attractive forces between individual nucleosomes have been observed in vitro by many groups (2), and effective attractions between monomers are seen in all polymer globules, arising when a polymer is immersed in a poor solvent (21). Therefore, in the present study, we incorporated attractive forces between monomers using the classic Lennard–Jones potential. This model of intermolecular attractions was originally developed to study van der Waals effects and is commonly used to describe the attractive forces between nucleosomes (15, 22). We examined a class of systems in which both intermonomer attractions and external crushing forces are present during condensation. The ratio of these forces is given by a single parameter, R, which we varied across eight orders of magnitude.

We probed these systems using Langevin dynamics simulations, in which random collisions between the solvent and the polymer are accounted for implicitly through parameters for viscosity and temperature. We ran our simulations using LAMMPS (23), accelerated using graphical processing units (24). Each monomer represented 1 kb (matching the above estimates of Kuhn length). The polymer chain contained up to 10,000 monomers, or 10 Mb. For each condition, we ran at least 100 simulations from randomized starting configurations.

Our simulations revealed a family of nonequilibrium states (Fig. 4A). When internal forces are weak (R << 1), polymer condensation closely resembles the isotropic dynamics observed with pure external forces (R = 0) and results in a fractal globule. Because intermonomeric attractions decay as monomers move apart, internal forces introduce a length scale into the system. When they are sufficiently strong (R >> 1), condensation transitions into an anisotropic regime. First, tiny globules form along an extended chain; tension along the chain then causes the globules to concatenate in a linear fashion (Fig. 4B). This model of polymer condensation was first postulated by de Gennes (21). The resulting state, which we dub a “tension globule,” is not scale-invariant. Instead, it contains long intervals in which genomic position is correlated with spatial position along a linear axis (SI Appendix, Fig. S8 D and E).

Open in a separate window

Fig. 4.

Value of obtained using Hi-C is consistent with a tension globule in which loops form by diffusion. (A) Value of for a polymer after condensation varies as the ratio of internal to external forces, R, is changed. (B) Condensation of a 10-Mb tension globule over time. (C) Contact probability vs. distance for 450 simulated tension globules with a length of 10 Mb. In each case, a power law is seen for distances between 20 kb and 800 kb. Mean = 0.73, SD = 0.07. (D) Simulation of a region of chromosome 4 in GM12878 (chromosome 4: 20.3–22.6 Mb). The experimental data exhibit four loop domains (Left) and can be recapitulated (Right) using a tension globule containing four loops (black arcs) which is tethered at both ends (SI Appendix). Contact domains form spontaneously.

Importantly, the values of differ depending on the regime. When R is small, (R) is slightly larger than unity, consistent with our earlier fractal globule simulations. When R is large, (R) is roughly 0.72, consistent with our observations for contact domains (Fig. 4C and SI Appendix, Fig. S8A). In-between, (R) exhibits intermediate values. Interestingly, all of the states in this family are dense and largely unknotted (SI Appendix, Fig. S8F).

Our findings were robust to variations in numerous simulation parameters: chain length and initial configuration, solvent temperature, viscosity, and simulation time (SI Appendix, Fig. S9 and Table S4). They were also robust to the mechanism underlying the internal forces. They did not change significantly when we replaced the Lennard–Jones potential with a Yukawa potential, a model of screened electrostatic forces, in which the attractions decay much more rapidly with distance (SI Appendix, Fig. S10 and Table S4).

Our simulations again suggest that the structure of nuclear chromatin inside contact domains is not consistent with a fractal globule. However, the value is consistent with a tension globule resulting from condensation amid strong internal attractions between monomers.

Condensation of a Tension Globule Results in Spontaneous Formation of Contact Domains Between the Anchors of a Loop.

One of the most surprising features of our in situ Hi-C maps is that contact domains often correspond to loops, that is, the domain boundaries lie at the loop’s two anchor loci.

We used our physical simulations to explore the effects of bringing together two anchor points followed by condensation into a tension globule. Notably, the formation of a loop led to enhanced contact frequency between all pairs of loci in the interval demarcated by the two loop anchors (i.e., a contact domain) (Fig. 4D). These contact domains exhibited values of (0.77 ± 0.08) that match our experimental observations (0.75 ± 0.05; SI Appendix, Fig. S11B).

Simulation Results for a Tension Globule Match Intradomain Distances Measured by 3D-FISH.

We next examined whether the tension globule model recapitulates the spatial distances observed experimentally. We studied four pairs of loci using 3D-FISH (8). Each pair lay in a single domain; the genomic distance between the loci ranged from 320 kb to nearly 1 Mb. We measured at least 50 3D distances for each locus pair. We compared these values with distributions obtained using our tension globule simulations (SI Appendix, Fig. S13B). The simulated distributions matched the experimental distributions as closely [Kolmogorov–Smirnov (K-S) statistic = 0.15] as experimental replicates matched one another (K-S statistic = 0.18).

Taken together, our findings show that the tension globule model accounts for three important features of genome folding observed in our 3D maps: the value of the contact probability scaling , the formation of contact domains between loop anchors, and the distribution of 3D distances between pairs of loci.

Tension Globule Model Fails to Explain Other Key Features of 3D Folding.

Although the tension globule model can explain the general scaling properties of our data, it does not explain many other facets, particularly those facets related to loop formation. In a tension globule model, loop formation would occur via diffusion. Looping proteins (e.g., CTCF) would initially bind to DNA anchor motifs. When diffusion brings two anchors into close spatial proximity, the proteins would dimerize, forming a chromatin loop.

This diffusive process inevitably leads to loops that overlap (i.e., a point in the interior of one loop is anchored to a point outside the loop), thus creating chromatin entanglements. By contrast, our experimental data show that pairs of loops rarely overlap. It is also hard to understand, in a diffusive model, why the CTCF motifs at pairs of loop anchors must lie in the convergent orientation. We therefore considered alternative models.

Loop Formation by Extrusion Complexes Would Explain Key Features of 3D Folding.

Nasmyth proposed a model based on an “extrusion complex” containing two DNA-binding subunits that are physically tethered together (25, 26). This complex is loaded onto chromatin at a single locus. Initially, its subunits are bound to nearby DNA elements, forming a tiny chromatin loop. Next, DNA is extruded through the subunits such that the two tethered subunits move in opposite directions with respect to the genome: one forward and one reverse. As a result, the loop continues to grow, without knotting. Eventually, the extrusion complex dissociates from DNA (Fig. 5A, i–iii).

Open in a separate window

Fig. 5.

Model based on loop extrusion makes it possible to recapitulate Hi-C maps accurately using only CTCF ChIP-Seq results. (A, i and ii) Extrusion complex loads onto the fiber at a random locus, forming an extremely short-range loop. (A, iii) As the two subunits move in opposite directions along the fiber, the loop grows and the extruded fiber forms a domain. (A, iv) When a subunit detects a motif on the appropriate strand, it can stop sliding. Unlike diffusion, extrusion cannot mediate co-location of motifs on different chromosomes. (B) Three-dimensional rendering of a 3-Mb extrusion globule from the ensemble described below. Convergent CTCF anchors (orange spheres) lead to an unknotted loop spanning a compact, spatially segregated contact domain (highlighted in blue). (C) Contact probability vs. distance for 12 domains with a length of 1 Mb, created in silico using loop extrusion, measured from a 100-kb locus at the center of the domain. In each case, a power law is seen for distances between 5 kb and 400 kb. Mean = 0.72, SD = 0.06. (D) We use loop extrusion to model a 2.3-Mb region on chromosome 4 of GM12878. CTCF ChIP-Seq signals are normalized and converted into binding probabilities for the simulated extrusion complex. Each peak is assigned a forward (green) or reverse (red) orientation based on the strand of the underlying CTCF motif. Extrusion simulations yield an ensemble of 3D polymer configurations; contact maps for the simulated ensemble (Top) recapitulate the features observed in our kilobase-resolution Hi-C experiments (Bottom), including the position of domains and loops.

We explored the behavior of extrusion complexes in our simulations as follows. As before, we began with a condensing polymer. Extrusion complexes are bound to the polymer at a density that depends on their concentration and dissociate at a rate that depends on their processivity. The complexes cannot pass through one another; if they collide, one must fall off.

We then added one novel feature to Nasmyth’s model (19), based on our observations about the role of CTCF motifs. We assume that each subunit of the extrusion complex recognizes the presence of a particular motif on a particular DNA strand, such as an appropriately oriented CTCF motif, by binding tightly and halting the extrusion process through the subunit. We therefore designated certain monomers in our simulated polymer as anchors, assigning each a forward or reverse orientation. In our simulations, the forward subunit’s progress may be halted by a forward anchor, but not by a reverse anchor; the reverse subunit may be halted by a reverse anchor, but not by a forward anchor (Fig. 5A, iv).

We began by simulating a condensing polymer containing a pair of convergent anchors 1 Mb apart. When an extrusion complex landed between the anchors, it began extruding a loop until its subunits eventually arrived at the anchor monomers. At this point, the extrusion came to a halt, yielding a “persistent loop” between the anchors (i.e., a loop that was present for a protracted period) (Fig. 5B and SI Appendix, Fig. S12C). This extrusion process did not prevent condensation and resulted in a globular state.

When we examined the resulting contact maps, we made three observations. First, a prominent peak was present between the two anchors, reflecting the formation of a persistent loop. Second, there was enhanced contact frequency between all pairs of loci lying between the two anchors, forming a contact domain. Third, the contact domains exhibited power-law scalings in contact probability with values of (0.72 ± 0.06) that match our experimental observations (0.75 ± 0.05; Fig. 5C). These findings reflect the equilibrium for a condensing polymer immersed in a solvent containing extrusion complexes and did not depend on the condensation forces (external or internal) or the initial condition (SI Appendix, Fig. S12 A and B).

Next, we examined whether the extrusion model recapitulates pairwise spatial distances observed using 3D-FISH. Repeating our previous analysis (SI Appendix, Fig. S13C), we found that the simulated distributions matched the experimental distributions almost as closely (K-S statistic = 0.16) as experimental replicates matched one another (K-S statistic = 0.18).

Thus, the extrusion model, like the tension globule model, is consistent with the contact probability scaling , the formation of contact domains between loop anchors, and the distribution of 3D distances between pairs of loci. The extrusion model explains a much wider range of observations, however. In particular, it provides a natural explanation for why loops overwhelmingly form between convergent CTCF sites. Because loops that occur in the same cell cannot overlap under the extrusion model, the model also explains why we rarely see overlapping loops in our experimental data.

Extrusion Model Can Recapitulate Hi-C Experimental Results, Given the Locations of CTCF Binding.

Next, we explored whether the extrusion model could recapitulate Hi-C experimental results in silico using CTCF ChIP-sequencing (ChIP-Seq) data alone.

We began by simulating the folding of a 2.3-Mb target region on chromosome 4 (20.3–22.6 Mb). Our algorithm created an in silico representation of the region as a uniform polymer, adding forward and reverse anchors at peaks of CTCF binding observed in experimental ChIP-Seq data for the region. The strength of each anchor (the likelihood that the corresponding subunit would halt) reflected the amplitude of the ChIP-Seq peak. Anchor orientation was assigned based on the strand of the CTCF motif associated with the peak. The inputs did not include Hi-C data.

We simulated the behavior of the model polymer in a solvent containing extrusion complexes. Strikingly, the resulting contact matrix closely resembled the contact matrix obtained using Hi-C experiments (Fig. 5D; Pearson’s r = 0.964). Peaks and contact domains were observed in the same positions, and appropriate values were obtained inside each contact domain. The results were similar for other regions (SI Appendix, Fig. S12D).

We also sought to simulate the target region using the tension globule model. As before, we identified peaks in CTCF ChIP-Seq data, and assigned each peak an orientation based on the strand of the CTCF motif associated with the peak. However, to obtain contact matrices that resembled our experimental results, we had to impose a number of ad hoc penalties that do not correspond to any natural processes in 3D diffusion: Loops were only allowed between pairs of convergent peaks, and the likelihood of such a loop depended on (i) the strength of the peaks, (ii) the distance between the peaks, and (iii) the number and strength of intervening CTCF peaks. Even so, the simulations did not provide as good a fit as the extrusion model (SI Appendix, Fig. S11E; Pearson’s r = 0.922).

Network of Loops Contains Hundreds of Isolated Cliques and Is Consistent with a Model in Which Consecutive Loops Can Form Simultaneously by Extrusion.

We then explored higher order relationships among the location of loops by constructing a “loop network” for GM12878 lymphoblastoid cells (8). The network’s nodes are genomic loci containing at least one loop anchor. Edges indicate loops.

We were particularly interested in “isolated cliques.” This network motif comprises a set of n ≥ 3 loci such that any pair in the set is connected by a loop (i.e., the set is a “clique”) but none of the loci are connected by loops with loci outside the set (i.e., the set is “isolated”) (Fig. 6).

Open in a separate window

Fig. 6.

Analysis of loop networks reveals many isolated cliques, consistent with a model in which consecutive loops can form simultaneously by extrusion. Cliques of size three (Left), four (Middle), and five (Right), shown as network representations (Above) and in the Hi-C contact map (Below). Nodes correspond to loop anchor loci; edges and open circles indicate a loop called in (8) (solid green), or using a more relaxed threshold (dashed blue). The five-clique exhibits an additional loop (gray) connecting a clique locus to a locus outside the clique. The loop anchor CTCF motifs are indicated; each middle clique locus contains a pair of CTCF motifs in the divergent orientation.

Isolated cliques were dramatically enriched in the loop network. For instance, in GM12878, we observed 69 isolated cliques with three nodes (sixfold enrichment relative to an ensemble of randomized control networks), 16 with four nodes (28-fold), and 1 with five nodes (161-fold). If we allowed for a small number of loops (≤ N − 2) between loci inside the clique and loci outside, these numbers increase: 145 cliques with three nodes (threefold enrichment), 86 with four nodes (12-fold), 5 with five nodes (14-fold), and 1 with six nodes (41-fold).

Typically, the clique loci were positioned one after another in the genome, with no other loop anchors intervening (in 64% of cases). The first clique locus (in genomic coordinates) typically contained a bound CTCF motif in the forward orientation (97%, fourfold enrichment). The last clique locus contained a bound CTCF motif in the reverse orientation (98%, fourfold enrichment). The middle clique loci typically contained CTCF motifs in both orientations (70%, 8.5-fold enrichment). These findings are consistent with the convergent rule for CTCF looping.

Isolated cliques may reflect simultaneous colocation of all clique loci at a single spatial hub within the same cells, forming a chromatin rosette (9, 27). Alternatively, it may be that some, but not all, clique loops are simultaneous. However, another possibility is that all of the loops in our isolated cliques tend to occur in distinct cells. We cannot distinguish between these alternatives using pairwise contact maps from cell ensembles because such maps do not reflect higher order colocation patterns.

However, we noticed that when a middle locus contained a unique motif in each orientation, the motifs were almost always arranged in a divergent configuration: The first motif was in the reverse orientation, pointing toward the preceding clique locus, and the second motif was in the forward orientation, pointing toward the subsequent clique locus (89% vs. 11% in the convergent orientation). Because the extrusion model predicts that the genomic intervals inside simultaneous chromatin loops cannot overlap, the overwhelming bias toward the divergent configuration would be expected only if consecutive loops in our cliques (e.g., loops A-B and B-C in a clique comprising A, B, and C) form simultaneously by extrusion.

Genome Editing of CTCF Motifs Disrupts Corresponding Loops, Consistent with the Convergent Rule.

We next sought to study the formation of loops experimentally. We used CRISPR/Cas9-based genome editing (14) to modify specific CTCF motifs and explore the effects on loop structure. To avoid allelic heterogeneity, we used HAP1, a human haploid cell line.

We generated an in situ Hi-C map of WT HAP1 cells, with 1.1B reads (SI Appendix, Table S6), in which we annotated 8,334 loops and 4,332 contact domains. We chose to study a target region containing three loci: A (chromosome 8: 133.9 Mb), B (134.2 Mb), and C (134.5 Mb). Each pair of these three loci forms loops with one another (A-B, B-C, and A-C). CTCF sites are present at each locus in accordance with the convergent rule: Locus A has a forward-oriented CTCF motif (dubbed A/Forward), locus B has a reverse-oriented CTCF motif (B/Reverse) followed by a forward-oriented motif (B/Forward), and locus C has a reverse-oriented motif (C/Reverse). All three loops are associated with contact domains.

If our loop anchor motif annotation is accurate, deleting the A/Forward site would disrupt the A-B and A–C loops. To test this hypothesis, we performed genome editing by using two guide RNAs designed to flank either side of the A/Forward motif (SI Appendix, Table S7). We grew a clonal population of cells carrying a 17-bp deletion spanning the A/Forward motif.

To study the effects of this and other genome editing experiments, we developed an inexpensive way to monitor Hi-C results only in the target region by performing HYbrid Capture on the in situ Hi-C library (28). We validated this method, dubbed “Hi-C²,” by applying it to WT HAP1, capturing a 2-Mb region (chromosome 8: 133–135 Mb). The results were equivalent to the results obtained using ordinary in situ Hi-C (SI Appendix, Fig. S14).

When we performed Hi-C² on our A/Forward mutant cell line, we observed that, as predicted, the A-B and A-C loops were disrupted (Fig. 7A).

Open in a separate window

Fig. 7.

Genome editing of CTCF motifs allows reengineering of loops in accordance with the convergent rule; the resulting contact maps can be predicted in silico using extrusion simulations. (A) Results of CRISPR/Cas9-based genome editing experiments at chr8:133.8–134.55 Mb in HAP1 cells. Extrusion simulations (Left) and experimental data (Right) are shown. (A, first row) Contact map for the WT locus, calculated using in silico simulations (Left), closely matches the map observed using Hi-C² experiments (Right). (A, second row) Deletion of A/Forward eliminates the A-B and A-C loops and the contact domain boundary at locus A. The predictions of our in silico simulations (Left) closely match the contact map observed using Hi-C² experiments (Right). All parameters in this and subsequent simulations of mutant regions use exactly the same parameters as the simulations of the corresponding WT contact map. The only difference in the mutant simulation is the modification of the appropriate CTCF-binding site (in this case, deletion of A/Forward). (A, third row) Deletion of B/Reverse eliminates the A-B loop. (A, fourth row) Deletion of B/Forward eliminates the B-C loop. (A, fifth row) Inversion of B/Forward eliminates the B-C loop. (A, sixth row) Simultaneous deletion of B/Reverse and inversion of B/Forward eliminates the B-C loop. (B) Similar series of results for chromosome 1 (180.3–181.3 Mb). Notably, the elimination of one loop anchor motif at the middle locus fails to eliminate either the D-E or E-F contact domain. When both loop anchor motifs are eliminated, both the D-E and E-F contact domains disappear. (C) We disrupted a forward CTCF motif by inserting a single base at chromosome 5: 31,581,788. Two loops are disrupted. The domain boundary moves to a nearby, weak CTCF site. Because the binding at this new site was weaker than the threshold value, this new boundary was not predicted by our extrusion simulations. (D) Our data suggest that the region shown in A is typically found in one of two states in wild-type cells. In the first state, both the A-B and B-C loop domains are present, but the A-C loop domain is absent. In the second, only the A-C loop domain is present. The data suggests a similar decomposition for the region in B. (E) Extrusion can explain the formation of exclusion domains. In this example, an extrusion complex forms a loop between adjacent motifs in the convergent orientation. Downstream, a second CTCF motif in the reverse orientation is unoccupied. Obstructed on both sides, extrusion complexes landing in the interval between the two reverse motifs tend to remain inside the interval. This leads to the formation of a domain.

Next, we examined the effect of disrupting the motifs at the B locus. We tested five predictions of the convergent CTCF rule: (i) Disruption of B/Forward (a 159-bp deletion) will eliminate the B-C loop alone (and not affect the loop between A and B, which, by the convergent rule, must be anchored at B/Reverse); (ii) disruption of B/Reverse (142 bp) will eliminate the A-B loop; (iii) inversion of B/Forward will eliminate the B-C loop; (iv) if B/Reverse is disrupted, inverting B/Forward will eliminate the B-C loop, but restore the A-B loop at a slightly offset anchor position (corresponding to the position of the inverted B/Forward rather than B/Reverse); and (v) inversion of both sites will not disrupt either loop, but will offset the positions of both anchors. The experimental data matched these predictions in every case (Fig. 7A and SI Appendix, Fig. S16).

We next examined a similar region on chromosome 1, spanning D (180.5 Mb), E (180.8 Mb), and F (181.1 Mb). We tested four predictions of the convergent rule: (i) Disruption of E/Reverse (a 16-bp deletion) will eliminate the D-E loop; (ii) disruption of E/Forward (10 bp) will eliminate the E-F loop; (iii) inversion of E/Forward will eliminate the E-F loop; and (iv) simultaneously disrupting E/Forward (7 bp) and E/Reverse (16 bp) will eliminate the D-E and E-F loops, leaving only the D-F loop. The experimental data matched these predictions in every case (Fig. 7B).

Finally, we targeted a similar region on chromosome 5, spanning G (31.6 Mb), H (31.8 Mb), and I (32.1 Mb). We inserted a single base pair into the G/Forward motif (chromosome 5: 31,581,788). Strikingly, both the G-H loop and the G-I loop disappeared, but the H-I loop remained (Fig. 7C).

The convergent rule, combined with our loop anchor motif annotation, correctly predicted the looping pattern in all 13 mutants (Fig. 7 and SI Appendix, Fig. S15). When we examined the effect of removing the corresponding CTCF motifs in our extrusion simulations, the contact maps predicted by the simulations closely matched those contact maps obtained experimentally (Fig. 7). These results show that it is possible to reengineer loops in a targeted fashion.

Our experiments also shed light on which loops in a clique occur in the same cells. For example, if all clique loops occurred simultaneously, then A-B would still be in close proximity even after deleting B/Reverse, because the A-C and B-C loops would sequester A and B near C/Reverse. The fact that we can disrupt the A-B loop without affecting the A-C and B-C loops suggests that the A-C and B-C loops do not occur in the same cells. Similarly, our ability to disrupt the B-C loop alone suggests that the A-C and A-B loops do not occur in the same cells. Taken together, our data are consistent with a model in which consecutive loops (i.e., A-B and B-C) tend to occur simultaneously in some cells, whereas the larger loop (A-C) tends to form in other cells (Fig. 7D).

We thank S. Koch, C. McMullen, A. Tripathy, S. Batra, and J. Onuchic for discussions; S. Knemeyer for figure assistance; N. Tarazi for animations; and I. Dodd for insights on flippase recombination. Horizon Genomics assisted with the generation of several engineered cell lines. This study was supported by National Science Foundation (NSF) Grant PHY-1308264 (to A.L.S.), NSF Grant PHY-1427654, NIH New Innovator Award 1DP2OD008540-01, Cancer Prevention Research Institute of Texas Scholar Award R1304, a McNair Medical Institute Scholar Award, and the President’s Early Career Award in Science and Engineering, and funding from the Welch Foundation, International Business Machines, and Nvidia (to E.L.A.).