We then investigate whether the predicted patterns can be detected in in vivo Hi-C data during S-phase. Although the presence of mild fountain-like patterns around origins in yeast was already reported in our previous study [25] using publicly available data in WT cells [52, 59, 61], it was not possible to directly establish a causal-effect relationship between these observations on Hi-C maps and the replication process, as well as to rule out possible confounding effects such as the concomitant presence of extruding cohesins.

For this reason, we generate new Hi-C datasets during early S-phase in WT and mutant yeast strains to isolate the contribution of replication and decouple the emerging signal from concurrent biological processes (see Methods, Additional file 1: Fig. S7 and S8). For this, cells were arrested in G1 using alpha-factor and then released in S-phase for 20 min before fixation.

In WT condition, analysis of the Hi-C data in early S-phase, shows the presence of a fountain-like pattern around early origins (Fig. 3A Bottom Left) with an enrichment around 10-30 kb. This structure appears specifically during S-phase as it is not visible in G1 (Fig. 3A Top Left) or G2/M-arrested (cdc20, Fig. 3A Top center) cells. Remarkably, the observed signal at this early-S time-point appears significantly stronger than our previous report [25] using publicly-available Micro-C data [61].

Fig. 3

Comparison with in vivo data and forks spatial distribution. A Average normalized (log₂ Observed over Expected) contact maps around early replicating origins for in vivo Hi-C data for different conditions. B log₂ ratio between Scc1-depleted and G1 Hi-C maps for chromosome 7. In this case, maps were plotted at 3 kb resolution showing the full chromosome (Left) and a 400 kb region (Right). Yellow and black dots indicate early and late replicating origins respectively. C Comparison between experimental aggregates of Scc1-depleted cells and simulation using the best matching time-interval from Fig. 2F (from t = 5 min to t = 7 min and 30, see Additional file 1: Fig. S9A). Aggregates were computed around all early, only non-centromeric (centromeres at a distance > 50 kb) and all late origins for the two scenarios of sister-forks organization. Black arrows highlight the qualitative differences between the two simulated signals with a depletion and enrichment at the loop basis in the interacting and non-interacting scenario respectively. The numbers of biological replicates of experimental Hi-C data shown merged in (A-C) are: G1-WT n = 2, early S-WT n = 3, early S Cdc45-depleted n = 2, early S Scc1-depleted n = 1, early S-Ctf4∆ n = 1, metaphase-arrested (cdc20) n = 2

To investigate the role of replication in the formation of the fountains, we perform similar experiments but in Cdc45-depleted cells using an auxin-degron system. Cdc45 being a component of the replicative helicase CMG, it is essential for replication bubble opening at origins [78]. Cdc45-depleted cells thus progress towards G2/M but without performing DNA replication and in the absence of SC cohesion. In particular, it was observed that, in this case, chromosomes still engage in a mitotic-like compaction via cohesin-mediated loop extrusion activity [56, 58].

Measuring the Hi-C maps, as in the WT, 20 min after release from G1 and depletion in Cdc45 (Additional file 1: Fig. S7C), we observe complete loss of the fountain signal around early origins (Fig. 3A Top Right) after Cdc45 depletion.

Overall, this result demonstrates that the progression along the cell cycle (Cdc45-depleted) nor the presence of a second SC (cdc20) are sufficient to retain the fountain signal, which requires ongoing replication and is observed only transiently during S-phase.

Despite such findings, it is still unclear whether the emerging pattern is a direct consequence of the 3D organization of sister-forks, or driven by other concurring processes also happening during S-phase. In particular, in budding yeast, cohesin-mediated SCs cohesion and intrachromatid loop-extrusion constitute two well-established processes that impact the 3D chromosome organization and that initiate during S-phase [49, 50, 51, 52, 53, 54–55]. In particular, cohesins start extruding during early S leading to the formation of loops between specific genomic regions called cohesin associated regions (CARs) [61] (Additional file 1: Fig. S9B).

The early-S Cdc45-depleted and metaphase-arrested maps presented in Fig. 3A demonstrate that the fountains observed in WT early-S are not just a simple by-product of some trivial correlations between CARs and early origins positioning along the genome, as a complete loss of fountains occurred in these cases while chromatin loops at CARs are still preserved (Additional file 1: Fig. S9B). However, this does not exclude that some cohesins might also be specifically recruited at replication forks and contribute to the signal.

