Introduction

The three-dimensional organization and dynamics of chromosomes are key regulators of gene expression in eukaryotes. Chromatin is hierarchically structured within the nucleus, from megabase-scale active and inactive (A/B) compartments to sub-megabase topologically associating domains (TADs) [1–6]. These features arise through distinct molecular mechanisms.

A/B compartments are strongly associated with epigenetic states and transcriptional activity, yet their formation appears to involve multiple, possibly redundant, processes that lead to phase separation of the different chromatin types [7, 8]. In both mammals (e.g., Mus musculus) and insects (e.g., Drosophila melanogaster), polymer models that incorporate preferential interactions between regions with similar epigenetic marks can robustly recapitulate compartmentalization patterns [9, 10].

While compartmentalization mechanisms appear broadly conserved, the formation of TADs differs between species. In Drosophila, TADs likely emerge from the same mechanisms that drive compartmentalization [10]. In contrast, mammalian TADs are primarily shaped by loop extrusion [11–13], a process in which Structural Maintenance of Chromosomes (SMC) complexes such as cohesin extrude chromatin loops until halted by convergently oriented CTCF binding sites, where loops are stabilized [14–16]. Experimental and computational studies support the coexistence of loop extrusion and compartmental segregation in mammals, with loop extrusion dominating at sub-megabase scales and compartmentalization prevailing at larger genomic distances [17–19].

In contrast, the evidence for loop extrusion as a mechanism of chromatin organization in Drosophila remains debated [20, 21]. In flies, insulating regions are often found enriched at TAD borders, but they appear to have locus-specific roles [22] and TAD borders are characterized by transitions between epigenetic states rather than by preferentially-bound CTCF [4, 23, 24]. Thus, an alternative mechanistic hypothesis proposes that the formation of TADs in Drosophila is likely driven by specific chromatin interactions based on epigenetic marks [21].

The three-dimensional (3D) organization of the genome plays a critical role in cell-type-specific gene regulation during development [25–27]. TADs often encompass genes and their regulatory elements, such as enhancers, facilitating their spatial proximity to target promoters while insulating them from non-target interactions [28–30]. This view of TADs as regulatory domains has helped explain how boundary disruptions can lead to ectopic gene expression [31–33], though exceptions have been described [34, 35]. To better understand the regulatory role of TADs, single-cell and time-resolved approaches have been instrumental. These studies suggest that TADs and compartments represent statistical rather than fixed structural units of genome organization [36], and that the internal dynamics of gene domains are key to modulating promoter–enhancer communication [37].

A recent study by Brückner and Chen et al. used live imaging in Drosophila embryos to simultaneously track enhancer–promoter (E–P) dynamics and transcriptional activity across genomic distances ranging from s = 58 kb to 3.3 Mb (Fig. 1A) [38]. These data enabled inference of chromatin dynamics, which were then compared to predictions from polymer physics models. As expected, the motion of individual loci followed subdiffusive behavior consistent with the Rouse polymer model, in line with previous observations in other eukaryotes [39–41]. However, the two-locus mean-squared displacement deviated markedly from model predictions. Surprisingly, loci separated by large genomic distances exhibited contact frequencies significantly higher than expected from standard polymer models [38]. These findings suggest that long-range enhancer–promoter interactions may be more prevalent than previously thought, but the underlying mechanisms remain unclear. In this study, we aim to address this gap by exploring potential physical and biological drivers of these unexpected long-range contacts.

Figure 1:

We modeled and simulated the dynamics of a large chromosomal region in silico, incorporating two well-established genome folding mechanisms – compartmental segregation and loop extrusion. We compared our simulation results to available experimental results including: enhancer–promoter (E-P) dynamics and physical distance, as well as Hi-C contact maps. We demonstrate that a simple heteropolymer model can reproduce all experimental data with a good agreement. Combining this model with loop extrusion further refines this agreement, suggesting the possible co-existence of these two mechanisms in Drosophila. Notably, the best agreement between simulations and experiments occurs near the critical coil–globule phase transition, where the chromatin fluctuates between extended and compact configurations [42, 43]. These findings suggest that compartmental segregation operating near this transition provides a mechanistic basis for the observed structure and dynamics of Drosophila chromosomes.

Results

