Our understanding of how chromosomes structurally organize and dynamically interact has been revolutionized through the lens of long-chain polymer physics. Major protein contributors to chromosome structure and dynamics are condensin and cohesin that stochastically generate loops within and between chains, and entrap proximal strands of sister chromatids. In this paper, we explore the ability of transient, protein-mediated, gene-gene crosslinks to induce clusters of genes, thereby dynamic architecture, within the highly repeated ribosomal DNA that comprises the nucleolus of budding yeast. We implement three approaches: live cell microscopy; computational modeling of the full genome during G1 in budding yeast, exploring four decades of timescales for transient crosslinks between 5kbp domains (genes) in the nucleolus on Chromosome XII; and, temporal network models with automated community (cluster) detection algorithms applied to the full range of 4D modeling datasets. The data analysis tools detect and track gene clusters, their size, number, persistence time, and their plasticity (deformation). Of biological significance, our analysis reveals an optimal mean crosslink lifetime that promotes pairwise and cluster gene interactions through “flexible” clustering. In this state, large gene clusters self-assemble yet frequently interact (merge and separate), marked by gene exchanges between clusters, which in turn maximizes global gene interactions in the nucleolus. This regime stands between two limiting cases each with far less global gene interactions: with shorter crosslink lifetimes, “rigid” clustering emerges with clusters that interact infrequently; with longer crosslink lifetimes, there is a dissolution of clusters. These observations are compared with imaging experiments on a normal yeast strain and two condensin-modified mutant cell strains. We apply the same image analysis pipeline to the experimental and simulated datasets, providing support for the modeling predictions.
The 4D Nucleome Project [1] proposes the integration of diverse approaches: increasingly powerful chromosome conformation capture techniques including high-throughput chromosome conformation capture (Hi-C); statistical and topological analyses of these massive Hi-C datasets; 3-dimensional (3D) and 4D super-resolution imaging datasets; and computational modeling approaches, both constrained by and independent of Hi-C datasets. The project aims to gain mechanistic understanding of 3D structure and dynamics of the genome within the nucleus, and to learn how the active chromosome architecture facilitates nuclear functions. In this paper, we contribute to these aims by combining three approaches: (i) live cell microscopy for experiments studying the effect on gene clustering for normal and condensin-modified mutant cell strains; (ii) first-principles-based, computational modeling based on the statistical physics of chromosome polymers coupled with transient gene-gene crosslinks formed by condensin proteins; and (iii) analysis of the dynamic chromosome architecture with temporal network community detection algorithms applied to 4D modeling datasets across four decades of crosslinking timescales. Our present understanding of basic principles that govern high-order genome organization can be attributed to incorporation of the physical properties of long-chain polymers [2–7]. The fluctuations of long-chain polymers, numerically simulated with Rouse-like bead-spring chain models of chromosomes confined to the nucleus, capture the tendency of chromosomes to self-associate and occupy territories [8–11] In addition, these models make predictions with regard to the spatial and dynamic timescales of inter-chromosomal interactions, a dynamic analog of topologically associated domains. The convergence of robust physical models with high-throughput biological data reveals the fractal nature of chromosome organization, namely an apparently self-similar cascade of loops within loops, or structure within structure, as one examines chromosomes at higher and higher resolution [12–14]. De novo stochastic bead-spring polymer models of the dynamics and conformation of “live” chromosomes, plus the action on top of the genome by transient binding interactions of structural maintenance of chromosome (SMC) proteins, e. g. condensin, provide complementary information to chromosome conformation capture (3C) techniques, genome-wide high-throughput (Hi-C) techniques, and restraint-based modeling [12,15–20]. 3C and Hi-C experiments rely on population averages of gene-gene proximity on all chromosomes over many thousands of dead cells whose chromosomes have been permanently crosslinked by formaldehyde; the restraint-based modeling approach then explores 3D chromosome architecture that optimizes agreement with the experimental data on gene-gene frequency and proximity across the genome. Many powerful inferences have been drawn from both Hi-C and polymer modeling approaches, using analyses of empirical and synthetic datasets encoding maps related to the pairwise distances between genes. A common major limitation for existing polymer models and whole-genome contact maps in mammalian cells is in mapping two essential regions of the chromosome, namely the centromere and the nucleolus. The centromere, essential for chromosome segregation, and the nucleolus, the sub-nuclear domain of ribosomal DNA, are comprised of megabases of repeated DNA (centromere satellites and nucleolus rDNA). Furthermore, these regions are not captured in methods used for generating contact maps. We have used single live cell imaging of the nucleolus in budding yeast coupled with whole genome polymer modeling to explore the minimal requirements for sub-compartmentalization. Implementation of protein-mediated cross-linking within the nucleolus is sufficient to partition this region of the genome from the remaining chromosomes. Furthermore, stochastic polymer models reveal that the relative timescales of crosslinking kinetics and fluctuations of the chromosome chains have a profound influence on nucleolar morphology [21]. In single cells, the positional fluctuations of tagged DNA sequences on specific chromosomes [22–24] through the lac operator/lac repressor reporter system validate the bead-spring models. Chromosomes fluctuate as predicted for the conformational dynamics of idealized Rouse chains [25]. Polymer simulations over the entire genome have revealed the ability of relatively fast binding and unbinding, and thereby short-lived (fraction of a second) protein crosslinks to concentrate the rDNA chain sequence in a smaller volume and increase the simulated fluorescent signal intensity variance when the model datasets were convolved with a point spread function to create two-dimensional, maximum intensity projections [21]. Visualizations of the monomers in the simulations revealed that the fastest kinetics explored, or shortest-lived crosslinks (∼. 09s), generated several clusters of high polymer density, and overall compaction of the nucleolus. In contrast, much slower kinetics (decades longer-lived protein crosslinks (∼ 90s) ) tended to homogenize the fluorescent signal intensity as evidenced in the decrease in simulated fluorescent signal intensity variance. These model visualizations were consistent with experimental results on live budding yeast. There is a growing interest to analyze Hi-C datasets and model chromosome interactions using network models [26–28], which has opened the door to study chromosomal datasets using network-based algorithms including centrality analysis [29,30] and community detection [31,32]. In this context, a ‘gene cluster’ is a set of genes that are in close physical proximity, and it is represented in a network by a community of nodes (i. e. , a set of nodes between which there is a prevalence of edges). These detection algorithms perform an unbiased search for robust structures (communities or clusters) at the scale they exist in an automated manner, quantifying how chromosome conformational changes can precede changes to transcription factors and gene expression [33,34] and leading to new approaches for cellular reprogramming [29,35]. Here, we apply temporal community detection algorithms including multilayer modularity [36,37] to simulated 4D datasets over four decades of SMC-binding kinetic timescales. This approach integrates both temporal and spatial information so that each community now represents a set of genes that are not only nearby one another, but they remain in close proximity for some duration. This approach allows us to detect, track and label transient gene communities (clusters) in the nucleolus. Simultaneously, we record summary statistics on the sizes and numbers as well as persistence times of communities, and the frequencies of community interactions leading to gene exchanges. We likewise record standard bead-bead summary statistics. In doing so, we detect spatial and temporal organization at the scales they exist, beyond two-point (gene-gene) spatial proximity statistics. We identify the timescales over which spatial organization persists, linking the timescales to the cluster identification algorithm. Since clusters can deform through the flux of genes into and out of clusters, we further are able to identify crosslink timescales for which spatial clustering persists over extended timescales, and whether individual clusters are relatively permanent or experience frequent interactions and gene exchanges. Perhaps the most striking prediction of our modeling and data analysis is that specific gene organization tasks (amplified below) are optimized at a relatively short crosslink timescale, on the order of. 19 sec. With these network tools applied to physics-based 4D nucleome simulated datasets, we explore the mechanistic basis for the experimentally observed variance in nucleolar morphology. From a high-resolution sampling of the timescales for crosslinking of 5k base pair (bp) domains, 4D model simulations of the yeast genome reveal the nucleolus on Chromosome XII undergoes a stark transition in dynamics and structure, and does so within a narrow “mean on” crosslink timescale range of. 09 − 1. 6 sec. A highly stable clustering regime exists with relatively short-lived crosslinks (. 09 sec), with relatively few cluster interactions and gene exchanges, as reported previously in [21]. At slightly longer-lived (. 19 sec) timescales, a novel “flexible” behavior is revealed. Gene clusters continue to self-organize, yet clusters are more mobile, frequently interact, and exchange genes. Indeed, there is a peak timescale, marked by highly mobile gene clusters, at which both pairwise and community-scale gene interactions are maximized. As the binding affinity of crosslinker proteins increases only slightly longer (1. 6 sec), the community-scale structure has dissolved, with no identifiable nucleolar sub-substructure. See Fig 1. From a methods perspective, our analysis of the 4D simulated datasets is based on network modeling and a temporal community-detection algorithm known as multilayer modularity [37]. From a biological perspective, this tunable dynamic self-organization reflects a powerful mechanism to coordinate gene regulation and the coalescence of non-contiguous genes into identifiable clusters (substructures). The transition shown in Fig 1 occurs within such a narrow crosslinker timescale regime (. 09 − 1. 6 sec), suggesting a relatively simple mechanism to control dynamic sub-organization of the genome; indeed we performed and report experiments below to support this prediction. Finally, we emphasize the counter-intuitive nature of this mechanism: clustering is most often associated with segregation, however we observe that the dynamic element of flexible clusters facilitates an overall increase in global gene interactions in the nucleolus.
We first focus in the relatively short crosslink timescale regime, extending the simulations of [21] at discrete values μ = 0. 09,0. 9,90. These will establish a basic understanding of how the kinetic timescale μ for crosslinking sensitively affects the organization of the nucleolus and the dynamics of the architecture. From our refined simulations across the above four decades, the essence of the story can be told with results for three selected values μ ∈ {0. 09,0. 19,1. 6}. In Fig 1 (A) –1 (C), we present visualizations, i. e. , “snapshots, ” of the beads’ 3D positions during the simulations. The nucleolus on Chromosome XII is highlighted in blue and all remaining chromosome arms are colored gray. In Fig 1 (D) –1 (F), we show only the nucleolar beads, which are colored according to the network community detection analyses that we describe in the following sections. We also show videos of the time evolution of the beads, along with a simulated microscope projection, for each timescale in S1, S2 and S3 Videos. Based on Fig 1 (A) –1 (F) and the videos, we identify three qualitative regimes for nucleolus clustering: We will continue to use this terminology when referring to these three clustering regimes. We note that [21] discovered the two extreme regimes: robust clusters for μ = 0. 09, and the lack of clusters for μ = 90. As they did not finely sample the decades of timescales in between, they did not discover that the transition from robust to no clustering is in fact non-monotone with respect to gene-gene interactions, nor that the transition is essentially complete already at μ = 1. 6, and that the most biologically interesting and relevant regime occurs at μ = 0. 19. Furthermore, without automated structure detection algorithms, they would not have been able to detect and dynamically track clusters of genes and their interactions that explain the peak in gene-gene interactions at μ = 0. 19. This transition behavior and the optimal properties that arise will be the focus of several sections to follow. In Fig 1 (G) –1 (I), we show heatmaps of the bead-bead distances associated with the bead positions of the snapshots in (D) – (F), identical to those in (A) – (C); construction of heatmaps is described in the Methods Section: Pairwise-distance maps for high-throughput chromosome conformation capture (Hi-C). Heatmaps are widely used in Hi-C to depict population averages of pairwise gene-gene proximity data [16,39–43] and in simulated data from polymer bead-spring models, both from 3D snapshots and time averages [9,11,21]. Comparing the second and third columns of Fig 1, we note the difficulty (false negatives and false positives) in detecting the presence of structure and sub-organization in column 2 from visual examination of heatmaps in column 3. As shown in [21], the time average of 4D simulated datasets, even in the strong clustering regime, wipes out the sub-structure of snapshots when averaging over the entire G1 phase. An alternative approach has been to use polymer modeling to generate chromosome conformations, and to select those conformations that best match Hi-C data, so-called restraint-based polymer modeling [1]. Simultaneously, there have been efforts to develop methodologies to identify gene clusters in a rigorous and automated way from Hi-C data [26–28]. Our conclusion is that there is a need for a more reliable and objective method to study the clustering of chromosome domains in the nucleolus, especially spatio-temporal methods that take into account how bead positions and sub-organization change with time, weighing both spatial proximity and temporal coherence in the detection method. In the following sections, we present a scalable and automated technique to identify and track the dynamics of clusters. First, however, we will present new experiments that provide empirical evidence for clustering in the nucleolus. We conducted experiments to qualitatively compare image-based cluster analysis between our model and empirical measurements obtained from live cell microscopy and demonstrate the effect that SMC protein mutation can have on clustering in the nucleolus, extending the results previously reported in [21]. Here, we study three yeast strains: wild-type (WT), fob1 and hmo1. Importantly, fob1Δ and hmo1Δ are mutations that lack key proteins reported to crosslink or loop segments of rDNA within the nucleolus. Fob1Δ is required for maintenance of the rDNA copy number and regulates the association of condensin with rDNA repeats [44,45]. The replication fork barrier within the rDNA is a binding site for Fob1Δ that, together with several other components (Tof1, Csm1 and Lrs4), are responsible for the concentration of condensin within the nucleolus [44]. Hmo1 is an abundant high mobility group protein that localizes to the nucleolus and has been proposed to share functions with UBF1, which is involved in rDNA transcriptional regulation within the nucleolus [46,47]. Fob1 and Hmo1 are non-essential genes and were deleted from the genome to allow us to study their effect on nucleolus morphology due to functional modifications of crosslinking. See Methods Section: Yeast strains for experiment for further details. In Fig 2, we present images and analyses of nucleoli of these strains using fluorescent, live-cell microscopy. To visualize nucleoli, we fused Cdc14 protein phosphatase to green fluorescent protein (GFP) [21]. Nucleolar protein fusions occupy a distinct region of the nucleus that is adjacent to the nuclear envelope and (typically) opposed to the spindle pole body. We describe the image acquisition and processing steps in Methods Section: Image acquisition and baseline processing, and we highlight a few details here. Following image acquisition, we construct maximum intensity projections (MIP) centered on the nucleolus. See top row of Fig 2 (A). Due to potential variation in CDC14-GFP protein copy number and nucleolar/rDNA size from cell to cell, we normalized the nucleolar CDC14-GFP signal after excluding all intensity values below an intensity threshold. To this end, we first selected a threshold using Otsu’s method [48], which we implemented using the MATLAB function multithresh. One can interpret the threshold as a binary mask, as shown in the second row of Fig 2 (A). After applying the mask, we normalized the nucleolar signal by subtracting all intensities by the minimum value and then dividing them by the new maximum intensity that is obtained after subtraction. The third row of Fig 2 (A) depicts normalized images. Mutations of Hmo1 and Fob1 were found to alter the area and signal intensity of nucleoli labeled with CDC14-GFP across a range of intensity thresholds, which we surmise is due to alterations in the architecture of, i. e. , clustering within, the nucleolus. In Fig 2 (B) and 2 (C), we provide results for an analysis of nucleolar morphology: (B) the area of nucleolar signal; and (C) the standard deviation of the normalized signal. This analysis was implemented using the numerical algorithms presented in [21], which we further describe in Methods Section: Image analysis. As shown in Fig 2 (B), null mutations of hmo1Δ significantly altered the area of the nucleolar signal, whereas null mutations of fob1Δ did not. This was assessed by a Student’s two-tailed T-test, which yielded p = 3 × 10−8 for the former and p = 0. 07 for the latter. As shown in Fig 2 (C), the standard deviations of the normalized images were significantly lowered for the fob1Δ null mutation, but this did not occur for the hmo1Δ null mutation p = 0. 01 versus p = 0. 2). The non-significant changes are labeled ‘NS’ in the figure. The error bars indicate standard errors across n cells, where n = 84,70 and 77 for the WT, fob1Δ and hmo1Δ strains, respectively. We note that [21] also studied the area and variance of the nucleolus using experimental and simulated images. They found, for example, that the distribution of areas occupied by the nucleolus displays a lognormal distribution for WT cells in G1. Also, recall that we implemented thresholding based on Otsu’s method; in contrast, [21] explored a range of threshold values and found qualitatively similar results to be consistent across a range of threshold values. They did not, however, explore the area and variance for simulated images for a wide range of μ, which is the focus of our next experiment. To explore whether varying the kinetic timescale μ for our simulations yields similar changes as those arising under the fob1Δ and hmo1Δ mutations, we applied the microscope simulator of [21] to our 4D simulated data and analyzed the images using the same image analyses as described in Fig 2. First, we converted our 4D simulated data into a timelapse sequence with 22 time points, i. e. , snapshots. Each nucleolus bead was convolved with a point spread function and a maximum intensity projection was created for each timepoint. We depict 11 such images in Fig 3 (A). In panel (B), we plot the area of the nucleolar signal (computed using Otsu’s threshold) versus μ. Note that the nucleolus area increases as μ increases. In panel (C), we plot the standard deviation of nucleolar signal versus μ, which has the opposite trend. In Fig 3 (D), we plot the standard deviation of images obtained after a normalization step that is identical to that implemented for the experimental images (see discussion for Fig 2. Interestingly, the dependence on μ of the signal’s standard deviation drastically changes depending on whether or not it is normalized. Given that normalization is required to control for cell-to-cell differences in CDC14-GFP and in nucleolar/rDNA size, we sought develop a metric to measure clustering in the CDC14-GFP signal that was independent of the absolute values of the intensities. Our final experiment studies cluster formation in the nucleolus and compares clustering observed in the experimental and simulated microscopy images. We developed a cluster detection algorithm written with MATLAB (see Methods Section: Image analysis) and applied it to both the experimental and simulated images. We have made the code available at [49]. In Fig 4 (A), we depict images of maximum intensity projections for WT, fob1Δ and hmo1Δ strains with (top row) and without (bottom row) visualizations of detected clusters, which are represented by green circles. In panel (B), we depict identical information as in panel (A) except we show simulated images for three values of μ. In Fig 4 (C) and 4 (D), we show the number of clusters for the experimental and simulated images, respectively. For the experimental images, we give results for WT, fob1Δ and hmo1Δ, whereas for the simulated images we present results for μ ∈ [. 09,90]. We observe that the number of clusters was significantly decreased in the hmo1Δ null mutation, but not the fob1Δ null mutation (p = 0. 04 for hmo1Δ versus p = 0. 3 for fob1Δ). We also observe that increasing μ yielded a general trend in which there were fewer clusters. Taken together, these data suggest that gene clustering can directly impact the size and shape of the nucleolus. This underscores the need for robust and objective tools for identifying gene clusters. A simple and previously used method for analyzing distances between beads is to create a histogram of all bead-bead pairwise distances. As explored in [21], this two-point statistic can provide evidence of clustering and can be used to query simple properties such as whether or not the clusters change over time. In this section, we repeat this analysis on our data and extend it by showing what effects averaging over time and averaging over populations have on the results. We show that averaging one cell over time prevents observing the flexible clustering through pairwise distances, and averaging over populations prevents observing any sort of clustering. We provide further details on the computation of these distances in Methods Section: Pairwise-distance maps for high-throughput chromosome conformation capture (Hi-C). In Fig 5 (A), we plot the distribution of all pairwise distances { d i j (t) } at a single time t for three kinetic timescales, given by the same values μ = 0. 09,0. 19,1. 6 as shown in Fig 1. For μ = 0. 09, the pairwise distance distribution is clearly a multimodal distribution [21]. The peak near d ≈ 50 represents a large number of very short pairwise distances between beads in the same cluster. For slightly larger d, the density drops to zeros, indicating a separation distance between clusters. Interestingly, we observe two more peaks near d ≈ 300 and d ≈ 600. The clarity of these peaks suggests that the clusters themselves are regularly spaced from one another, reminiscent of a lattice structure. This shows the three layers of the multiscale structure of the nucleolus for μ = 0. 09: its existence as a dense, secluded section of the nucleus, the self-organization of intra-nucleolar clusters, and the individual beads within each cluster. For μ = 0. 19, one can also observe in Fig 5 (A) three peaks in the empirical probability density for bead-bead distances, but these peaks are much less pronounced. This shows a gradual transition in the degree of clustering as we increase μ. There is also a smaller gap between peaks. Together, these observations recapitulate our observations in Fig 1 (E), wherein the clusters can be observed to be less compact. Finally, for μ = 1. 6, there is no multimodal structure in the bead-bead distance plot. This is consistent with our expectation that there is no clustering structure present for this range of μ. The rigid and flexible clustering cases differ not only in how strong the clustering is at any given time, but also in how stable the structure is in time. We investigate this by considering how averaging pairwise-distances either across across time (Fig 5 (B) ) or over multiple simulations (Fig 5 (C) ) influences pairwise-distance probability densities. In Fig 5 (B), we plot the empirical probability densities for pairwise distances averaged across our 20 minute simulations. Note for μ = 0. 19 that the density is no longer multimodal, implying that aggregating the data across a large time range inhibits the detection of flexible clusters, which by definition change with time. Note that the rigid clusters, which are very stable across time, remain discernable as the pairwise probability density remains multimodal. Unsurprisingly, the slow crosslinking appears qualitatively very similar in the long time average, as there was no apparent structure in the first place. In Fig 5 (C), we plot the empirical probability densities for pairwise distances at a single time but averaged across 10 simulations with different random initial conditions. Note for all μ that there is no longer any multimodal structure for these densities, highlighting that averaging across heterogeneous cell populations obscures the detection of clusters. Next, we study how the kinetic time scale μ (i. e. , and thus the presence of clusters) affects the properties of pairwise gene interactions. A pair of beads is said to be interacting if they are in very close proximity and the distance between them drops below d*. As discussed in Fig 5, we choose d* = 100nm unless otherwise noted. In the following experiment, we show that increasing μ not only inhibits the formation of clusters, but that there exists a particular range of μ that optimizes gene mixing, or the overall interaction frequency of all pairs of genes. These experiments illustrate how clustering—which inherently describes multi-way relationships—can be studied through pairwise distances—which inherently describe two-way relationships —, and how there remain important open problems related to the time series signal processing of 4D chromosome conformation datasets. We study the following summary statistics for gene mixing: In Fig 6 (A) –6 (D), we plot these summary statistics across a wide range of μ. We identify three regimes that optimize different attributes. μ ≈ 0. 1 yields a self-organized structure that maximizes the number and the duration of gene-gene interactions (see panels (B) and (D) ). Recall from Video 1 that μ = 0. 09 yields many large clusters that are stable (i. e. , do not change) over time. This is reflected in a high number of interactions with beads in the same cluster and low number of interactions with beads not in the same cluster. With 0. 15 ⪅ μ ⪅ 1, we see flexible clustering behavior from Video 2. Notably, we find here that this flexible clustering has interesting properties beyond simply being a weaker version of the strong clustering from the rigid clustering regime. Namely, Fig 6 (A) shows that these μ values maximize the fraction of pairs of beads that interact at least once over the simulation, and Fig 6 (C) shows that these values minimize the waiting times between subsequent interactions. Thus, we can say that flexible clustering promotes the number of both simultaneous and overall distinct pairwise gene interactions in the nucleolus. This behavior arises from a balance between the number of intra-cluster gene-gene interactions, which is still elevated due to the moderate clustering as shown in (B), and the ability for genes to frequently switch between clusters during cluster interactions, as indicated by the reduced waiting time in (C). SMC proteins with such crosslinking timescales will thereby promote collective interactions among all active genes. These circumstances could accelerate a homology search, for example, to facilitate DNA repair, if the sister chromosomes were suddenly activated by a family of SMC proteins whose binding affinity was near this “sweet spot”. Finally, μ ⪆ 1. 5 is associated with a non-clustering regime, as shown in Video 3. The lack of clustering is reflected by a low number of gene-gene interactions, and the freely diffusing nature of the beads is reflected by short interaction duration and high interaction fraction. Having found that flexible clustering maximizes interesting properties of gene interaction, we seek to develop tools to identify and label the spatiotemporal clusters. In the rigid clustering regime, the clusters are so well-defined that any reasonable algorithm will detect them, but this is not the case for the flexible clustering. To detect and track flexible clustering, we utilize both spatial and temporal information to identify and track clusters. While we have access to 4D bead position time series data, we begin by transforming this into a multilayer network problem as described in Methods Section: Gene-interaction networks from pairwise-distance data. This is motivated by the fact that the most similar data available in biology, the Hi-C dataset, does not measure true distances between genome regions, but rather a notion of similarity based on average proximity [12]. The result of this transformation is a time sequence of weighted, undirected networks whose edge weights represent how near two beads tend to be to each other at that point in time. We refer to this sequence of networks as a temporal network. Given a gene-interaction network, we identify communities using an approach based on multilayer modularity [37]. See [31,32] for examples where community detection was applied to network models derived from Hi-C data. We present the algorithm in detail in Methods Section: Spatiotemporal gene clusters revealed by community detection in temporal networks, and we briefly describe it here. The modularity measure was originally introduced [36] to detect communities in a single, non-temporal network; it is a scalar that quantifies—as compared to a null-model lacking communities—the extent to which a network’s nodes can be partitioned into disjoint sets (i. e. communities) so that there is a prevalence of edges between nodes in the same community and relatively few edges between nodes in different communities. By searching over different possible ways to partition nodes into communities, one seeks to find an optimal partition that maximizes the modularity score [38]. Because each community contains a prevalence of edges, and edges only exist between pairs of genes that are in close proximity, a modularity-optimizing partition equivalently assigns genes into disjoint clusters so that each gene is nearer to genes in its cluster than to genes in other clusters. Our analysis is primarily based on an extended version of modularity that allows one to detect time-varying communities in temporal networks and is called the multilayer modularity measure [37]. In contrast to a community in a time-independent network (which is defined by a set of nodes), to specify a time-varying community one must also identify for each node the time-steps for which it is in the community. Using a variational technique [38], we optimize the multilayer modularity measure to simultaneously assign every node to a community at every time step. Each time-varying community in the network corresponds to a time-varying gene cluster, which is a set of genes that remain in close proximity for some duration. A key feature of the multilayer-modularity approach for community detection is that the framework involves two parameters, γ and ω, which provide “tuning knobs” [37,50] to identify, respectively, the appropriate sizes and temporal coherence of communities/clusters. Parameter γ is a resolution parameter [38] and allows one to select whether modularity-optimizing partitions involve many small communities or just a few, very large communities. Similarly, ω is a coupling parameter and allows one to choose if the communities can change drastically from one time step to the next or if they are restricted to changing slowly over time. We explored a range of choices to select appropriate values. Finally, we highlight that this approach significantly contrasts traditional clustering algorithms such as k-means clustering [51], which specifies the number of clusters a priori (i. e. , rather than allow the appropriate resolution to be dictated by the data) and which does not naturally extend to time-varying data. In S4, S5 and S6 Videos, we present equivalent videos to S1, S2 and S3 Videos, respectively, except in the new videos we color the beads according to the community labels that we detect. These new videos provided qualitative evidence supporting the ability of our algorithm to find clusters that have appropriate spatial and temporal scales. Looking at the rigid clustering in S4 Video, we see the coloring strongly agrees with our visual perception in the clusters. A similar but less decisive conclusion can be made from observing S5 Video. These videos indicate good agreement between visual perception of clusters and the clusters that are detected by the multilayer modularity algorithm—when beads visually appear to be clumped together, they tend to also be the same color in the videos, which reaffirms the validity of our choice of clustering algorithm. However, especially when looking at S6 Video depicting the non-clustering regime with slow crosslinking, we identify a key and common issue with clustering algorithms—they typically identify the “best” clusters, even when no clusters actually exist. In this section, we provide further evidence that flexible clustering is the mechanism that is responsible for the optimality observed in Fig 6. To this end, we will revisit and modify our definitions for gene interactions and gene mixing, which were defined at the “bead level” (i. e. , for pairs of beads). We now define similar, but slightly different, concepts that are defined at the “community-level” in that they reflect only community-membership information and do not require the precise bead locations. We say that two beads are “communicating” if they are in the same community. That is, all beads in the same cluster are communicating with each other, and beads in different clusters are not communicating. With this modified definition in hand, we define summary statistics for gene mixing at the community level that are analogous to the 2-point summary statistics for pairwise gene interactions that we previously defined in Section: Flexible clustering regime maximizes bead mixing. Analogous to gene mixing at the bead level, we now define “cross communication” at the community level. In Fig 7, we present summary statistics for cross communication that are analogous to our results in Fig 6 for pairwise gene interactions. Note that the results in Fig 7 are qualitatively identical to those in Fig 6, supporting our hypothesis that gene-level mixing is determined by community-level cross communication. That is, the formation of gene clusters and exchanges of genes between them governs the timescale at which nucleolar domains come in close proximity of one another. Using temporal community detection, we are now finally able to quantitatively support our first observations made in Section: Transient crosslinking timescale influences nucleolus clustering that there are three distinct clustering regimes: rigid clustering; flexible clustering; and no clustering. We support these observations by studying the properties of the detected clusters. In Fig 8 (A) –8 (C), we plot the average lifetime of clusters as a function of their average size (averaged over time). Panels (A) – (C) indicate the three clustering regimes with μ ∈ {0. 09,0. 19,1. 6}, respectively. For the rigid clustering regime, in Fig 8 (A), we see that most of the clusters are large, with an average size of 10 or more beads, and also have a long lifetime of over 100 seconds. This is consistent with our prior observations (e. g. Video 1) that showed large clusters that appeared very stable in time. We also see that the large clusters survive for much longer than the small clusters. For the flexible clustering regime, in Fig 8 (B), we can see the same general trend that larger clusters tend to have a longer lifetime than smaller clusters, but there is a much wider spread of cluster sizes, with only a moderate number of large clusters. For the non-clustering regime, there appears to be little relationship between cluster size and stability beyond an average size of approximately 3 beads. The clusters also tend to be much smaller, with almost no clusters with an average size over 10 beads. In Fig 8 (D) –8 (E), we plot the probability that a bead remains in the same community upon the next timestep, again as a function of cluster size. Panels (D) – (E) indicate results for μ ∈ {0. 09,0. 19,1. 6}, respectively. In agreement with panels (A) – (C), one can observe that larger clusters are more stable. Note also that the communities exhibit more plasticity for μ = 0. 19 than for μ = 0. 09 since beads have a higher average probability for changing the community to which they belong.
The dynamic self-organization of the eukaryote genome is fundamental to the understanding of life at the cellular level. The last quarter century has witnessed remarkable technological advances that provide massive datasets of both the spatial conformation of chromosomal DNA from cell populations (3C and Hi-C generalizations) and the dynamic motion of domains in living cells (GFP tagging and tracking of specific DNA sequences), from the yeast to the human genome. Data mining of this massive data has likewise witnessed remarkable advances in understanding the hierarchical packaging mechanisms of DNA that act on top of the genome, e. g. , histones and structural maintenance of chromosome (SMC) proteins, the topology of individual chromosome fibers, their topologically associated domains, and the territories they occupy in the nucleus. The third wave of advances has come from 4D modeling of chromosomes based on stochastic models of entropic, confined polymers, and the coupling of SMC proteins that either bind and crosslink genes on chromosomes or generate loops on individual chromosomes. As these three approaches continue to mature and inform one another, at an ever-increasing pace, insights into the structure and dynamics of the genome continue to deepen. The motivation for this paper lies in the information that can be inferred from these massive datasets, from Hi-C, live cell imaging experiments, and polymer physics modeling. In previous studies, cf. [21] and references therein, we showed that heterogeneity in experimental images derive from substructures that are formed within the nucleolus. We decided to use the array of tools from network community detection analysis, modeling, and associated fast algorithms, to automate the search for dynamic architectures in the 4D datasets generated by polymer modeling. We note a similar network analysis approach has been applied to Hi-C datasets [26–28], whereas our datasets have the added feature of highly resolved temporal information. Our aim was to infer organization beyond 2-point, time-averaged or population-averaged, gene-gene proximity statistics and heat maps generated from the statistics, and to remove the bias of an individual’s visual determination of structure. To do so, we used the advances in network-based models, their temporal generalization, data analysis, and algorithms, and applied this arsenal of tools to 4D datasets across four decades of crosslinking timescales to: (i) robustly identify clusters, or communities, of genes (5k bp domains, or beads, in our model), i. e. , to directly detect gene sub-organization at the scale it exists rather than attempt to seek larger scale organization from gene-gene statistics and heat maps; (ii) determine the size distribution (number of genes) in such sub-structures; (iii) determine the persistence times of communities; and, (iv) determine the interaction frequency of communities and corresponding gene exchanges, which are the drivers of gene-gene interaction statistics. In this way, network algorithms automate gene community detection and persistence, with robustness built in by enforcing insensitivity to algorithm tuning parameters. We elected to build and implement these network tools on the 4D datasets generated in house, from simulations of our recent polymer modeling of interphase budding yeast [21]. In this model, a pool of SMC proteins transiently and indiscriminantly crosslink 5k bp domains within the nucleolus on Chromosome XII. The kinetics of the cross-linking anchors relative to the substrate is a major driver of sub-nuclear organization. If the crosslinkers bind and release more rapidly than the chains can relocate, the non-intuitive consequence is that the chains explore less space. When dense clusters of crosslinker/binding sites arise, they persist for extended time periods when the crosslinking kinetics is sufficiently fast. We previously showed in [21] that very short-lived (μ = 0. 09sec) binding kinetics provided closer agreement with experimental results (highest degree of compaction of the nucleolus into a crescent shape against the nuclear wall). It was also shown via visualization of the simulated 4D datasets that this timescale induces a decomposition of the nucleolus into a large number of clusters each consisting of many 5k bp domains, and these clusters were persistent over time. On the other hand, with long-lived crosslinks (μ = 90sec) the clusters disappeared. These results reveal that the timescales of the crosslinkers relative to entropic fluctuations of the chromosome polymer chains are a fundamental contributor to genome organization. For the present paper, the sample set of binding kinetics in [21] was expanded to 4D datasets of interphase, sampling over four decades of crosslinking duration timescales. We applied standard, distance-based, 2-point statistical metrics and visualization tools, and then analyzed the full range of 4D datasets with the fast, automated network models and tools. The network algorithms search for and detect time-varying communities (clusters, sub-structures) at the scales they exist, not at a prescribed scale of two or more genes; the spatial and temporal scales are identified with the criteria that they are robust to the algorithm parameters. We then use this information to label and color-code communities, using community-level description to understand the persistence and interactions (merging and division marked by gene exchanges) between communities. With the above community-scale information and statistics, we generalize standard gene-gene interaction statistics across the four decades of bond duration timescale. As a generalization of waiting times for 2 distant genes to come within a specific distance of one another, we calculate waiting times for genes in the same community to leave and then re-enter another common community, and calculate the fraction of all genes that were in the same community at least once during interphase, which we call the community cross-communication fraction. From these analyses, we discovered a novel dynamic self-organization regime, wherein the rigid, persistent communities at relatively short-lived crosslink timescales (μ =. 09sec) transition at slightly longer-lived crosslink timescales (μ = 0. 19sec) to more mobile (literally, the clusters diffuse faster) communities that interact far more frequently, each interaction corresponding to merger, subsequent division, and an exchange of genes. We refer to this regime as flexible community structure with enhanced cross-communication. Furthermore, we discovered non-monotonicity in the dynamic self-organization behavior: the community cross-communication fraction is maximized, coincident with a minimum waiting time between genes departing and returning to common communities, with a crosslink timescale of μ = 0. 19sec. Both properties fall off, albeit in different ways, for shorter and longer timescales. We emphasize that these network tools and fast algorithms are amenable to any 4D dataset from polymer models. While we restricted the analysis in this study to the nucleolus during G1 of budding yeast where SMC proteins are allowed to transiently crosslink 5k bp domains, the same analyses can be applied to data with tandem SMC crosslinking and loop generation, for any cell type and for any phase of the cell cycle. Moreover, given the growing interest in network-based analyses for Hi-C data, network modeling is well-positioned to provide a fruitful direction for data assimilation efforts aimed at connecting simulated and empirical 4D chromosome conformation data. An important challenge facing this pursuit is the development of improved data pre-processing and community-detection methodology for temporal and multimodal network datasets [52,53].
Chromatin dynamics within confined yeast nuclei has been widely modeled using Rouse polymer bead-spring chains; see earlier citations. We employ the identical model and code in [21], resolving the entire yeast genome into 2803 total beads, each of which represents approximately 5k base pairs (bp). As such, each bead is interpreted as a chromosome domain or a tension blob, as explained in the polymer physics literature. The beads are arranged on 32 chromosome arms having lengths that reflect their experimentally identified lengths. Each arm is tethered at both ends to the nuclear wall: all emanating from the centromere at one end, with the other end tethered to one of six telomeres. Along each arm, there are entropic, nonlinear springs (the wormlike chain model is used here) that connect neighboring beads. The beads also experience Brownian, entropic, repulsive and hydrodynamic drag forces, and are physically confined to the nucleus. We simulate approximately 20 minutes of G1 during interphase, and each simulation is initialized with 32 chromosome arms tethered at the centromere and one of six telomeres on the nuclear wall, otherwise randomly located within the nucleus (idealized as a spherical domain). We model the effect of SMC proteins by transiently crosslinking pairs of non-neighboring beads in the nucleolus, represented by a contiguous chain of 361 beads on chromosome XII. (We note in passing that [21] also studied, experimentally and computationally, when the nucleolus is split onto separate chromosome arms, showing this to have a negligible affect on the clustering and interaction results). A transient crosslink is modeled by a pair of stochastic events—a binding and unbinding of two non-adjacent nucleolar beads—and the timescale of these events is tuned by a single parameter μ (measured in seconds). We defer to [21] for the details, but elaborate here on the role of μ. To implement crosslinks, all nucleolar beads are assigned a state of “active” or “inactive, ” and crosslinks are allowed only between active beads. Each bead’s state fluctuates stochastically as follows: an active bead becomes inactive after a duration that is a random number drawn the normal distribution N (μ, (μ/5) 2); and an inactive bead becomes active after a duration drawn from N (μ/9, (μ/45) 2). If two nucleolar beads are both active and the distance between them is less than d~ = 90nm, then a crosslink is formed between them (i. e. , they bind). Each bead may be involved in at most one crosslink at a time, decided on the basis of pairwise proximity among all active beads at each timestep. If either bead becomes inactive, then the crosslink is broken (i. e. , they unbind). Hence, crosslinks are established and broken stochastically in the nucleolus, and the single parameter μ dictates the kinetic timescales for crosslinking. See [21] for additional model details. Each simulation of the Rouse-like polymer model yields time-series data { x i (t) } ∈ R 3 that defines the 3D location of each bead i ∈ {1, …, N} at each discrete time step t = 0,1, …. We establish a connection between our simulated data and the state-of-the-art in chromosome imaging—namely, high-throughput conformation capture (Hi-C) —by constructing and analyzing pairwise-distance maps. Hi-C “images” the conformation of chromosomes using a combination of proximity-based ligation and massively parallel sequencing, which yields a map that is correlated with the pairwise distances between gene segments. While the actual pairwise distances between gene segments cannot be directly measured, Hi-C implements spatially constrained ligation followed by a locus-specific polymerase chain reaction to obtain pairwise count maps that are correlated with spatial proximity: the count between two gene segments monotonically decreases as the physical 3-dimensional distance between them increases. To provide an analogue to Hi-C imaging, we construct pairwise-distance maps for our simulated data {xi (t) }. Let F: { x i } i = 1 i = N ↦ R N × N (1) define a map (used here in the mathematical sense) from a set of N points { x i } ∈ R 3 to a matrix such that each entry (i, j) in the matrix gives the distance between point i and point j. Whereas Hi-C imaging aims to study the positioning of chromosomes using noisy measurements that are inversely correlated with pairwise distances, for our simulations we have access to the complete information about the chromosome positioning. We therefore define and study several variations for pairwise distance maps, which will allow us to also study artifacts that can arise under different preprocessing techniques, such as averaging the time series data across time windows and/or averaging across multiple simulations with different initial conditions. We define the following pairwise distance maps: These pairwise distance maps represent the data that is sought after, but cannot be directly measured, by Hi-C imaging. Moreover, by defining several distance maps we are able to study “averaging” artifacts that can arise due to various limitations of Hi-C imaging. For example, Hi-C imaging obtains measurements that are typically averaged across a large heterogeneous distribution of cells that are subjected to nonidentical conditions and exist at nonidentical states in their cell cycles.
The spatiotemporal organization of the genome plays an important role in cellular processes involving DNA, but remains poorly understood, especially in the nucleolus, which does not facilitate conventional techniques. Polymer chain models have shown ability in recent years to make accurate predictions of the dynamics of the genome. We consider a polymer bead-chain model of the full yeast genome during the interphase portion of the cell cycle, featuring special dynamic crosslinking to model the effects of structural maintenance proteins in the nucleolus, and investigate how the kinetic timescale on which the crosslinks bind and unbind affects the resulting dynamics inside the nucleolus. It was previously known that when this timescale is sufficiently short, large, stable clusters appear, but when it is long, there is no resulting structure. We find that there additionally exists a range of timescales for which flexible clusters appear, in which beads frequently enter and leave clusters. Furthermore, we demonstrate that these flexible clusters maximize the cross-communication between beads in the nucleolus. Finally, we apply network temporal community detection algorithms to identify what beads are in what communities at what times, in a way that is more robust and objective than conventional visual-based methods.
chemical bonding
chromosome structure and function
applied mathematics
simulation and modeling
algorithms
chromosome mapping
mathematics
materials science
network analysis
genome analysis
molecular biology techniques
cell nucleus
macromolecules
cellular structures and organelles
nucleolus
physical chemistry
research and analysis methods
polymers
polymer chemistry
computer and information sciences
genomics
gene mapping
chromosome biology
chemistry
cross-linking
molecular biology
cell biology
genetics
biology and life sciences
physical sciences
materials
computational biology
chromosomes
Transient crosslinking kinetics optimize gene cluster interactions