To explore this hypothesis, we analyze early-S phase Hi-C map from cells partially depleted of Scc1, a cohesin subunit, using an auxin-degron system (Additional file 1: Fig. S7D). We first verify that such a depletion reduces strongly the cohesin loop extruding activity and removes the loops between CARs in the Hi-C map, while their signature remains present in WT, Cdc45-depleted cells (Additional file 1: Fig. S9B). Around early origins on the contrary, fountains are persistent to cohesin depletion (Fig. 3A Bottom Middle). Unexpectedly, the average fountain signal around origins appears enhanced in the absence of cohesin compared to WT. We speculate that the average decrease in contact frequency (P(s) curve) due to the depletion of loops on chromatin, may play a role when normalizing the maps by the expected value (see Methods). As a result, replication-driven fountains are likely enhanced due to less interplay with loops occurring at similar genomic scale (10s of kb). Interestingly, due to better contrast, fountains can be visualized not only at the population level, but also at the single origin level when we normalize the Hi-C map for Scc1-depleted cells in early-S by the one in G1. In this case, most early origins exhibit a fountain-like pattern (Fig. 3B, Additional file 1: Fig. S10), that are absent when repeating the analysis using metaphase-arrested cells (cdc20, Additional file 1: Fig. S11).

Overall, this suggests that the fountain signal does not depend on cohesin loading and loop extrusion activity.

Given experimental evidence of its ability to oligomerize in vitro [79], Ctf4, a component of the replisome, has been proposed as a possible mediator of sister-fork interactions. We thus perform Hi-C experiments as before but in Ctf4-KO cells that are still progressing through S-phase (Additional file 1: Fig. S7E). In this case, we still observe loops around CARs (Additional file 1: Fig. S9B) and more importantly, the characteristic enrichment around early origins is still present (Fig. 3A, Bottom Left), suggesting that Ctf4 might not be a driver of the fountain pattern. Interestingly, Ctf4 depletion leads on average to smaller replicon sizes (Additional file 1: Fig. S12). We hypothesize that a decrease in the fork progression speed may thus explain the differences in size and strength between fountains in WT and Ctf4-KO.

The Scc1-depleted maps discussed above offer the best comparison to our polymer model, which does not account for cohesin-mediated processes such as cohesion or loop-extrusion and only implements the formation and growth of replication bubbles.

To make the most meaningful comparison possible, we first have to infer the effective time in our simulations corresponding to the early-S stage in the experiments where cells do not start replication immediately after release. To do so, we estimate and compare the simulated and experimental average number of Hi-C reads around early origins, which reflects the local probability of having been replicated (see Additional file 1: Fig. S9A for a more detailed description of the method used). We find a corresponding time interval in our simulations between 5 and 7.5 min (Fig. 2E, top Right).

We then compare the contact signals observed around early origins (Fig. 3C). Overall (Top panels), experimental data are qualitatively more similar to the fountains observed in the interacting fork scenario with common characteristic features: (1) a well-localized enrichment perpendicular to the main diagonal, corresponding to an average replicon size of 30 kb, and (2) a depletion of the signal at the loop basis (see black arrows in Fig. 3C Top).

Since early origins are more densely located around centromeres, to avoid any possible structural artifacts related to the contact enrichment observed around centromeres (Fig. 1D), we also compare specifically the signal for noncentromeric, early origins (Fig. 3C Middle). They are characterized on average by greater inter-origin distances and thus are likely to have experienced fewer merging events between nearby replicons in early-S. In vivo, we still observed a fountain-like pattern which also qualitatively resembles the predictions of the interacting fork scenario.

Overall, the interacting fork predictions are thus more compatible with the fountain shape observed in the experimental data.

However, the strength of the enrichment signal is much stronger in the predictions than in the experiments (see also Additional file 1: Fig. S9, S13). This may suggest an intermediate scenario where not all the sister-forks are paired. Indeed, our working hypothesis of strictly-bound sister forks is likely to be a strong approximation of the real biological system, as it has been experimentally observed that replication with independently moving forks is still viable [26].