Polymer simulations are widely used alongside experiments to investigate the structure and dynamics of eukaryotic chromosomes [44–47]. Here, we simulated a large region of chromosome 2R in Drosophila (Fig. 1A) to reproduce the experimental observations reported by Brückner et al. Specifically, we implemented three models: an ideal chain, a block copolymer, and a loop-extruded polymer.

We compared simulation outputs to experimental data using three observables: (i) the Hi-C contact frequency map from the same developmental stage (nuclear cycle 14) [48] (Fig. 1B, right), (ii) the mean enhancer–promoter (E–P) physical distance R(s) as a function of genomic separation s (Fig. 1C), and (iii) the two-locus mean-squared displacement MSD₂(t) for varying s (Fig. 1C) [38] (see Methods). The first two observables capture structural properties, while the third reflects chromatin dynamics.

To evaluate model performance, we quantified deviations between simulations and experiments using three metrics: σ_Hi-C for the Hi-C map, σ_str for R(s), and σ_dyn for MSD₂(t) (see Methods). For each model, we identified the parameter set minimizing these deviations and visually assessed the fits to determine consistency with experimental data.

Ideal chain model

We first applied our approach to the simulation of an ideal chain, a basic polymer model with well-established theoretical predictions. A 4.2 Mb segment of chromosome 2R was modeled as a bead-spring polymer without excluded volume (Fig. 1D) and simulated as described in Methods. The observed scaling exponents for R(s) and MSD₂(t) were both close to 0.5 (Fig. 1F,G), consistent with ideal chain behavior [49].

We then scaled both time and spatial units to match experimental data for the three observables and computed the deviations σ between simulation and experiment. The simulated Hi-C map, using a contact threshold of k = 5 bead diameters (simulation spatial units, s.u.), showed a decay in contact frequency with increasing genomic distance s, but lacked the structured patterns observed in the experimental map (σ_Hi-C = 0.982; Fig. 1E). The mean spatial distance R(s) deviated from the experimental trend (σ_str = 0.796; Fig. 1F), and while the fit for MSD₂(t) was reasonable for genomic separations up to 190 kb, it diverged at larger distances (σ_dyn = 3.204; Fig. 1G), highlighting the limitations of the ideal chain in capturing long-range chromatin behavior.

To test whether excluded volume effects could improve the model, we added steric interactions to the ideal chain. However, this modification further increased the discrepancies with experimental data (σ_Hi-C = 0.989, σ_str = 1.029, σ_dyn = 4.863), as confirmed by visual inspection of the fits (Fig. S1A–C). These results indicate that the ideal chain, even with excluded volume, fails to capture the structural and dynamic features observed in Drosophila embryos, motivating the exploration of more complex models.

Loop extrusion model

Extensive simulation-based and experimental studies have shown that loop extrusion by chromatin-associated factors—such as cohesins—affects both chromatin compaction and its diffusive scaling behavior [50–52]. We hypothesized that loop extrusion could promote globular chromatin organization while preserving near-ideal chain dynamics. To test this, we implemented loop extrusion using a well-established two-step framework [13, 53]: first, the positions of loop extruding factors (LEFs) were determined on a one-dimensional (1D) lattice, and then these loop constraints were applied to a three-dimensional (3D) bead-spring polymer simulation (Fig. 2A; see Methods).

Figure 2:

We explored a range of LEF occupancies ⟨Ω⟩ from 0.19 to 1.9 LEFs per 100 kb, consistent with previous estimates [13, 54]. LEF dynamics were governed by four parameters: the loading probability p_L (varied to control ⟨Ω⟩), unloading probability p_U = 0.01, sliding probability p_S = 0.5, and crossover probability p_C = 0.5, which allows LEFs to bypass one another along the chromatin fiber (see Methods for details). By fixing p_U and p_S, we maintained a mean loop size of approximately 100 kb, consistent with in vivo estimates in other organisms [41].

Although TAD boundaries in Drosophila are often associated with insulator proteins [20], there is no direct evidence that these elements block LEFs in vivo. Therefore, we did not impose boundary constraints in our simulations; LEFs were allowed to move freely unless stalled by collisions with other LEFs, with the possibility of crossover. To monitor polymer compaction across simulations with varying LEF occupancy, we measured the mean-squared radius of gyration for all enhancer–promoter segments (Fig. 2B). As expected, decreased with increasing mean LEF occupancy, indicating progressive chromatin compaction.