We address this hypothesis by mixing predicted contact maps of the interacting and non-interacting fork cases at different ratios (Additional file 1: Fig. S9F). Interestingly, a small percentage of interacting forks signal (≥ 20%) is already sufficient to retain the localized enrichment perpendicular to the main diagonal (fountain pattern) while lowering the absolute value of the enrichment.

We also explore alternative confounding factors that may affect the predicted fountain strengths. While our model captures the experimental 1D dynamics and hierarchical firing of origins, all of the simulated trajectories are perfectly synchronized and enter replication at the same time. This is likely not the case the experiments in which the beginning of the S-phase in cycling cells or in cells released from G1 after synchronization is likely to be more heterogeneous [57, 61]. If we mix predictions for early-S with some for G1 (Additional file 1: Fig. S9C), we also observe a reduction of the signal around early origins. Interestingly, the fountain pattern is very robust even at strong dilution, suggesting that such replication-dependent signals may also be captured in a non-ideally synchronized system.

Finally, we investigate the impact of the model parameters used in estimating the in silico contact maps: the radius of capture (Additional file 1: Fig. S9E) which is the maximal distance between two monomers used to define a contact; and the length of the time intervals (Additional file 1: Fig. S9D) used to compute maps. While these parameters have a quantitative impact on the enrichment, they both maintain the overall qualitative difference between the two fork scenarios.

Overall, experimental data around early origins are consistent with the presence of a significant population of interacting sister-forks in vivo.

We then investigate if our model also well captures the signal around late replicating origins (Fig. 2B, Top). In the in vivo Scc1-depleted maps, we record only a faint signal around late origins, both in terms of size and strength (Fig. 3C Bottom, lower triangle). Predictions in both sister-fork scenarios show strong differences. In the interacting fork case, we continue to record a mild fountain-like pattern that is not present in the experiment (Fig. 3B Bottom Right). On the other hand, in the non-interacting case (Fig. 3C Bottom Left), only a faint signal is present, resembling more the experimental observation. Mixed scenarios with a majority of non-interacting cells are also consistent with the experiments (Additional file 1: Fig. S14). This is consistent with what is observed locally at late origins in the differential plots (Scc1-depleted over G1) of Fig. 3B and Additional file 1: Fig. S10 (black points). While we do observe enrichment around some of the late origins, in most cases no specific pattern is detected. Overall, this suggests that interactions between sister-forks is less likely when fired from a late origin.

In the previous section, we observed a complex and heterogeneous distribution of replication forks within the nucleus driven in part by the Rabl architecture. Here, we investigate whether this peculiar spatial organization can explain by itself the presence of replication foci (RFi), where forks colocalize, as observed experimentally [20].

To define and detect RFi in each simulation snapshot, we consider the independent components of the corresponding adjacency matrix A_i,j, where A_i,j = 1 if forks i and j are distant by less than a threshold distance θ (A_i,j = 0 otherwise) (see Methods). θ may represent the typical spatial resolution of any experiments that allow to “visualize” the replication fork. For example, for fluorescence-based microscopy, we might expect θ to vary from 60 nm (super-resolution-like) to 200 nm (confocal-resolution-like). For both scenarios of sister-forks association, we computed the distribution of RFi during early S-phase, considering all frames between 2.5 and 5 min (Fig. 5A,B) for several θ values.

Fig. 5

Detection of RFi within the nucleus during S-phase. A Distributions of the number of distinct RFi detected at a given resolution θ in early-S simulations (between t = 2 min and 30 s and t = 3 min and 45 s) for non-interacting (Top) and interacting (Bottom) sister-forks. The dashed lines correspond to the mean value of each corresponding distribution. B Percentage of clusters containing a given number of forks for simulations in early-S (between t = 2 min and 30 s and t = 3 min and 45 s) for non-interacting (Top) and interacting (Bottom) sister-forks and for various resolutions θ. The inside panels show the distributions in logarithmic scale, including the probability for higher number of forks per cluster. C Percentage of clusters containing a given number of forks for different time windows in the non-interacting (Top) and interacting (Bottom) sister-fork cases. The inside panels show the distributions in logarithmic scale. D Left: percentage of clusters containing a given number of forks pairs. Gray squares indicate the experimental estimates by Saner et al. (see Methods). Simulated data in early-S (between t = 2 min and 30 s and t = 3 min and 45 s) was analyzed to strictly obtain even number of forks per cluster (see Methods). Right: schemes and color relative to the three models explored in panel D

In the non-interacting scenario (Fig. 5A, Top), at high resolution (θ = 60 nm), most of the individual forks can be resolved. Consequently, the average total number of detected RFi in the system is approximately equal to the average total number of replicating forks at that corresponding moment of S-phase. In the interacting sister-fork case (Fig. 5A, Bottom), as expected, at θ = 60 nm, the total number of RFis now corresponds mainly to the number of replisomes. As θ increases and the resolution is reduced, additional links are incorporated into the adjacency matrix A_i,j, leading to a progressive shift towards lower values for the distribution of the number of RFi (Fig. 5A). Interestingly, for low θ values (≤ 160 nm), the interacting scenario always results in fewer detected RFi, each containing more forks on average due to the forced colocalization of sister-forks. However, for low resolution (θ ≥ 160 nm), observed RFi in the non-interacting case are less numerous (Fig. 5A) and larger on average. This is likely due to the more homogeneous distribution of forks in the non-interacting scenario (Fig. 4).

To further characterize the cluster population, we compute the probability distribution of the RFi sizes (Fig. 5B) at different resolutions θ. In the non-interacting case and at lower resolutions (θ < 125 nm), most detected RFi consist of a single diffusing fork (Fig. 5B, Top). Increasing θ leads to the detection of some larger clusters. In particular, the probability of detecting sister-forks within the same cluster is enhanced, thus biasing the distribution to even numbers of forks (Fig. 5B, Top, lighter green points). In the interacting forks scenario, such a bias to even number is already present at high resolution (Fig. 5B Bottom, darker green points) and the presence of clusters containing 4, 6 or 8 forks is not rare. Since we do not introduce interactions between the two remaining diverging forks after a termination event (see Methods), even in presence of sister-forks interactions, odd numbers of forks can be detected in our simulations (Fig. 5C Bottom) Interestingly, at very low resolution (θ ≥ 160 nm), in both scenarios, very large RFis containing a large fraction of the active forks may be observed (secondary peaks in the inset panels Fig. 5B).

Fixing θ = 125 nm to mirror the optical resolution in the experimental work by Saner et al. [20] (see below), we analyze how the sizes of RFis change during S-phase. In the non-interacting case, at very early S, the distribution shows again a bias in favor of even values (Fig. 5C Top panels, dark violet points) due to the small average replicon size and the subsequent proximity of sister-forks that are detected colocalized at such a resolution. The high density of forks in centromeric regions (as shown in Fig. 3C,D) also increases the probability of observing clusters (over 40% with 4 forks or more). This is rapidly lost as replication proceeds, leading to the detection of a high percentage of single-forks (> 60%), and to the gradual diminution of large RFis (Fig. 5C top panels). In the presence of sister-fork interactions, the bias toward even cluster sizes persists until 5-7 min (Fig. 5C Bottom panels) after the beginning of S-phase. After, termination events where nearby replication bubbles merge lead to the loss of sister-fork interactions and thus to an increasing number of RFi with single fork (Fig. 5C bottom panels, yellow and green points).

In summary, our in silico predictions indicate that random collisions of diffusing forks or sister-fork pairs may indeed lead to the significant detection of RFi containing several forks or replisomes. Clusters of large sizes are more abundantly detected in early-S when the spatial distribution of forks is enriched around the SPB (Fig. 4).

We test whether the predicted, resolution-dependent colocalization of multiple freely diffusing forks is consistent with experimental observations made by Saner et al. [20]. Using super-resolution microscopy (SIM technology) on GFP-tagged PCNA in early S, a key component of replisomes, the authors identified RFi at different time points along S-phase. For one time point where the total number of RFi is maximal, they inferred the number of sister-fork pairs inside a given RFi based on its overall fluorescence intensity (gray squares in Fig. 5D). Note that the inferred distribution is based on an in silico estimation of the total number of forks [80] and—motivated by their previous live-cell imaging observation [19]—on the assumption that sister-forks always co-localize [19] (see Methods). Such a minimal approach does not account that termination events (encounters of convergent forks) can occur even in early S-phase. Aware of such limitations in the experimental estimates, we mirror the data analysis performed by Saner et al. by strictly casting the numbers of forks within RFi to even numbers (see Methods).