We then evaluated how each model’s agreement with experimental data changed with LEF occupancy by tracking the three deviation metrics (Fig. 2C). The structural deviation σ_str decreased with increasing ⟨Ω⟩ and reached a minimum at 1.9 LEFs per 100 kb. In contrast, the dynamic deviation σ_dyn was minimized at an intermediate occupancy of ⟨Ω⟩ = 1.14 LEFs per 100 kb. The Hi-C deviation σ_Hi-C remained relatively constant across the tested range. These results suggest that loop extrusion contributes to chromatin compaction and can partially reconcile structural and dynamic features observed in vivo, though no single occupancy value optimally fits all observables.

The Hi-C map obtained from the simulation with the lowest deviation revealed a lack of long-range contacts, in contrast to the experimental map (Fig. 2D). The fit to R(s) was improved relative to the ideal chain, and the experimental scaling exponent was reproduced up to s = 190 kb. However, R(s) increased markedly at larger genomic distances, reflecting the limited influence of loop extrusion at these scales in our simulations (Fig. 2E). Similarly, the model accurately captured the MSD₂(t) curves for s = 58 to 190 kb, but failed to reproduce the rapid plateauing observed experimentally for s = 595 and s = 3327 kb (Fig. 2F), indicating incomplete agreement with the dynamic data.

By fitting the scaling factors to convert simulation units to experimental spatial and temporal units (s.u. and t.u., respectively), we estimated a mean LEF sliding speed of approximately 0.8 kb/s, consistent with in vitro measurements of cohesin velocity [55]. The LEF unloading probability p_U corresponded to a mean residence time of ~ 63 s, which is substantially shorter than experimental estimates of ~ 360 s obtained under CTCF-depleted conditions in mouse embryonic stem cells [41, 56].

Taken together, these results suggest that while loop extrusion contributes to chromatin organization, it is insufficient on its own to account for all structural and dynamic features observed in the experimental data.

Chromosome structure and dynamics driven by compartments

Next, we investigated the effect of compartmental segregation on chromatin structure and dynamics. To identify chromatin compartments, we computed the first principal component (PC1) of the experimental Hi-C map at 16 kb resolution (Fig. 1B, left), assigning each genomic bin to one of two compartments, C0 or C1 (see Methods). This analysis revealed alternating compartment blocks of varying genomic lengths within our region of interest (Fig. 1A).

We then assigned each bead in our polymer model to the compartment corresponding to its genomic position (Fig. 3A, left). To model compartmental interactions, we introduced a distance-dependent attractive potential between beads belonging to the same compartment, following established block copolymer approaches [10, 45] (Fig. 3A, right). For simplicity, we assumed a uniform interaction energy ε for both C0–C0 and C1–C1 bead pairs.

Figure 3:

While systematically varying the strength of the attractive potential ε in our simulations, we observed that each enhancer-promoter segment underwent a characteristic coil to globule transition (i.e., larger to smaller mean , Fig. 3B) [42], consistent with previous studies of self attracting polymers [43, 57]. To identify the value(s) of ε that best matched the experimental data, we computed the deviations σ_str, σ_dyn, and σ_Hi-C across a range of ε values. Both σ_str and σ_dyn reached minima at ε = 0.29, while σ_Hi-C was minimized at a slightly lower value of ε = 0.275 (Fig. 3C). Since the contact threshold used to compute Hi-C maps from simulations is not precisely known and depends on experimental crosslinking conditions, we tested a reduced threshold of k = 4 s.u. (from k = 5 s.u.). This adjustment shifted the minimum of σ_Hi-C to ε = 0.285, aligning more closely with the minima of the other two metrics. Visual inspection of the contact maps revealed a strong compartmentalization at both values ε (Fig. S2A,B). Overall, this model yielded significantly better agreement with experimental data than previous models (σ_Hi-C = 0.552, σ_str = 0.314, σ_dyn = 2.544).