In Fig. 5D, we estimate the distributions of the number of fork pairs inside our detected RFis (for θ = 125 nm, the typical spatial resolution of 3D-SIM [81]) at a time (2-4 min after the start of S-phase simulations) where the total number of replicating forks in silico is maximal (Additional file 1: Fig. S5) and comparable to the estimation made by Saner et al. at the time of their measurements (∼ 240 excluding rDNA-related forks). Interestingly, in our simulations, under these conditions, the average number of detected RFis is ∼ 80 (respectively ∼ 100) in the interacting sister-fork scenario (resp. non-interacting scenario), very close to the experimental observation (∼ 75). Both scenarios of sister-forks association exhibit similar trends with a peak at 1 pair per cluster, as observed experimentally, and follow by a rapid decrease for higher numbers of pairs, failing to capture the significant probability of having RFi of intermediate sizes (2-4 pairs) inferred by Saner et al.

We thus wonder if such an underestimation of detected clusters of intermediate sizes in silico may be explained by the presence in vivo of aggregating forces that might stabilize larger clusters. Therefore, in addition to sister-forks interactions, we included non-specific pair-wise contact interactions between all the forks (see Methods, Additional file 1: Fig. S16). At the Hi-C level, we predict an enhancement of contacts between early origins (Additional file 1: Fig. S17) as the strength of non-specific interactions J_f increases, which becomes significant for J_f ≥ 1.75 kT. At the level of the distribution of RFi sizes (Additional file 1: Fig. S18), we observe only an enrichment of large clusters (5-15 pairs) for J_f = 1.75 kT, and intermediate sizes (3-4 pairs) become more frequent for stronger J_f.

Experimental Hi-C data in early S (Additional file 1: Fig. S17D) do not exhibit significant contacts between early origins, suggesting that J_f < 1.75 kT. However, at such interaction strength, the predicted frequency of intermediate RFi sizes is not compatible with microscopy data (Fig. 5D, Additional file 1: Fig. S18).

In summary, the colocalization of replicating forks into RFi predicted by our modeling, either just by random collisions or mediated by aggregating non-specific forces, are not compatible with the super-resolution microscopy data of Saner et al. (see Discussion and Additional file 1: Fig.S19).

Similarly to previous analogous implementations [64, 66, 67–68], we model the nuclear envelope by introducing a spherical barrier to monomers in the system. In the lattice framework used here, this is accomplished by occupying all sites above a certain radius from the box center, now inaccessible to the polymer chains due to excluded volume interactions (Additional file 1: Fig. S1A). Based on microscopy experiments [41, 108], we set the diameter of the sphere to 2 µm which corresponds to ≈ 70 lattice units given a bead size of 20 nm (Additional file 1: Fig. S1B). We then introduce several distinct polymer chains into the system to model each individual chromosome. Since we aim to model the haploid budding yeast genome, the number of chains is set to 17 where the 12th and 13th chains both model chromosome 12 (see below). Each chain is built using the HedgeHog algorithm [103, 109]: the chain is recursively grown to the desired size (different for each chromosome), starting from a simple backbone of length L/4 where L = 70 is the box size. Arbitrarily, we choose the backbone to be divided into two equally long linear branches (each L/8), in random directions (See Additional file 1: Fig. S1A for an example). The starting position of each backbone (corresponding to the first monomer of the chain) is placed at one of the 17 points inside the equatorial plane. These points are produced using the sunflower algorithm [110] to efficiently pack n points inside a circle (Additional file 1: Fig. S1A, small panel). To randomize the process and avoid any bias, chromosomes are built in a random order, ensuring that the two random directions used to build the backbone of a new chain do not cross the ones already initialized. Note that such initial configuration does not aim to mimic any biological scenario (such as chromatids organization upon Mitotic exit). The choice of a such randomized V-shaped backbone allows for simple and efficient construction of non-crossing chromosomes in the lattice, whose chain ends (telomeres) are randomly oriented. Given the strong forces that will be introduced into the system to recreate the Rabl configuration, this initial choice is likely to have a minimal impact. In fact, different initialization strategies in other studies [70] eventually converge to similar results once those forces are established.

To mimic centromere attachment to the SPB, we introduce a spring-like potential between the centromere of each chromosome and position in XYZ coordinates (Additional file 1: Fig. S1B). This forces the centromeres in the first times of the simulation to move toward the SPB (spring constant k = 100 kT). We complement our potential to allow the centromeres to freely diffuse in a spherical shell located between 250 and 325 nm from the SPB. The Hamiltonian describing the kinetics of centromeric monomers is as follows:

2

To simulate the tethering of telomeres to the nuclear envelope (NE), we implement an outward force acting on the terminal monomers of the chains, with a negative spring potential connecting them to the centers of the spheres to push them towards the NE. Again, we allow the telomeres to freely diffuse within a distance of 50 nm from the spherical wall (Additional file 1: Fig. S1B). Therefore, the Hamiltonian describing the kinetics of telomeric monomers is as follows:3Since the nucleolus consists of a large physical barrier for non-rDNA, we introduce an additional impermeable wall positioned at 800 nm in the z-direction (Additional file 1: Fig. S1B,C). Spherical confinement and the nucleolus wall lead to a final volumic fraction Φ = 3%, compatible with the chromatin concentration standardly used for modeling yeast chromosomes [103].

To simulate chromosome 12 and account for its role in the nucleolus assembly, we split chromosome 12 into two polymer chains anchored by one side to the nucleolus wall (Additional file 1: Fig. S1C). The first chain models the first 460 kb of chromosome 12 where its rightmost monomer corresponds to a boundary of rDNA (Additional file 1: Fig. S1C). The remaining portion of the chromosome, without centromere, is modeled with a second chain whose leftmost monomer constitutes the second boundary of rDNA, while the rightmost one is the second telomere. The two rDNA beads are pulled to the nucleolus wall with a spring potential only in the z-direction (downward) enabling their localization at the nucleolus wall periphery:

4

Similar to other implementations [68, 69], no constraint was applied to the 3D distances between the two rDNA boundaries, which are allowed to diffuse independently on the nucleolus surface.

We model explicitly replication in 3D using a specialized class of self-replicating polymers to simulate each chromosome. On-lattice chain duplication in KMC simulations was physically characterized in detail in our previous study [25]. Although, in this study, we include multiple chains and the Rabl-like organization, the local realization of replication is analogous since each polymer chain, while governed by a global underlying 1D replication dynamics (see below), duplicates “independently” in 3D. Briefly, a given unreplicated monomer can be used as an origin of replication, leading to the formation of a replication bubble. In particular, origin firing consists of introducing a new monomer into the system, connecting it to the two neighbors of an origin along the maternal chains. This process effectively creates two polymer branches that locally simulate the two nascent SCs. Monomers with triple connectivity at the extremities of the replication bubbles are defined as replication forks with a given directionality based on their position with respect to the origin. With a constant speed v = 2.2 kb/min, every fork in the system can trigger further replication by introducing additional monomers at the fork positions and moving the triple connectivity. In particular, we convert this physiological speed into lattice units which gives a probability of ∼ 3 · 10⁻⁶ for a given fork to be replicated per MCS. Both in the case of origin firing and fork progression, the newly added monomers are positioned at the same site of their homologous beads, thus temporarily breaking the excluded volume rule. As monomers diffuse rapidly in 3D, such a constraint is quickly reinstated. A more detailed description of the local realization of monomer duplication in the lattice can be found in Fig. S1 of our previous study [25]. Replication bubbles merge once two forks of opposite directionalities converge along the chain. In such a configuration, the two converging forks, now neighbors along the chain, are always replicated in pairs, leading to the formation of a single larger bubble. Similarly, when a fork is adjacent to a telomere, the bubble is opened by replicating both the fork and the chain-end monomers.

We introduce implicit 1D replication dynamics by adapting the formalism of Arbona et al. [62, 63] to our KMC framework. Replication is modeled as a bimolecular reaction between a limited number of N_f firing factors and the assigned origins of replication and a firing trial move is introduced in our scheme. Under the assumption of a well-mixed system, at a given time, a set of N_origin origins is randomly selected (see below) and each origin of this set is fired with a probability:5