To further validate this improvement, we visualized the fits to the three observables. Since compartment annotations were derived from experimental data, the simulated Hi-C map closely resembled the experimental one, though with generally lower inter-compartment contact frequencies (Fig. 3D). The R(s) curve from simulations followed the experimental trend across all scales (Fig. 3E), although the scaling exponent for s = 58 to 190 kb was approximately 0.5, differing from the experimental value. The fit to MSD₂(t) was notably improved, with the curves for s = 595 and s = 3327 kb correctly transitioning into the plateau, although discrepancies remained at shorter time scales for smaller s values (Fig. 3F).

Given the overall improvement of the compartmental segregation model over previously tested models, we extended our analysis to a fourth observable from Brückner and Chen et al.: the scaling of the relaxation time τ (s) of chromosomal segments. The relaxation time τ (s) characterizes the timescale at which the two-locus mean-squared displacement MSD₂(t) reaches its plateau, and was computed as described in Methods. For an ideal chain, theory predicts τ (s) ~ s², whereas experimental data yielded a scaling exponent of 0.7 ± 0.2 [38]. From our best-fit simulation to dynamics (i.e., the one minimizing σ_dyn), we measured a scaling exponent of τ (s) ~ s^0.88, which falls within the experimental uncertainty (Fig. 3G).

Taken together, our results show that the compartmental segregation model captures multiple experimental observables—including structure, dynamics, and relaxation behavior—more accurately than the ideal chain or loop extrusion models. We therefore conclude that compartmental segregation is a major contributor to the structural and dynamic organization of chromatin in the Drosophila embryo.

Combined compartmental segregation and loop extrusion

Genome compartments emerged as the principal driver of chromatin structure and dynamics in Drosophila. However, discrepancies remained, particularly at shorter genomic distances and time scales. We hypothesized that additional mechanisms could influence chromatin behavior at these scales. Given its dynamic nature and known interplay with compartmental segregation in mammals—both synergistic and competitive [18, 19, 58]—we considered loop extrusion as a complementary mechanism. Notably, our earlier loop extrusion simulations showed good agreement with R(s) and MSD₂(t) at ~100 kb scales, with deviations primarily at larger distances.

We therefore combined compartmental segregation and loop extrusion in a unified bead-spring polymer model with excluded volume (Fig. 4A). To identify optimal parameter combinations, we performed a grid search over a restricted range of values: 0.275 ≤ ε ≤ 0.300 for compartment strength and 0.19 ≤ ⟨Ω⟩ ≤ 0.95 LEFs per 100 kb for loop extrusion. For each parameter pair, we computed the deviations from experimental data for all three observables.

Figure 4:

The fit to Hi-C favored low LEF occupancy and slightly reduced compartment strength relative to the best-fit compartment-only model, with a minimum at ε = 0.280 and ⟨Ω⟩ = 0.19 LEFs/100 kb (Fig. S3A). In contrast, the best fit for R(s) occurred at ε = 0.300 and ⟨Ω⟩ between 0.57 and 0.76 LEFs/100 kb (Fig. S3B). The dynamic observable MSD₂(t) followed a similar trend, with a minimum at ε = 0.290 and ⟨Ω⟩ = 0.76 LEFs/100 kb (Fig. S3C). Each deviation metric showed improved minima compared to the compartment-only model.

To identify a global optimum, we normalized the deviations (see Methods) and visualized the combined deviation landscape across the parameter grid. The global minimum occurred at ε = 0.290 and ⟨Ω⟩ = 0.76 LEFs/100 kb (Fig. 4B), corresponding to the best-fit compartment model supplemented by approximately one freely extruding LEF every 130 kb. We summarized the results in Fig. S4, which shows that this combined model provides the best overall agreement with all three experimental observables.

The fit to Hi-C indicated a much higher strength of compartmentalization compared to the experiment in the best-fit simulation case (Fig. 4C). R(s) showed overall visual improvement and a significant change to the scaling exponent observed up to s = 190 kb (β = 0.41). This change in the exponent marks a reduction from that of the block copolymer best-fit from the previous section and aligns better with experimental observations (Fig. 4D). Lastly, the fit for MSD₂(t) was also in better agreement with experimental data (Fig. 4E). The fit at lower s and t showed improvement, the relaxation of the curves for s = 595 kb and 3.3 Mb into the plateau was reproduced and we observed better global overlap between the two sets of curves (Fig. 4E). Thus, we conclude that the combination of compartmental segregation and loop extrusion displays the best agreement with experimental observations.