where (N_f −N_forks(t)/2) represents the current number of available freely diffusing firing factors and OriginRate is a constant rate, equivalent of k_on in [62, 63]. To model the heterogeneous time patterns of eukaryotic replication, the set of N_origins is sampled from the “Initiation Probability Landscape Signal” (IPLS) computed in [63]. The IPLS is a signal that defines origin priming efficiency and was inferred in [63] through a machine learning approach to accurately reproduce experimental Mean Replication Timing and Replication Forks Directionality data [44]. Practically, the sampling of origins is performed by randomly picking a fixed number of monomers (N_origins = L_chromosome/5 as inferred by Arbona et al. in [63]), the weighted probability of each monomer being given by the binned IPLS at the same resolution of our model (1 kb). Note that the drawing is performed with replacement, which means that each monomer can be selected multiple times to account for that a certain coarse-grained genomic position may contain multiple (fine-scale) origins, all of which have the same probability of firing (OriginRate).We then set N_f= 120 and OriginRate = 1 · 10⁻⁹MC⁻¹ which are the KMC parameters equivalent to the ones inferred by Arbona et al. in [63]. Firing factors are rapidly introduced in the system once the S-phase starts according to the formula N_f(1 − 0.13exp^−t/τ) with τ = 1 min and 30 s.

Co-localization of sister-forks is maintained by adding an elastic-like potential between the forks to the Hamiltonian of the system on which the Metropolis criterion is applied:6

With

7

where the sum is performed over all the pairs of still active sister forks localized at monomers i and j, d_i,j is the 3D Euclidean distance between i and j, Θ(x) is the Heaviside step function and th = 50 nm is the distance below which the energy saturates to J_sister = 100 kT. Springs between forks are kept until one of the two interacting partners is lost due to the merging of replication bubbles or the replication of the chain-ends. In other words, the two remaining forks after the encounter freely diffuse in 3D any binding energy connecting them.

We include contact interactions with strength J_f between all replication forks in the system. Similarly to previous copolymer lattice models [102, 103–104, 106], we include an additional term on the Hamiltonian of the system:8

where σ_i = 1 if monomer i is a replication fork, otherwise σ_i = 0. f_i,j = 1 if monomers i and j are first or second nearest-neighbors otherwise f_i,j = 0. In practice, during the KMC algorithm, we keep track for each lattice site on the number of forks within a 40 nm distance. This number is locally updated after each trial move (see scheme in Additional file 1: Fig. S15). This strategy allows for a rapid scanning for nearby forks and facilitates Hamiltonian evaluation. Similarly to other implementations [102, 103–104, 106], we do not impose any constraints on the number of interactions for a given fork.

The value of J_f = 1.7 kT was established by visually inspecting Hi-C maps (Additional file 1: Fig. S17A) and offdiagonal aggregate plots (Additional file 1: Fig. S17B) between early origins for increasing values of J_f. In particular, we selected an energy value that lies in the transition between non-aggregated and aggregated states [102, 104].

In Fig. 1B, we illustrate the main features of the Rabl organization at the full genome (16-kb resolution, top panel) and single chromosome (1-kb resolution, Bottom panel) levels. To plot together in vivo and in silico maps, we rescale the simulated contact matrices by multiplying the signal by a factor . The Pearson correlation between predicted and experimental maps at the 16-kb resolution is 0.93.

We also compute Pearson correlations for the 16 individual chromosome maps at 1 kb resolution (Additional file 2: Table S1), giving an average of 0.88.

We compute aggregate plots around centromeres (on diagonal plots, see below for technical details) normalizing by the intra-chromosome P(s) (Fig. 1D). Our model correctly predicts the average “flame” shape around centromeres in experimental maps even if we observe a too strong depletion of contacts between centromeres and chromosome arms, compared to in vivo maps (blue stripes in Fig. 1D).

We then quantify the enrichment of centromere-centromere (trans) contacts arising from the imposed colocalization of centromeres in the confining shell around SPB by computing off-diagonal aggregate plots. To compare experiments and simulations, we normalize both by dividing by the respective P_intra(s = 10 kb) (Fig. 1C). This strategy allows to recover the relative enrichment with respect to the intra-chromosomal contacts. In particular, we choose a scale above the polymer persistence length where our simulations match well the in vivo signal (see Additional file 1: Fig. S2). In the simulations, we estimate that the average centromere-centromere contact frequency is ∼ 11% of the P_intra(s = 10 kb) while the experiment exhibits a more moderate enrichment ∼ 8%.

In both cases, these quantitative differences are likely due to the strong modeling assumptions made. For instance, chromosomes are strictly bound by the strong potential to localize between 250 and 325 nm from the SPB. Relaxing these conditions should decrease the centromere-centromere contact frequency. Similarly, the strong potential brings telomeres and centromeres strictly to different parts of the nucleus (NE and SPB respectively). This polarity is enforced in all trajectories, leading to a very strong depletion of contacts between centromeres and chromosome arms (blue stripes in Fig. 1A). The segregation between centromeres and telomeres might be more heterogeneous in vivo and other biological processes (such as condensin activity) may hinder this effect.

To compare simulated and experimental P(s) in Additional file 1: Fig. S2B,C and S13B, we multiply simulated curves by a factor α, which minimizes the Mean Squared Error (MSE) within a given genomic interval [s_min,s_max]:9

Practically we use s_min = 10 kb and s_max = 1 Mb when comparing P(s) averaged among different chromosomes (Additional file 1: Fig. S2B and S13B) while from s_min = 10 kb and s_max = 50 kb in the case of individual chromosomes (Additional file 1: Fig. S2C).

When comparing on-diagonal aggregate plots in different simulated conditions or between in vivo and in silico data, no additional normalization was used. In fact, for these particular analyses (such as in Fig. 1D or Figs. 2 and 3), we present Observed over Expected maps where each pixel (i,j) is normalized by the average contact frequency at the corresponding genomic distance s =|i − j| (P(s) signal).

We define RFi in silico by constructing a graph based on the 3D positions of forks, , at each simulation frame. For a given threshold 3D distance θ, we define an association matrix A_i,j, where A_i,j = 1 if |− |≤ θ and A_i,j = 0 otherwise. An RFi is defined as an independent component of A_i,j (i.e. a set of forks where each fork in the set is separated by a distance > θ from any forks outside the set). The total number of RFis corresponds to the number of graph components. To efficiently compute this quantity, we use the Python library networkx [114].

To compare our predictions with the experimental estimates of forks within RFi by Saner et al. [20] (Fig. 5D), we adapt our cluster detection to mirror the experimental protocol. In the original paper [20], the authors analyze 3D microscopy images of replicating yeast nuclei where PCNA (forks subunit) has been labeled with GFP proteins. Each distinct bright spot detected is then classified as an individual RFi and the number of forks inside is computed on the basis of an in silico estimation of 302 forks in the nucleus at that specific replication stage. In practice, the relative intensity of each focus (signal of the focus divided by total signal) is multiplied by the total number and cast strictly to an even number. This assumption is motivated by early microscopy evidence that sister-forks are associated [19]. Experimental estimate of forks pairs within each RFi (gray squares) was extracted from Fig. 3E of the original paper [20]. In our in silico predictions, we roughly follow a similar approach by dividing the number of forks of each RFi by the total number of forks in the system in that specific frame to obtain a relative signal. For each RFi, we then multiply this value by the average number of forks between t = 2 min and 30 s and t = 3 min 45 s (229 forks). Finally, the resulting number was made even. In Additional file 1: Fig. S19, we explore confounding factors that may impact the classification.

In order to simplify the analysis, all the observables described below are computed for only one chromatid for each chromosome. In particular, once a given monomer is replicated, each measurement is computed only for its copy belonging to the Sister Chromatid 1.

In Fig. 6, Additional file 1: Fig. S20C,S21 and S22, we compute the Mean Squared Displacement with the formula:10

In Additional file 1: Fig. S20B, we average over the monomer (mon) in the centromeric or telomeric regions, defined as being at a genomic distance less than 20 kb of a centromere or telomere respectively. Similarly, in Additional file 1: Fig. S22 the MSD is computed restricting the average to monomers belonging to individual chromosomes. In Additional file 1: Fig. S22B, we impair the replication of chromosome 8 by manually deleting all the potential origins before the onset of replication. In contrast, “instantaneous” replication was achieved by setting a homogeneous high firing rate for all origins and high replication speed (Additional file 1: Fig. S22C) for all chromosomes or exclusively for chromosome 8 (and hindering replication for all the others, Additional file 1: Fig. S22D).

The average Squared Displacement at 15 s in Fig. 6B, is computed with the formula:11

To rescale the time, we recorded each individual SD trajectory along with its replication status and aligned all trajectories at the moment of replication (t = 0 min) (see Additional file 1: Fig. S23A for an example of a single monomer in a single simulation).