Discussion

Here, we leveraged coarse-grained polymer simulations to test three mechanistic models of chromatin organization and reconcile structural and dynamical properties of Drosophila chromosomes. First, we confirmed that the ideal chain model, with or without excluded volume, fails to reproduce experimental observations. We then tested a mechanistic model of loop extrusion by systematically varying the mean occupancy of loop extruding factors (LEFs) while maintaining their processivity. Although this model captured some features of chromatin dynamics, it was partially inconsistent with structural data. A recent study [52] similarly showed that the two-locus dynamics reported by Brückner and Chen et al. could be explained by either fast-sliding LEFs or LEFs with long residence times. However, in the absence of LEF blocking and loop stabilization, loop extrusion alone cannot reproduce the structural features observed in Hi-C maps. Further experimental work is needed to clarify the relationship between cohesins and insulators in Drosophila.

We next tested compartmental segregation using a block copolymer model and established its central role in shaping chromatin structure and dynamics. The best agreement with experimental data was achieved near the coil–globule phase transition, a regime known to produce polymers with globular structure and fast dynamics characterized by the Rouse exponent [42, 43]. This near-critical state has also been observed in other eukaryotic systems using imaging-based analyses [42, 43]. We propose that this criticality facilitates enhancer–promoter interactions by reducing their dependence on genomic distance and enabling rapid spatial exploration within the nucleus.

Other recent studies have also sought to reconcile chromatin structure and dynamics using diverse modeling approaches [52, 59, 60]. For example, Kadam et al. [52] used simulations of a mouse embryonic stem cell locus to show that loop extrusion could explain both Hi-C and live imaging data. Shi et al. [60] used polymer models to generate structural ensembles from contact maps and analyzed them to interpret the results of Brückner and Chen et al. Our work complements these efforts by offering a mechanistic interpretation of Drosophila chromatin organization based on physical modeling.

Finally, our combined model of compartmental segregation and loop extrusion revealed that compartments dominate genome organization at 100 kb to Mb scales, while loop extrusion may act at shorter length scales. Although the role of loop extrusion in Drosophila remains debated [20, 21], our results—together with the conservation of cohesin and CTCF—suggest that chromatin extrusion may indeed occur in flies and contribute to their chromatin architecture.

Conclusion

In this study, we used coarse-grained polymer simulations to systematically evaluate mechanistic models of chromatin folding and reconcile structural and dynamical properties of a Drosophila chromosome. Our results showed that neither an ideal chain nor loop extrusion alone could simultaneously reproduce all experimental observables. In contrast, a block copolymer model incorporating experimentally derived genome compartments, particularly near the coil–globule phase transition, significantly improved agreement with experimental data. This agreement was further enhanced by combining compartmental segregation with loop extrusion.

Drawing on recent developments in polymer physics, we propose that Drosophila chromosomes fold in a near-critical regime, enabling rapid spatial exploration and frequent long-range enhancer–promoter interactions. This criticality may represent a general principle of genome organization. Future work should investigate whether similar near-critical folding regimes operate in mammalian systems and how additional factors—such as boundary elements, epigenetic modifications, and transcriptional activity—modulate chromatin structure and dynamics across scales.

Materials and Methods

Computing observables from the experiments and the simulations

Two of the four observables, i.e., the end-to-end distance R(s) and the two-locus mean-squared displacement MSD₂(t), were computed from the raw trajectory data of Brückner & Chen et al. [61] using the code originally published with the article [62]. From the simulation trajectories comprising the coordinates of nearly 7 × 10⁴ polymer conformations sampled over time, R(s) and MSD₂(t) were calculated for s = 58, 82, 88, 149, 195, 595, 3327 kb using Eq. 1 and Eq. 2 respectively.

The effect of the localization error in [38, 63] was tested on R(s) from the simulation trajectories to ensure that this factor could not affect the fits significantly. Indeed, the effects of localization error were negligible (Fig. S5).

The relaxation time τ (s) was measured from the simulations by finding the intercept between the initial regime and the plateau of each MSD₂(t) curve; the scaling exponent was then estimated by computing the slope of τ(s) plotted in the logarithmic scale.

The experimental Hi-C map (16 kb resolution) was obtained from another study [48] performed during the same developmental stage of the Drosophila embryo as in Brückner & Chen et al. The Hi-C map for each simulation was generated by first coarse-graining polymer conformations from the trajectory to 16 kb resolution, then calculating pairwise distances between beads and determining contacts with a contact threshold of k bead radii. These contact maps were then averaged and normalized using Sequential Component Normalization [64]. The contact threshold was selected to ensure consistently high Spearman’s correlation [65] between simulated and experimental Hi-C maps across all simulations (i.e., for all values of ε; Fig. S2C). This analysis identified a suitable range of thresholds, 5 ≤ k ≤ 8 simulation units (s.u.). We selected k = 5 as a representative value, as it provided strong correlation and corresponded to a reasonable physical distance after conversion to nanometers (see next section).

Fitting observables between simulations and experiments

To fit the experimental R(s) to the simulated one, the curve was translated along the y-axis (i.e., the axis corresponding to length) in the logarithmic scale (Eq. 3, equivalent to scaling by a numerical factor) until the deviation σ_str between the experimental and the simulated curves reached a global minimum. No translation was necessary along the x-axis since the genomic distance s was common to both sets of data. The L-BFGS-B optimization algorithm [66, 67] was used for this fit.

This fit was then used to convert simulation spatial units s.u. (i.e. bead-bead bond length) into real spatial units (i.e. nm) using Eq. 4).

This correspondence allowed us to convert the contact threshold used to compute Hi-C maps in real spatial units. Depending on the simulation case, the contact threshold k = 5 varied roughly between 132 and 308 nm based on Δ_R computed from the corresponding R(s) fits. In the case of the best-fit for the block copolymer, , and the contact threshold was estimated to be ≈ 250 nm.

Fitting MSD₂(t) required a more elaborate procedure, as the simulated and experimental datasets contained different numbers of points. Each simulated MSD₂(t) curve was first fit to a fifth-degree polynomial ℱ_s in logarithmic scale to interpolate the data and enable functional evaluation. Then, using the same optimization algorithm, the experimental data points for each genomic separation s were translated along both space and time axes in log–log space—equivalent to fitting both squared displacement and time simultaneously—until the deviation σ_dyn was minimized (Eq. 5).

As with spatial units, real time units (in seconds, s) can be converted to simulation time units (timestep of the Langevin integrator, t.u.) using the temporal scaling factor Δ_t obtained from the fitting procedure (Eq. 6). For the best-fit case of loop extrusion, we obtained , implying that 1 s = 1278.71 t.u.

Lastly, to estimate the deviation σ_Hi-C between the simulated and the experimental Hi-C map, we compute one minus the Spearman’s correlation (ρ) between the corresponding observed-over-expected maps (see Eq. 7).

where M_x denotes a Hi-C map, SCN is the normalization function [64], and OE the observed-over-expected function, consisting in calculating for each value the log ratio with an expected value corresponding to the diagonal on which this value is located.

To compute a global minimum for all observables, the individual deviations were normalized and combined into a single metric. Specifically, the global deviation was defined as , where each represents the min–max normalized deviation for observable x, computed across all parameter combinations for that simulation case.

Modeling and simulations

The 4.2 Mb stretch of chromosome 2R (chr2R:6,457,453–10,658,620), encompassing all enhancer–promoter segments, was modeled as a fully flexible homopolymer bead-spring chain—i.e., an ideal chain—with a total of 4208 beads, each representing 1 kb of genomic distance. For simplicity, the MS2 and parS insertions used in the experiment were modeled as point insertions along the polymer.

Based on data from [38], promoter loci were assigned to beads 3563, 3439, 3593, 3372, 3695, 2926, and 193. The enhancer was positioned at bead 3505 for downstream promoter loci and at bead 3521 for upstream ones, preserving the experimental enhancer–promoter genomic separations of 58, 82, 88, 149, 190, 595, and 3327 kb.

The conformational dynamics of this ideal chain are described by the Rouse model [49, 68]. In all models except the Rouse model, excluded volume interactions were introduced using a parabolic repulsive potential, as defined in Eq. 8.

where the maximum height of the potential was ρ = 5 kT and its width was r_rep = 1

For all models, the thermal fluctuations of the polymer chain over time were simulated under over-damped Langevin dynamics [69] using the OpenMM-based polychrom package [70, 71]. The simulations of the Rouse model with and without excluded volume were first equilibrated for 10⁷ t.u., then the chain dynamics were recorded for 10⁹ t.u., with one conformation being sampled every 10³ t.u.. For all other simulations, these values were 8 × 10⁶ t.u., at least 5.6 × 10⁸ t.u. and 8 × 10² t.u. respectively. For the simulations of the Rouse model, the only parameters used were the equilibrium bead-bead harmonic bond length (1 s.u.) and its standard deviation (≈ 0.28 s.u.); the parameters for the other models are discussed in the sections below.

Loop extrusion

A common two-step framework was used to simulate loop extrusion on the chromosome [13, 53]. For each case, first, a stochastic simulation of LEF dynamics was performed on a one-dimensional lattice of the same length as the polymer. LEFs could stochastically load onto and unload from the lattice. They were allowed to slide across the lattice, and while doing so, could cross one another or collide and stall. The loading, unloading, sliding, and crossing of LEFs were dictated by a set of loop extrusion parameters (Table 1). No insulators or boundary elements were introduced in any of the simulation cases, and therefore, the LEFs moved freely unless they were stalled by another LEF. During the simulation, the positions of the LEFs on the lattice were recorded over time to keep track of the loops formed and dissolved. Next, during the 3D Langevin dynamics simulation of the polymer, beads corresponding to the LEF positions from the lattice-based simulation trajectory were brought together and anchored by harmonic bonds to mimic the bridging action of LEFs, thus forming loops in the chain. By constantly moving the loop anchors to beads that matched the positions on the lattice-based trajectory, dynamic loop extrusion was simulated on the polymer.

Table 1.

The stochastic model was implemented such that the ratio between the probabilities of LEF loading p_L and unloading p_U determined the steady state mean occupancy of LEFs ⟨Ω⟩ on the lattice (Eq. The steady state was dependent on the size of the lattice and the sliding probability of the LEFs p_S, and was achieved for time t >> 10⁴. Therefore, to vary ⟨Ω⟩, p_L was varied while keeping the sliding probability p_S and p_U constant. p_S = 0.5 and p_U = 0.01 were assigned to attain a mean LEF p_U processivity

Steady state achieved when , therefore

From Eq. 6 and the average time taken for a LEF to slide 1 kb (800 t.u.,1), LEF parameters were estimated in real units. Depending on the simulation case, taking the best-fit case, yielding a mean LEF speed of Similarly, the mean LEF residence time was

Compartmental segregation

The compartments of the chromosomal segment were predicted by computing the first principal component of the genomic distance normalized experimental Hi-C map at 16 kb resolution [48]. A large region of the chromosome harboring the loci of interest but excluding the centromere (chr2R:6,464,000– 25,286,936) was used for the computation. Each 16 kb segment of this stretch was associated with one compartment (names C0 or C1). Each bead in the bead-spring chain with excluded volume was assigned a type based on the compartment predicted, resulting in blocks of alternating C0 and C1 beads. Attractive interactions were introduced only between beads of the same type during Langevin dynamics simulations. The strength of attraction between the beads was controlled by the parameter ε, which corresponds to the depth of the attractive well in a custom parabolic potential (Eq. 10). This potential was devised to induce strong attraction between beads at a short range.

where 0.18 ≤ ε ≤ 0.32 and the width as measured from r = 0 was r_att = 1.8

Supplementary Figures

Figure S1:

Figure S2:

Figure S3:

Figure S4:

Figure S5:

Acknowledgements

We thank Maria Barbi, Timothy Földes, Thomas Gregor, David Lleres, and Marie Schaeffer for critical reading of the manuscript. We are grateful to Alessandro Barducci, Maria Barbi, Timothy Földes, Andrea Parmeggiani, and Antoine Coulon for discussions. Funding from the doctoral school EDCBS2 Montpellier (G.G.). The PCIA computing platform at the MNHN. This project was funded by the European Union’s Horizon 2020 Research and Innovation Program (Grant ID 724429) (M.N.) and the French National Research Agency (ANR-23-CE12-0023-01) (M.N.).