{'fraction_non_alphanumeric': 0.029043558359269236, 'fraction_numerical': 0.019399589677743505, 'mean_word_length': 4.5978125, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heterogeneous Spatial Distribution of Transcriptional Activity in Budding Yeast Nuclei

Recent microscopic and simulation studies have shown that the genome structure fluctuates dynamically in the nuclei of budding yeast Saccharomyces cerevisiae. This genome-wide movement should lead to the fluctuations of individual genes in their territorial regions. This raises an intriguing question of whether the resulting distribution of genes is correlated to their transcriptional activity. An effective method for examining this correlation is to analyze how the spatial distribution of genes and their transcriptional activity are modified by mutation. In this study, we analyzed the modification observed in a budding yeast mutant in which genes necessary for anchoring telomeres to the nuclear envelope, yku70 and esc1, are silenced. Taddei et al. reported that 60 genes are clearly misregulated by this mutation, with 28 and 32 genes downregulated and upregulated, respectively. We calculated the probability density maps of the misregulated genes using a model of dynamical movement of the yeast genome in both wild-type (WT) and yku70 esc1 mutant and showed that the density of downregulated genes is larger near the nucleolus, whereas the density of upregulated genes is larger at the opposite side of the nucleus. By comparing these genes with those highly (top 200 of transcriptome) and lowly (bottom 200) expressed, we showed that the simulated distribution of 28 downregulated (12 out of 32 upregulated) genes has a distinctly larger overlap with the distribution of lowly (highly) expressed genes in the mutant than in the WT. The remaining 20 upregulated genes are localized near the nuclear envelope both in the WT and in the mutant. These results showed that the transcriptional level of genes is affected by their spatial distribution, thus highlighting the importance of the structural regulation in the yeast genome.


INTRODUCTION
The eukaryotic genome is compactly packed in the nucleus. For example, the haploid genome of budding yeast, with a total length of~4 mm, is packed in the nucleus with a diameter of 2 mm. Understanding how the genome is organized in such a narrow space has been a long-standing problem in cell biophysics (1)(2)(3). In particular, the quantitative understanding of the genome structure and dynamics is important because they should regulate the access of proteins to DNA sequences (4). Recently developed experimental techniques shed new light on this problem. Particularly, the Hi-C method, which combines next-generation sequencing and chromosome conformation capture techniques, is important for inferring the average three-dimensional (3D) genomic structure (5)(6)(7)(8). The dynamics of chromosome structural fluctuation has been monitored for budding yeast using the LacO/LacI or TetO/TetR system marked with green fluorescent protein (9)(10)(11). Here, we focus on haploid budding yeast Saccharomyces cerevisiae as a model organism. Using the Hi-C results (12) and the results from fluorescence monitoring (10,11,13) as reference data, we develop a computational model to simulate the 3D structure and dynamics of the yeast genome to improve upon the model we previously reported (14). With this improved version of the model, we show that the largely fluctuating genome determines the spatial distribution of genes in the nucleus, and we discuss how the spatial distribution of genes is correlated to their transcriptional regulation.
Experimental and computational evidence suggests relationships between the spatial distribution of genes and their transcriptional regulation. For example, telomeres and subtelomeric regions in yeast are distributed around the nuclear envelope to coexist with silent information regulators (SIR proteins; Sir2, Sir3, and Sir4), which form complexes with other factors at around the nuclear envelope (4,15). The effect of SIR proteins to repress genes located next to telomeres is referred to as the ''telomere position effect'' (16)(17)(18). Some part of the repeats of ribosomal DNA (rDNA) in yeast is repressed through localization of rDNA near the nuclear envelope, which is mediated by inner nuclear membrane proteins including Heh1 and others (19,20). RNA Pol III-transcribed genes, such as transfer RNA (tRNA) genes in yeast, are localized around the nucleolus (21), and RNA Pol II-transcribed genes nearby the tRNA genes are repressed, which is a phenomenon termed ''tRNA gene-mediated silencing'' (22). In mammalian cells, actively transcribed genes are localized in discrete foci as transcription factories, which are associated with the localization of RNA Pol II (23-25). For fission yeast (26,27) and budding yeast (27), a similar localization of highly expressed genes was suggested using computational models of 3D genome structures. Moreover, it was suggested that functionally related genes belonging to particular gene ontology groups tend to colocalize (26,27). Localization of the functionally related genes was also assessed through statistical analyses of the Hi-C data (28-30).
One approach to examine the hypothesis that the spatial distribution of genes is correlated with transcriptional regulation is to analyze mutants in which the genome structure is arranged in a different way from the wild-type (WT). In the WT budding yeast nucleus, telomeres are localized around the nuclear envelope. Taddei et al. (31) generated a mutant in which Yku70 and Esc1, which are proteins necessary for anchoring telomeres to the nuclear envelope, are knocked out. They showed that the pattern of gene expression in the microarray data was significantly changed by this mutation (31). In particular, 60 genes were misregulated by >1.5-fold change in the transcription level, wherein 28 genes were downregulated and 32 genes were upregulated (31). They considered that genes located near telomeres, which detect the telomere position effect in the WT, are upregulated because they do not interact with the SIR proteins when telomeres are not anchored to the nuclear envelope in the mutant. They showed that seven upregulated genes are indeed located within 20 kb from telomeres, and argued that this portion (7/32 ¼ 22%) is higher than the probability of finding randomly selected genes near telomeres. However, the reason for the upregulation of the remaining 25 genes is not clear. The authors also demonstrated a statistical tendency for the downregulation of genes whose promoters contained RNA polymerase A and C promoters, ribosomal RNA processing elements, or Abf1 motifs. They considered that SIR proteins are sequestered in the SIR complex in the WT, which is formed together with Yku70 and Esc1, and that SIR proteins are released to the nucleoplasm and operate on these target promoters to downregulate the corresponding genes when the SIR complex is lost in the mutant. However, not all the downregulated genes have the above RNA polymerase A and C promoters, ribosomal RNA processing elements, or Abf1 motifs. A possible explanation is that many genes are misregulated by the action of genes that directly interact with SIR proteins; however, there is as yet no evidence for such an indirect mechanism. Taddei et al. (31) also argued that there is no evidence for the bias in the one-dimensional (1D) positions of misregulated genes distributed along chromosomes. In this way, the reason for the observed mis-regulation has not yet been fully elucidated, and other mechanisms than the one suggested in Taddei et al. (31) may also occur in the mutant. In this article, we consider the possibility that the shifted spatial distribution of genes in the mutant may lead to the shift in the pattern of gene expression, and this hypothesis is tested using a computer simulation.
When initiating a computer simulation, some basic problems must be considered: One is about the approach using Hi-C data in modeling. In the model of Duan et al. (12), the Hi-C data were used as major constraints to simulate a genome structure together with some minor constraints representing localization of centromeres and confinement of rDNA in the nucleolus. In other models, the Hi-C data were not used for modeling, and structures were simulated under the constraints to represent interactions between chromosomes and nuclear structures such as the spindle pole body (SPB), nuclear envelope, and nucleolus (27,32,33). In our previous model, constraints generated from the Hi-C data were used as moderate forces to restrict structural fluctuation. The Hi-C based interactions to represent these constraining forces and interactions to represent nuclear structures were used in a combined way in our previous model (14). Because of the small size of the yeast nucleus, interactions between chromosomes and nuclear structures have significant effects on the constraining of the genome structure; therefore, these interactions should be carefully considered during the construction of a model. Even without using the Hi-C data, the simulated genome structures were in reasonable agreement with the observed experimental features when interaction potentials between chromosomes and nuclear structures were suitably introduced (27,32,33). However, interactions between chromosome loci mediated by proteins binding on chromatins are not represented by the interactions between chromosomes and nuclear structures. We will show later in this article that interactions between tRNA genes and other chromosome loci, which were estimated from the Hi-C data (12), play important roles in determining the genome structure in yeast; therefore, the careful consideration of interactions between chromatin loci appears indispensable. Therefore, in this study, we use the combined interactions including the Hi-C based interactions by extending our previous model (14). The other problem to consider is the approach to represent fluctuations in the genome structure. Microscopic observations using of the fluorescent markers have shown that the budding yeast genome fluctuates dynamically in the nucleus (9,11,34). An equation of motion of the genome to simulate these dynamics has not yet been established; however, we here assume that a method analogous to the one used to describe protein folding dynamics can be used in this problem of genome fluctuation, although the physical scale of the two problems differs largely. In protein dynamics, potentials that stabilize the observed average structure (i.e., the native protein structure) can describe the fluctuating structural dynamics of proteins even for the process of large structural change, including unfolding/folding. Similarly, we expect that in genome folding, potentials that stabilize the observed average structure are useful to describe the large-amplitude liquidlike fluctuation of genome structure around the average structure inferred from the Hi-C data. We use the Langevin dynamics method to describe fluctuations under these potentials and to calculate the spatial distribution of genes in the nucleus.
In this simulation of Langevin dynamics, 16 chromosomes in haploid budding yeast are described as 16 polymer chains, where each chain is composed of beads representing 3 kb segments of chromatin connected by springs. These coarse-grained chains of chromosomes are assumed to be subjected to the Langevin equation of motion under the potentials representing the Hi-C based interactions and interactions between chromosomes and nuclear structures. This model is an improved version of the previous model (14). Three aspects of this version represent major improvements to the previous version: (1) interactions between tRNA genes and 5S rRNA genes were introduced to explain the localization of tRNA genes around the nucleolus as was observed by monitoring using fluorescence markers (21) and shown by the Hi-C data analysis (12,30); (2) the previous semilunar shape of potential function to represent the constraint for confinement of rDNA in the nucleolus was changed to a more realistic crescent shape; and (3) we removed biases from the Hi-C data by filtering out the noise from the data. The details of the filtering method are described in the Methods section. The yku70 esc1 mutant in which telomeres are not anchored to the nuclear envelope is represented in the model by turning off the attractive interactions between telomeres and nuclear envelope while these interactions are turned on in the WT model. In this way, we describe the yku70 esc1 mutation as a perturbation to the WT.
With this model, we simulate the genome structure and fluctuation in both the WT and the yku70 esc1 mutant and calculate the spatial distribution of the misregulated genes. We also calculate spatial distributions of highly and lowly expressed genes by extracting 200 genes showing the highest and lowest transcription levels, respectively (top 200 and bottom 200 genes, respectively) from the ArrayExpress database (Accession number: E-TABM-630) (31,35). By comparing these distributions, we discuss how the spatial distribution of genes is correlated to their transcriptional activity.

MATERIALS AND METHODS
In a previous study, we developed a model to simulate dynamic genome movement in interphase haploid budding yeast (Tokuda-Terada-Sasai (TTS) model) (14). Here, we use the improved version of the model to better reproduce the observed spatial distributions of centromeres, telomeres, rDNA, and other genes.
We describe 16 chromosomes as polymer chains that are composed of beads connected by finitely extensible springs. One bead corresponds to a 3 kb segment of chromatin, so that each chromosome has 78 (chromosome I) to 806 (chromosome XII) beads and the entire genome comprises 4460 beads. We simulate movement of these 16 chains by numerically solving the following Langevin equation: Here, m is the mass of a bead; r m i is the position of the ith bead of the mth chain with m ¼ 1-16; U is the potential representing interactions between chromosomes and those between chromosomes and inner structures of the nucleus; z is the friction coefficient; and the vector w m i is Gaussian white noise with where a and b represent the x, y, or z component of the vector, and T is the effective temperature with the unit of k B ¼ 1. Although the movement of chromosomes in the nucleus is out of equilibrium and the rate of structural change depends on the ATP consumption (4,15), we expect that the rapidly fluctuating movement of local segments to be much faster than the global movement of chains. Therefore, we use white noise to represent such rapid fluctuation and its strength is represented here by the effective temperature T. We use T as a unit to define interaction strengths in U. The potential U in Eq. 1 comprises three terms: where the first term describes interactions that define the basic physical features of the chain, the second term represents interactions derived from the Hi-C data (12), and the third term represents interactions between chromosomes and nuclear structures.
Interactions to define the physical features of the chain The first term in Eq. 3 is composed of two terms as U chain ¼ U ex þ U spring , where U ex represents the exclusive interactions between beads, and U spring represents a finitely extensible spring potential that connects the neighboring beads in chains. Explicit functional forms of U ex and U spring are explained in the Supporting Material. In addition to the above terms, a kinkable bending potential was included in U chain of the previous TTS model. This bending potential constrains the angle fluctuation between neighboring segments, which results in the chain of a 20-30 kb persistent length. We used this bending potential in the previous model because the fluorescence-in-situ-hybridization observation showed that chromatin chains in yeast cells should be kept straight for this distance (36). This persistent length was also consistent with the assumption that chromatin forms a regular 30-nm fiber. However, we should note that in the yeast nucleus, chromosomes tend to show stretched forms because centromeres are anchored to the SPB and telomeres are anchored to the nuclear envelope (4,15). Therefore, chains in the in situ condition are far from free polymer chains in solution, and chromosomes in the nucleus should be kept rather straight for 20-30 kb even when the persistent length, which should be defined in free solution, is sufficiently small. Because the regular 30-nm fiber is unstable against perturbations (such as the heterogeneous length distribution of linker DNA between nucleosomes, partial depletion of histone H1, or histone acetylation (37-40)), it is as yet uncertain whether the regular 30-nm fiber is indeed formed in the yeast nucleus. With the recent high resolution Micro-C analyses, no evidence for the regular nucleosome arrangement was found in yeast chromatin (41,42). Existence or absence of regular 30-nm fibers is under debate because the accumulating evidence shows that irregularly folded 10-nm fibers, rather than the regular 30-nm fibers, dominate the nuclei of various organisms (43). Here, we do not use the bending potential to constrain angles between segments, which corresponds to the assumption that the persistence length is <3 kb. Even with this assumption, the simulation results are not significantly different from those obtained with the bending potential because chromosomes tend to be stretched through interactions with nuclear structures. To define U ex , we use 30 nm as a measure of thickness of chromosomes by assuming that an irregularly folded 10-nm fiber should typically have the effective thickness of 30 nm (44) particularly in a coarse-grained representation, although the chain is more flexible than the regular 30-nm fiber.

Interactions derived from the Hi-C data
The potential U HiC in Eq. 3 represents interactions derived from the Hi-C data obtained by Duan et al. (12). In the previous TTS model, the potential representing the Hi-C based interactions was termed the G o-like potential because a similar structure-based potential termed the G o-like potential was used in protein folding studies as the potential to stabilize the observed average structure. In this study, we emphasize the effectiveness of the use of such structure-based potential to describe structural fluctuations. However, unlike the protein folding problem, the simulated genome does not fold into a unique structure, but rather continues fluctuating largely around the average structure. In this study, to avoid confusion, we refrain from referring to our potential as G o-like potential, and instead refer to Hi-C based potential.
We first remove biases from the observed Hi-C data. An interaction frequency obtained by the Hi-C measurement represents frequency demonstrating that two chromatin fragments in the genome are in physical proximity in an ensemble of cells. However, it has been shown that Hi-C data contain some biases because various features such as the fragment length, the guanine-cytosine content of the fragment, and the mappability of fragment affect the interaction frequencies (45,46). These biases depend on the definition of fragments used, and the fragments are determined by which restriction enzyme is used to excise the fragments from the genome. Therefore, we filter out these biases by comparing two sets of the Hi-C data obtained from two experiments using different restriction enzymes, HindIII and EcoRI, and extracting the common signals from them. From the libraries of fragments defined by the Tables S5 and S6 for HindIII and Tables S7 and S8 for EcoRI in Duan et al. (12), we first map each fragment to a bead in the model by assigning a bead nearest to the midpoint of the fragment. From this mapping, the matrix of interaction frequency f mn ij between the ith bead in the mth chromosome and the jth bead in the nth chromosome is defined for both HindIII and EcoRI, and we select the matrix elements that satisfy the condition where f total ðHindIIIÞ ¼ P mn P ij f mn ij ðHindIIIÞ and f total ðEcoRIÞ ¼ P mn P ij f mn ij ðEcoRIÞ are total frequencies. Using the threshold value x ¼ 4 for both intra-and interchromosome frequencies, we use f mn ij that satisfies Eq. 4 for estimating the average distance between beads by maintaining the value of f mn ij as it is, and set the other matrix elements that do not satisfy Eq. 4 as f mn ij ¼ 0. As discussed by Yaffe and Tanay (45), the effects of the filtering can be assessed by comparing the 1D coverage. Here we use the normalized 1D coverage, g m intra i;j f mm ij for the intrachromosome comparison in the chromosome m, and g m inter If biases are removed through filtering, two values of 1D coverage, g m a i ðHindIIIÞ and g m a i ðEcoRIÞ, where a ¼ ''intra'' or ''inter'' should show a large correlation. Fig. S1 shows the 1D coverage g m intra i with m ¼ 1 (chromosome I) as a function of i. We find that data for two 1D coverages, g 1 intra i ðHindIIIÞ and g 1 intra i ðEcoRIÞ, show a small correlation before filtering, as indicated by a Spearman correlation coefficient of r ¼ 0.2; however, this correlation coefficient increases to r ¼ 1.0 after filtering. In Fig. S2, g m intra i ðHindIIIÞ and g m intra i ðEcoRIÞ are compared for m ¼ 1-16 (the normalized intrachromosome 1D coverage for chromosomes I-XVI) for the data after filtering. Also for the interchromosome interactions, comparing two g m inter i s, it shows that the correlation is largely increased from r ¼ 0.5 to r ¼ 1.0 as shown in Figs. S3 and S4.
From the filtered interaction matrix, f mn ij ðHindIIIÞ, we estimate the average spatial distance rðm; i; n; jÞ between beads ðm; iÞ and ðn; jÞ. The relation between the intrachromosome frequency f mm ij ðHindIIIÞ and the average distance rðm; i; m; jÞ was estimated as summarized in Supplementary Fig. 17 of Duan et al. (12). Here, we use the same relation in this study. In addition, for interchromosome interactions, as was proposed by Duan et al. (12), the same relationship is used to estimate the average interchromosome distances, rðm; i; n; jÞ. We expect that the estimation could be further improved if we used the interchromosome fluorescence-in-situ-hybridization data for budding yeast as a reference, as was used for fission yeast by Tanizawa et al. (26).
Here, the Hi-C based potential has the following terms: The first term in Eq. 5 represents the interchromosome interactions, whereas the second term represents the intrachromosome interactions. The third term represents the effective interactions between tRNA genes and 5S rRNA genes. The first and second terms have the same functional form as U HiC À r mn ij Á ¼ À xr 0 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2pðrðm; i; n; jÞÞ 2 q Â exp À À r mn ij À rðm; i; n; jÞ Á 2 2ðrðm; i; n; jÞÞ 2 ! : Because the number of interacting pairs becomes smaller after filtration, we use a larger x than used in the previous TTS model before filtration. Here, we use r 0 ¼ 10 nm as a unit of distance, and x=T ¼ 50. Notably, the overall features of the simulated genome structure, such as the distributions of telomere positions, are not altered greatly from those obtained before filtration because they are mostly determined by interactions between chromosomes and nuclear structures; however, the spatial distributions of individual genes are largely modified because they are affected by the interactions between chromatins. Further attention should be paid to the interactions of tRNA genes. Although the sequence positions of tRNA genes in budding yeast are scattered over the genome, some are spatially localized around the nucleolus and the region where centromeres are gathered (4,12,15,21). Such localization would be natural when we consider that Pol III-transcribed genes such as tRNA and 5S rRNA genes, where the latter is confined in the nucleolus, may share a concentration of RNA polymerase III and enter close proximity. Because the frequency of interactions of 5S rRNA genes is not explicitly included in the Hi-C data, we additionally introduce the following attractive interaction between the ith bead corresponding to the tRNA gene and the jth bead corresponding to the 5S rRNA gene as where we use x t =T ¼ 10. We consider that the potential Eq. 7 does not represent interaction in the atomic scale, but rather a mean-field-like effective interaction for a statistical tendency of gene gathering; therefore, we use a relatively long distance as a parameter to define the potential as r t ¼ 200 nm. By analyzing the interacting frequencies in the Hi-C data, Duan et al. showed that tRNA genes are clustered into three groups; one localized around the nucleolus, one localized around the centromeres, and one not showing a distinct localization, Details of the clustered results of the Hi-C data are found on the web site (47). For the beads ðm; iÞ in Eq. 7, tRNA genes that were identified as genes localized around the nucleolus by the Hi-C analysis (47) are assumed to be targets of the potential in this model; however, the tRNA genes in the vicinity of telomeres are not considered as targets because they should be more attracted to the nuclear envelope in the WT. For ð12; jÞ in Eq. 7, we use beads representing 5S rRNA genes, which are located in chromosome XII. We assume that 5S rDNA has 150 repeats of sequences as j ¼ 154 þ 3k with k ¼ 0-149.

Interactions between chromosomes and nuclear structures
In budding yeast, 16 chromosomes are packed in the nucleus, and their motion is constrained by interactions between chromosomes and nuclear structures such as the nuclear envelope, nucleolus, and SPB. We represent these interactions by U nucleus in Eq. 3, which consists of the following three terms; where the first term U envelope represents interaction between chromosomes and the nuclear envelope, and the second term U cen represents interactions between centromeres and the SPB. The last term U nucleolus represents interactions that confine rDNA in the nucleolus and exclude other parts of chromosomes from the nucleolus. In this model, the semilunar-shaped functional form for U nucleolus in the previous TTS model is replaced by the more realistic crescent-shape functional form. To describe interaction terms in U nucleus , we define the coordinate system as shown in Figs. S5 and S6. With this coordinate system, the center of the nucleus is r center ¼ (1000, 1000, 1000) in units of nanometers, the nucleus is a sphere of 1000 nm radius, and the SPB is placed on (1000, 1000, 10). U envelope in Eq. 8 is composed of three terms, where the first term U telÀenv represents interactions between telomeres or subtelomeric regions and the nuclear envelope, the second term U rDNAÀenv represents interactions between the rDNA region and the nuclear envelope, and the last term U chrÀenv represents interactions between the other parts of chromosomes and the nuclear envelope. Telomeres and subtelomeric regions attractively interact with the nuclear envelope in the WT because they form a complex at the nuclear envelope in association with several factors: the inner membrane-associated protein Esc1 (48), the Sad1-UNC-84 domain protein Mps3 (49), the end-binding complex Yku70/80 (50), and the Rap1 (51) and SIR proteins (52,53).
In particular, SIR proteins in this complex promote heterochromatin formation near telomeres and repress genes in that location (52,53). Here, we represent the attractive interactions arising from the formation of this complex as follows: As a unit of interaction strength, we use e ¼ T. The values c 0 and c 1 are determined to smoothly connect the parts of U telÀenv . For simulating the yku70 esc1 mutant, we set y ¼ 0 and c 1 ¼ 0 (turned off), whereas y ¼ 1 and c 1 s0 (turned on) for simulating the WT. Parameters are chosen to set the telomere behaviors simulated with y ¼ 1 to be consistent with the experimental observations of the WT (10,13) as u 0 ¼ 800 nm, The rDNA repeats in budding yeast can be tethered to the nuclear envelope by several proteins such as Heh1, Lrs4, and Sir2 (4,15,19). Here, we consider that beads i ¼ 151-606 in the chain m ¼ 12 represent the rDNA repeats. We represent the attractive interactions due to this tethering as is a sum over the randomly chosen 30% of beads in the rDNA repeat region, and we use the same function for U rDNAÀenv ðR m i Þ as U telÀenv ðR m i Þ by maintaining y ¼ 1 and c 1 s0 in both the WT and the mutant.
The potential U chrÀenv ¼ P m P i U chrÀenv ðR m i Þ represents interactions between other parts of chromosomes and the nuclear envelope, which confine chromosomes inside the nucleus, where u 0 ¼ 800 nm and u ¼ 1000 nm are used. Explicit functional forms of other terms, U cen and U nucleolus in Eq. 8 are explained in the Supporting Material.

Simulation
We set m ¼ T ¼ 1. With this unit, t has the dimension of length L and the friction constant z has the dimension of L À1 . The friction constant z was set to be z ¼ 10 À5 to allow efficient sampling in the allowed computation time.
Values of other parameters are explained in the Supporting Material and summarized in Table S1. Starting from the structure inferred by Duan et al. (12), the structure was relaxed under the potential U with the steepest descent method. Then, using the Langevin equation of motion of Eq. 1, the high temperature run with T 0 ¼ 2T was performed, and the last structure of this high-temperature run was used as the initial structure of our simulation. Equation 1 was discretized with the step length Dt ¼ 0.01. We calculated up to 1 Â 10 5 steps, and from the last 5 Â 10 4 steps we sampled data in every 200 steps to obtain 250 structures. By calculating 40 simulation runs with different random number realizations, 250 Â 40 ¼ 10 4 structures in total were sampled for each analysis.

Simulated genome organization in the WT
To assess the performance of this model, we compare the simulated and observed (13) probability density maps of centromeres, rDNA, and the left telomere of chromosome VII in the WT. As shown in Fig. 1 a, the 3D distribution is projected onto a two-dimensional (2D) plane by defining a grid of 0.1 Â 0.1 mm spacing. The probability density maps are then calculated using the method explained in the caption of Fig. 1. Hereafter, this method of projection is used throughout the Results to show the calculated probability density maps of genes on the 2D plane. The simulated probability density maps in Fig. 1 b show an improved agreement with those of the observed data when compared with the previous version of the model (14). In Fig. S7 and Table S2, we compared the simulated and observed (11) distributions of the telomere-telomere distance in the WT, which shows that the calculated results explain the observed features of distributions. Two telomeres of a short chromosome, such as chromosome III, tend to reside near to each other, whereas two telomeres of a long chromosome, such as chromosome XIV, tend to be separate from each other. Therizols et al. (54) reported the median distances from subtelomeric regions of chromosomes IV, VI, and X to the other subtelomeric regions, which are compared with the simulation data in Fig. S8. The Pearson correlation coefficient is 0.52, which is similar to the value obtained by the model of Wong et al. (33) created without using the Hi-C based potential, suggesting that the distribution of telomeres is almost determined by the interactions between chromosomes and the nuclear envelope.
The simulation results for the WT model show intriguing features. An example of the simulated movement of the genome structure is shown in Movie S1. This movie shows that the genome largely fluctuates in the nucleus. Dependences of the motion on the friction coefficient z and the other parameters that can affect fluctuation of chains are discussed in the Supporting Material and Table S2. Through this genomewide movement, individual genes fluctuate in their territorial region. Therefore, the probability density maps of genes should provide useful information. The simulated probability density maps of 5S rRNA genes and several groups of tRNA genes are shown in Fig. 2. Near the nucleolus, we find the localized population of these tRNA genes, which is consistent with the cluster analysis of the Hi-C data (12) and the observed fluorescence localization of tRNA FIGURE 1 2D representation of probability density maps of various parts of the yeast genome. (a) A locus on an illustrated chromosome is represented by the cylindrical coordinate (R, a). The nucleus is represented as a sphere of radius 1 mm, and the z axis is defined as the axis running through the SPB (spindle pole body) and the center of the nucleus. A grid of 0.1 Â 0.1 mm spacing is defined on the 2D plane and a cylindrical tube running around the z axis is defined. This tube has the radius R sin a and the 0.1 Â 0.1 mm 2 square intersection. The probability density of beads at the position (R cos a, R sin a) on the 2D plane is calculated by the number of beads found in this cylindrical tube divided by the volume of the tube, which is further normalized by dividing by the total number of corresponding beads within the nucleus, and then averaged over genes (21). Here, each group comprises 9-16 copies of genes, and the probability distribution of the summed density of these multiple genes are shown in each panel of Fig. 2. Among these genes, only a few were identified by the cluster analysis of the Hi-C data (12) to localize around the nucleolus. Therefore, the attractive interactions between tRNA and 5S rRNA genes were assumed in the model only for those small number of genes. The ratio of the number of genes attracted to 5S rRNA genes to the number of copies of genes in the genome is 2/10, 3/14, 2/16, 2/9, and 4/14 for Fig. 2, b-f, respectively. Therefore, the major portion of the density was not attracted to the region near the nucleolus but concentrated around the center of the nucleus or around the region centromeres are gathered. It would be interesting to see whether the localized density around the nucleolus increases when a larger number of tRNA genes would have the attractive interactions; however, we use the assumption derived from the clustering analysis of the Hi-C data in a straightforward way. Thus, the study of various attractive interactions is not addressed in this article.
The simulated localization of the tRNA genes around either the nucleolus or the center of the nucleus is weakened when the attractive interactions between tRNA genes and 5S rRNA genes are turned off by setting U tRNA ðr m12 ij Þ ¼ 0 in Eq. 7; however, the density around the center of nucleus remains (Fig. S9), which is in clear contrast to the results of the model in which the Hi-C based potential is not used (27); in the latter model, the tRNA genes strongly concentrate around the region where the centromeres are gathered. Thus, the interactions between chromosome loci in the model do not affect the distributions of telomeres but considerably affect the distribution of genes.
Probability density maps of the highly expressed genes (top 200 genes) and the lowly expressed genes (bottom 200 genes) were calculated in the WT by extracting those genes from the microarray data deposited in ArrayExpress (E-TABM-630) (31,35). Notably, five upregulated and one downregulated genes were lowly expressed in the WT, and two upregulated and four downregulated genes were lowly expressed in the yku70 esc1 mutant; these genes were removed from the list of the above bottom 200 genes, so that there is no overlap of the gene names between the misre-  Table S3 for both the WT and the mutant. Fig. 3 shows maps of differential probability density for the top 200 genes and bottom 200 genes, DP ¼ (probability density of the genes) À (probability density of randomly selected  Heterogeneous Distribution in Yeast Nuclei sites). Here, the density of randomly selected sites was generated by selecting 100 beads randomly from every 10 4 structures; therefore, the density of randomly selected sites represents the density of chromatins in nucleus. In Fig. S10, the probability density maps before this subtraction are shown together with the density map of random sites. As shown in Fig. 3, highly expressed genes show a higher density in the inner regions of the nucleus and the lowly expressed genes show a higher density at the nuclear periphery. This is similar to the simulation results of Gong et al. (27), and is consistent with the localization of repressive factors such as SIR proteins near the nuclear envelope. It should be noted that 67 of the top 200 genes are ribosomal protein genes, which may be regulated in a different way from remaining 133 genes; many of the latter genes are related to basic metabolism. As evident in Fig. 3, it is interesting that the ribosomal protein genes are localized around the center of nucleus, whereas the other highly expressed genes are more localized around the region near the nucleolus and the outside part of the region where the centromeres are gathered with some additional population around the center of nucleus. This result suggests the existence of a mechanism of gene regulation that depends on the spatial position of genes.

Spatial distributions of transcriptional activities in the yku70 esc1 mutant
The yku70 esc1 mutant was simulated by setting y ¼ 0 and c 1 ¼ 0 in Eq. 10 to turn off the attractive interactions between telomeres or subtelomeric regions and the nuclear envelope. The validity of this calculation can be confirmed by comparing the simulated and observed radial distributions of telomere positions. In Hediger et al. (10) and Taddei et al. (55), the peripheral-most zone with 0.82 mm < R < 1 mm was defined as zone I, where R is the radial distance from the center of the nucleus as defined in Fig. 1 a. The probability to find the left telomere of chromosome XIV in zone I in the mutant is 40% in cells observed in the G1 phase (55) and 19% in the simulated data, and probability to find the right telomere of chromosome VI in zone I in the mutant is 38% in the observed data (55) and 44% in the simulated data. Therefore, the model has the ability to reasonably reproduce the observed radial distributions of telomeres in the mutant. Further discussion on the telomere residence probability in zone I is given in the Supporting Material. Fig. 4 is a snapshot of the simulated genome structure in the yku70 esc1 mutant, and the misregulated genes are drawn as spheres. We get the impression that the downregulated 28 genes (blue spheres) and the 32 upregulated genes (red spheres) do not appear to be homogeneously mixed, rather each group of genes shows some clustered features with heterogeneous density distributions. These misregulated genes are listed in Table S4. An example of the simu-lated movement of the genome structure in the yku70 esc1 mutant is shown in Movie S2. In this motion, the misregulated genes fluctuate largely based on their characteristic distributions.
The tendency of the distributions found in the snapshot and the movie is confirmed by plotting the probability density maps obtained from the average of 10 4 structures. In Fig. 5, the differential probability density maps are shown for upregulated (Fig. 5 a) and downregulated (  Fig. 5 e, we show the differential probability density map of highly expressed genes by subtracting ribosomal protein genes from the top 200 genes. From Fig. 5, we get the impression that the distribution of downregulated genes is similar to the distribution of the lowly expressed genes in the mutant.
The similarity or difference between spatial distributions can be quantified by calculating the overlap of 3D distributions of genes in different groups. We calculate the following degree of overlap between the spatial distributions of the bottom 200 genes and the 28 downregulated genes as where N bottom ¼ 200, N down ¼ 28, and d ¼ 100 nm; r k i is the position of the ith bead of the bottom 200 genes in the kth structure sampled from simulations; and r k j is the position of the jth bead of the downregulated genes in the kth structure. The term I bottomÀdown k is calculated for each of 10 4 structures with k ¼ 1-10 4 . In Fig. 6 a, the distribution of the calculated 10 4 values of I bottomÀdown k is shown by plotting the fraction, f bottomÀdown ðIÞ, which is the probability to find a structure k in the range I À DI=2 < I bottomÀdown with a suitable bin size DI. In addition, similar quantities, f topÀdown ðIÞ and f randomÀdown ðIÞ, are plotted together. Fig. 6 a shows that the overlap between the lowly expressed genes and the downregulated genes in the mutant are distinctly larger than the overlap between the random sites and the downregulated genes (p value < 10 À4 ) and that the overlap between the highly expressed genes and the downregulated genes tends to be slightly smaller than the overlap between the random sites and the downregulated genes (p-value < 10 À4 ), where the p values were calculated with the Wilcoxon signed-rank test. A similar argument holds for the 32 upregulated genes. In Fig. 6 b, distributions, f topÀup ðIÞ, f topðex:rib:ÞÀup ðIÞ, f bottomÀup ðIÞ, and f randÀup ðIÞ, are plotted. Fig. 6 b shows that the upregulated genes in the mutant tend to have a larger overlap with the highly expressed genes than with the random sites (p value < 10 À4 ), and they tend to have a smaller overlap with the lowly expressed genes than with the random sites (p value < 10 À4 ).
Comparisons among these degrees of overlap suggest that there are spatial actively transcribing region and inactively transcribing region in the yeast nucleus. The genes should be downregulated when they show higher probability of remaining in the inactive region in the mutant than in the WT. Similarly, the genes could be upregulated when they show higher probability of remaining in the active region in the mutant. This hypothesis, however, should be carefully examined because the difference between f topðex:rib:ÞÀup ðIÞ and f randomÀup ðIÞ shown in Fig. 6 b is not large. The reason for this small difference becomes evident when the 32 upregulated genes are classified into the following three groups: group 1, which comprises seven genes located within 20 kb from any of the telomeres; group 2, which comprises 12 genes with a tendency to reside in the active region; and group 3, which comprises 13 genes with a tendency to FIGURE 5 Simulated differential probability density maps, DP, of misregulated, highly expressed, and lowly expressed genes in the yku70 esc1 mutant yeast nucleus. Differential probability density maps of (a) 32 upre-  Table S5. Distributions of degree of overlap between the group-i genes and the top 200 genes except for the ribosomal protein genes, f topðex:rib:ÞÀupðiÞ ðIÞ with i ¼ 1-3, are shown in Fig. 7. We find that group 2 genes largely overlap with the highly expressed genes, whereas group 1 and group 3 genes show less overlap with the highly expressed genes in the mutant. The distribution f topðex:rib:ÞÀup ðIÞ in Fig. 6 b is obtained by averaging three distributions, f topðex:rib:ÞÀupðiÞ ðIÞ with i ¼ 1-3, each of which has the large difference from f randomÀup ðIÞ in opposite directions, resulting in the small difference between f topðex:rib:ÞÀup ðIÞ and f randomÀup ðIÞ.
The difference among three groups of upregulated genes is visualized with 2D maps of differential probability densities (Fig. 8). In the WT, the group 1 genes, which are located near telomeres, are strongly localized near the nuclear envelope, particularly around the SPB. Group 3 genes are also localized at the region near the nuclear envelope, whereas the group 2 genes are in various locations in the nucleus. In the mutant, group 1 and group 3 genes largely remain near the nuclear envelope. Therefore, the activity of group 1 and group 3 genes may be repressed in the WT due to their interactions with SIR proteins or other factors localized on the nuclear envelope, and the resolution of SIR complexes or loss of interactions with other repressive factors in the mutant should upregulate group 1 and group 3 genes. Group 2 genes do not show extensive interactions with the nuclear envelope in the WT, and largely shift toward the active region in the inside of the nucleus in the mutant, leading to their upregulation.
Such effects of spatial distributions can be further analyzed by examining changes in the overlap induced in the mutant. In Fig. 9, distributions of the degree of overlap in the WT and in the mutant are compared. Fig. 9, a-c, shows that the overlap between the downregulated genes and the inactive region is distinctly larger in the mutant than in the WT (Fig. 9 b), whereas the overlap between the downregulated genes and the active region is smaller   in the mutant than in the WT (Fig. 9 a). This result indicates that the spatial distribution of downregulated genes shifts toward the inactive region in the mutant. Fig. 9, d-l, shows the shift of upregulated genes. Group 1 genes, which are located near telomeres, show a small shift toward the active region ( Fig. 9 d) but not the inactive region ( Fig. 9 e). The overlap between group 1 genes and the random sites is somewhat increased (Fig. 9 f) because of the shift of group 1 genes from the nuclear periphery in the WT toward the inner region in the mutant due to the loss of telomere anchoring by the mutation. The transcriptionally repressive SIR complexes are formed near telomeres in the WT; thus the release of group 1 genes from the SIR complexes in the mutant can be sufficient for the upregulation as suggested by Taddei et al. (31). Group 2 genes, on the other hand, show a distinct shift toward the active region (Fig. 9 g) and a shift away from the inactive region (Fig. 9 h). Therefore, we hypothesize that these genes are upregulated because of their spatial shift toward the active region. Group 3 genes show a different behavior; they shift away from both the active (Fig. 9 j) and inactive (Fig. 9 k) regions. As shown in Fig. 8, group 3 genes are kept localized around the region near the nuclear envelope in the mutant, which separates these genes from both the active and inactive regions. The release of SIR proteins or other factors from the envelope in the mutant may explain the upregulation of these genes, but further examination is necessary to identify the interactions responsible for the upregulation of group 3 genes.
We showed that the simulated 28 downregulated genes shifted toward the inactive region and 12 out of the 32 upregulated genes shifted toward the active region in the yku70 esc1 mutant. Seven genes located near the telomeres were strongly localized around the nuclear envelope, and 13 genes were also localized near the envelope. These 20 genes did not show a distinct shift toward either the active or the inactive regions in the mutant; however, their localized distributions suggest that the loss of interactions with suppressive factors on the envelope leads to upregulation. Combining these simulation results, we hypothesize that the spatial distribution of genes and its shift induced by the yku70 esc1 mutation explain the observed misregulation.

DISCUSSION
We developed a computational model of a dynamically fluctuating genome structure for simulating the genome of the WT and the yky70 esc1 mutant of budding yeast. With this model, we showed that the loss of the anchoring of telomeres to the nuclear envelope changes the spatial distribution pattern of genes. We calculated the distribution of the genes misregulated by the mutation and the distribution of lowly (the bottom 200 genes in the transcriptome) and highly (the top 200 genes) expressed genes in both the mutant and the WT models and found that the downregulated genes are localized around the nucleolus, and the overlap of distributions of the downregulated genes and the lowly expressed genes is larger in the mutant when compared with the WT. We also found that the overlap of distributions for 12 out of 32 upregulated genes and the highly expressed genes is larger in the mutant when compared with the WT. The remaining 20 upregulated genes are localized near the nuclear envelope; therefore they may be affected by the change of interactions with factors on the nuclear envelope. These simulation results suggest that the transcriptional activity of genes is regulated by their spatial position in the nucleus. This mechanism of regulation could arise from the localization of RNA polymerases, transcription factors, or repressive factors such as SIR proteins in the nucleus. Further computational and experimental tests are required to verify this hypothesis. Particularly, the genomewide test for the generic genes is needed. This type of analysis is discussed in the Supporting Material.
To develop efficient computational tests, it is important to further improve the simulation model. In particular, a more accurate treatment of interactions between chromatin loci is important. As shown by previous computational models (27,32,33), simulations that do not consider the specific interactions between chromosome loci reproduce many aspects of the overall genome structure such as the distribution of telomere positions. Indeed, those simulated features do not differ significantly in this model, which explicitly considers the specific chromosome interactions. However, the distribution pattern of individual genes is more sensitive to the interactions between chromosome loci, which can, for example, be found in the simulated distribution of the tRNA genes. Given the interactions between tRNA genes and 5S rRNA genes in this model, tRNA genes tend to localize around the nucleolus consistently with the observation (12,21), and without these interactions, tRNA genes distribute more around the center of the nucleus or the centromere region. As shown in Fig. S12, the distribution of other genes is also affected by these interactions; the distribution of the bottom 200 genes shifts, particularly in the mutant, from the region near the nucleolus to the opposite side of the nucleus if the interactions between tRNA genes and 5S rRNA genes are turned off.
The results of our simulation showed that the budding yeast genome fluctuates dynamically in the nucleus, which determines the regions in which individual genes are localized. An interesting finding is that the ribosomal protein genes are highly expressed; however, their position is localized in a region different from other highly expressed genes. As shown in Fig. 3, the lowly expressed genes show the smallest density around the region of ribosomal protein gene accumulation, suggesting that the repressive regulatory factors are excluded or RNA polymerases are accumulated in this region. This result is a significant feature of the position-dependent regulation in the yeast nucleus, and should be examined by monitoring the distribution and dynamics of these genes and factors in the nucleus using fluorescent marker techniques.
In this study, we used the yku70 esc1 mutant genome as a perturbed genome from the WT. A similar analysis can be used for other mutants or perturbed cells. In particular, an artificial manipulation to chemically change the genome architecture is important (56); such genome engineering methods should provide a useful test platform to analyze the effect of perturbations. We consider that computational simulation models will play an important role in the elucidation of genome organization and function as tools for integrating results from different approaches, including a biochemical approach like the Hi-C method, a biophysical approach like the fluorescence monitoring of chromatin, and the genetic and chemical methods used for the genome organization.

SUPPORTING MATERIAL
Here, m is the mass of a bead, and r µ i is the position of the ith bead of the µth chain with µ = 1-16. U is the potential representing interactions between chromosomes and those between chromosomes and inner structures of the nucleus, and ζ is the friction coefficient. The vector w µ i is a Gaussian white noise with where α and β represent the x, y, or z component of the vector. T is the effective temperature with the unit of k B = 1. We use T as a unit to define interaction strengths in U . The potential U in Eq.1 is composed of three terms; where the first, second, and third terms describe interactions that define the basic physical features of the chain, interactions derived from the Hi-C data (1), and interactions between chromosomes and nuclear structures, respectively. In the following subsections, we explain these terms one by one.

1
Interactions to define the physical features of the chain The first term in Eq. 3 is comprised of two terms as U chain = U ex + U spring , where U ex represents the exclusive interactions between beads, where r µν ij = |r µ i − r ν j |. The first and second terms of Eq. 4 represent interand intra-chromosomal interactions, and the second term does not include interactions between neighboring beads. Here, except for interactions of the rDNA sites, we use the soft-core repulsive potential for U ex instead of the Lennard-Jones (LJ) type potential used in the previous version of the TTS model, so as to allow chains to cross each other to represent the action of topoisomerases that relieve DNA knots in the cell. This soft-core repulsive potential was used by Zhang et al . (2) and its use in simulation increases the sampling efficiency; with where = T and a = 30 nm are the same as in the previous TTS model (3), and will hereafter be used as units of energy and length. E cut = 4 and r s is determined to satisfy U ex (r s ) = 0.5E cut . For U ex representing interactions of the rDNA sites, we use the LJ-type potential; however, the range of its repulsive force is modified to be larger than in the previous TTS model by assuming that rRNA and other nucleolus proteins are associated with rDNA to result in a larger effective volume of rDNA as was assumed by Wong et al . (4). We use with a r =150 nm for the interactions between rDNA sites, and with a r = (30 + 150)/2 nm for the interactions between rDNA sites and other sites. Here, we regard beads i = 151 ∼ 606 in the chain µ = 12 as the rDNA sites. For the spring interaction connecting neighboring beads, we use the same potential as was used in the previous TTS model: where r µ i = |r µ i − r µ i+1 |, and θ(r µ i ) = 1 when r µ i < R 0 < R 0 and θ(r µ i ) = 0 when r µ i ≥ R 0 . The first term U FENE represents the finitely extensible nonlinear elastic (FENE) potential (5-7); where we use the same parameters as was used by Rosa et al . (6): R 0 = 1.15R 0 , and we assume R 0 = a = 30 nm. The parameter k was calibrated to be k = 6.0 /a 2 to result in a mean-segment length of approximately 23 nm, which is the value expected when the packing density of chromatin is 130 bp/nm (8); the calibration was performed by monitoring the simulated results for chromosome IV. The potential in the second term of Eq. 9 was introduced to stabilize the numerical behavior around r µ i ∼ R 0 with U l (r µ i ) = (0.93k/a 3 )(r µ i ) 5 . The third term U r (r µ i ) in Eq. 9 represents the excluded volume interaction between neighboring beads; Interactions derived from the Hi-C data The potential U HiC in Eq. 3 represents interactions derived from the Hi-C data obtained by Duan et al. (1). The interaction matrix obtained by the Hi-C method was filtered to remove biases as explained in the Method section of the present paper. From the filtered interaction matrix, f µν ij (HindIII), the average spatial distance r(µ, i; ν, j) between beads (µ, i) and (ν, j) was estimated. The Hi-C based potential was then defined as in the following way; where the first term represents the inter-chromosome interactions, and the second term represents the intra-chromosome interactions. The third term represents the effective interactions between tRNA genes and 5S rRNA genes. The first and second terms have the same functional form as with ξ/T =50, and r 0 =10 nm. The attractive interaction between the ith bead corresponding to the tRNA gene and the jth bead corresponding to the 5S rRNA gene is with ξ t /T = 10 and r t =200 nm. For the beads (µ, i) in Eq. 14, tRNA genes that were identified as genes localized around the nucleolus by the Hi-C analysis (9) are assumed to be targets of the potential in the present model; however, the tRNA genes in the vicinity of telomeres are excluded from the target because they should be more attracted to the nuclear envelope in the wild-type (WT). For (12, j) in Eq. 14, we use beads representing 5S rRNA genes located in chromosome XII. We assume that 5S rDNA has 150 repeats of sequences as j = 154 + 3k with k = 0, . . . , 149.

Interactions between chromosomes and nuclear structures
Interactions between chromosomes and nuclear structures are represented by U nucleus in Eq. 3, which has the following three terms; where the first term represents interactions between chromosomes and the nuclear envelope and the second term represents interactions between centromeres and the SPB. The last term represents interactions that confine rDNA in the nucleolus and exclude other parts of chromosomes from the nucleolus. The coordinate system is defined as shown in Supporting Figure  5, with the center of nucleus on r center = (1000, 1000, 1000) in units of nm, the nucleus is a sphere of 1,000 nm radius, and the SPB is placed on (1000, 1000, 10). U envelope in Eq. 15 is comprised of three terms, where the first term U tel−env represents interactions between telomeres or subtelomeric regions and the nuclear envelope, the second term U rDNA−env represents interactions between the rDNA region and the nuclear envelope, and the last term U chr−env represents interactions between other parts of chromosomes and the nuclear envelope. Telomeres and subtelomeric regions in the model attractively interact with the nuclear envelope in the WT with the potential U tel−env = µ i∈telom, subtelom U tel−env (R µ i ) as where R µ i = |r µ i −r center |. The values of c 0 and c 1 are determined to smoothly connect the parts of U tel−env . For simulating the yku70 esc1 mutant, we set y = 0 and c 1 = 0 (turned off), whereas y = 1 and c 1 = 0 (turned on) for simulating the WT. Parameter values are set based on the results of the telomere behaviors when simulated with y = 1 consistent with the experimental observations of the WT (10, 11) as u 0 = 800 nm, u 1 = 400 nm, u 2 = 800 nm, u 3 = 950 nm, u = 1, 000 nm, and c = 3.
The attractive interactions due to the tethering of rDNA to the nuclear envelope are described as U rDNA−env = i∈rDNA, µ=12 U rDNA−env (R µ i ), where i∈rDNA, µ=12 is the sum over randomly chosen 30% of beads that represent the rDNA region, and we use the same function for U rDNA−env (R µ i ) as U tel−env (R µ i ) by keeping y = 1 and c 1 = 0 for both the WT and the mutant. The potential U chr−env = µ i U chr−env (R µ i ) represents interactions between other parts of chromosomes and the nuclear envelope, which confine chromosomes inside the nucleus: where u 0 = 800 nm and u = 1, 000 nm are used.
Interactions between centromeres and the SPB are represented by spring potentials, where h/T =0.3 and s=100 nm. l µ is distance from the SPB to the centromere of the µth chromosome, and l µ 0 is the same distance in the structure inferred by Duan et al . (1).
In Eq. 15 the potential U nucleolus was introduced to represent the nucleolus of a crescent shape and occupying roughly one-third of the nuclear volume: Here, U rDNA nucl (R 12 i ) is the potential to confine rDNA in the nucleolus and U nucl (R µ i ) is the potential to exclude other parts of chromosomes from the nucleolus. Supporting Figure 6 shows how the nucleolus is defined in our simulation model, namely a sphere with the radius R n = 1, 670 nm is used to fit the confinement region to the observed distribution of the rDNA (11) by placing its center at r center = (1000, 1000, −340) in units of nm. Potentials are then defined as where R µ i = |r µ i − r center |.

Parameters
The parameters used in the model are summarized in Supporting Table 1.
Among the parameters used to characterize the physical features of chromatin chains, values of E cut and R 0 were chosen from other models (2,6). The value of k was chosen to reproduce the observed packing density of chromatins (8), and the value of a r was chosen to reproduce the microscopic observation of distribution of rDNA in the nucleolus (11). The value a = 30 nm was chosen by assuming that the irregularly folded 10 nm chromatin chain has an effective thickness of approximately 30 nm. Distributions of telomere and gene positions are largely determined by the balance among different types of forces in the model; attractive forces between telomeres and the nuclear envelope, forces acting between genes, which were derived from the Hi-C data, and forces between genes transcribed by RNA Polymerase III. Here, the forces between telomeres and nuclear envelope were characterized by u 0 , u 1 , u 2 , u 3 , c, and (see Eqs. 17 and 18 in the Supporting Text), the Hi-C based forces were characterized by the parameter ξ (Eq. 13 in the Supporting Text), and the forces between genes transcribed by RNA Polymerase III were characterized by ξ t and r t (Eq. 14 in the Supporting Text). The values of these parameters were chosen to reproduce the observed microscopic data of distribution of telomeres (11,12) as discussed by using Figure 1 of the main text and Supporting Figure 7, and to explain the localization of tRNA genes near the nucleolus (13) as shown in Figure 2 of the main text. The microscopic data of the radial distribution of telomeres in the nucleus of the yku70 esc1 mutant (14) were also used as constraints to choose values of ξ, ξ t and r t as explained in the Results section of the main text. The values of other parameters representing the nuclear constraints, R n and h, were chosen to reproduce the microscopic observation of the centromeres and nucleolus (11).
Notably, various features of the simulated results depend on the value of the strength of the Hi-C based interactions, ξ. Agreement between the calculated and observed distributions of the telomere-telomere distances shown in Supporting Figure 7 was improved when the larger value of ξ was used. However, we here used a moderate value of ξ/T = 50 to yield results consistent with a variety of different observations. Because the Hi-C based interactions act as attractive interactions between chromatin loci, an increase in ξ reduces the effective volume of the simulated mutant genome, which lowers the probability to find telomeres in zone I in the mutant nucleus, leading the probability to be ∼ 0 when ξ/T = 500. We used the moderate value of ξ/T = 50 to explain the experimental data for the zone I probability of telomeres in the mutant as discussed in the Results section, and also to explain the microscopic observations discussed with Figure 1, Figure 2 in the main text, Supporting Table 2, and Supporting Figure 7 in a balanced way.
The movement of chromatin chains depends on both the friction coefficient ζ and other parameters that can affect the rigidity of the simulated chain. When ζ is too small, the movement is exaggerated and the distributions of telomeres are modified to diminish the agreement with the microscopic observation. On the other hand, when ζ is too large, the movement is slowed down and the time length to reach a steady state is elongated. In the present paper, the value of ζ was chosen for the simulated trajectories to reach a steady state during the available computation time, and to reasonably reproduce the observed distribution of telomeres in the WT nucleus. To examine the effects of the other parameters on the chromosome movement, we calculated the dynamical component of the standard deviation, σ 0 , for the simulated telomere-telomere distance. We numerically followed the Langevin dynamics of 40 trajectories, which started from different initial structures, and 250 structures were sampled from each of these trajectories. σ 0 is the average of 40 standard deviations, each of which was calculated from 250 structures that appeared along a trajectory. Therefore, σ 0 is a measure of the fluctuation of the chromosome structure along individual dynamical trajectories.
In Supporting Table 2, the σ 0 of several telomere-telomere distances are shown for models having different parameters. The standard deviation, σ, and the average, d, of distances in 250 × 40 = 10 4 structures are compared with σ 0 . Different models are compared using a nonzero value of k φ , and a larger ξ and smaller c than the WT model used in the main text. Here, k φ is the elastic constant of a kinkable potential for the angle between successive segments in the chain, which was defined in Eq.8 of Ref. (3). k φ was set to be k φ /T = 2 in Ref. (3) to make the chain as stiff as the regular 30 nm fiber, and it was set to be 0 in the present WT model. The larger ξ implies the stronger Hi-C base interactions, and the smaller c represents the weaker telomere-nuclear envelope interactions than in the WT model. From Supporting Table  2, we find that the simulated σ and d are smaller than the observed ones in all models, but the simulated values show a consistent tendency that accounts for the difference among telomere-telomere distances. σ 0 is about 12-25% of d, showing the liquid-like large dynamical fluctuation in all models. It is interesting that σ 0 is large when ξ is large for the telomere-telomere distances having large d, showing that the large Hi-C based interactions enhance the dynamical fluctuation of the long distance structures; this effect should be due to the shallow Gaussian functional form of potentials we used for the long-range Hi-C based interactions in Eq. 3 in the Supporting Text.

Shift of spatial distributions of generic genes
Genes that show small changes in the expression level by the mutation also exhibit changes in the overlap with the active or inactive regions. The spatial shift of generic genes was examined by calculating the average overlap between the distribution of the bottom 200 genes and the distribution of the ith gene,Ī bottom−i , and the average overlap between distributions of the top 200 genes excluding the ribosomal protein genes and the ith gene, I top(ex.rib.)−i . Among the 5,643 genes available from ArrayExpress datafile (E-TABM-630), we found an increase inĪ bottom−i andĪ top(ex.rib.)−i through the mutation in 2,706 and 2,462 genes, respectively. In other words, in the yku70 esc1 mutant, 2,706 genes moved toward the inactive region and 2462 genes moved toward the active region. In addition, we define the relative expression level, E 0 (i), of the ith gene by E 0 (i) = E(i)/Ē, where E(i) is the expression level of the ith gene described in the datafile andĒ is the average of E(i) over all genes. Then, the change of the relative expression level by the mutation is δE 0 (i) = E 0 (i, mut)/E 0 (i, WT), where E 0 (i, mut) and E 0 (i, WT) are E 0 (i) for the mutant and the WT, respectively. Among 2,706 genes that moved toward the inactive region, 655 genes (24.2%) showed an increased δE 0 and 2,051 genes (75.8%) showed a decrease in δE 0 . On the other hand, among 2,462 genes that moved toward the active region, 630 genes (25.6%) showed increase in δE 0 and 1,832 genes (74.4%) showed a decrease in δE 0 .

Supporting Movie 1
An example motion of the genome in the WT budding yeast is shown for a sequence of 250 images taken every 200∆t of the simulation. 60 misregulated genes in the mutant are shown with spheres; 32 up-regulated genes (red) and 28 down-regulated genes (blue). The movie was constructed using Visual Molecular Dynamics (VMD) software.

Supporting Movie 2
An example motion of the genome in the yku70 esc1 mutant of budding yeast is shown for a sequence of 250 images taken every 200∆t of the simulation. 60 misregulated genes in the mutant are shown with spheres; 32 up-regulated genes (red) and 28 down-regulated genes (blue). The movie was constructed using Visual Molecular Dynamics (VMD) software. The differential probability density maps shown in Figure 3 of the main text were obtained by subtracting the probability distribution of randomly selected sites as a background from other probability density maps. The probability density is shown in gray scale; the lighter shade represents the higher density.  Figure 5 of the main text were obtained by subtracting the probability distribution of randomly selected sites as a background from other probability density maps. The probability density is shown in gray scale; the lighter shade represents the higher density. Supporting Figure 12: Effects of interactions U tRNA (r µ12 ij ) between tRNA genes and 5S rRNA genes are examined by comparing the differential probability distribution maps of lowly expressed genes. Differential probability density maps of (a) the bottom 200 genes in the wild-type (WT) simulated with U tRNA (r µ12 ij ) = 0, (b) the bottom 200 genes in the WT simulated with U tRNA (r µ12 ij ) = 0, (c) the bottom 200 genes in the yku70 esc1 mutant with U tRNA (r µ12 ij ) = 0, and the bottom 200 genes in the yku70 esc1 mutant with U tRNA (r µ12 ij ) = 0. The differential probability density is shown in gray scale; the lighter shade represents the higher density.
Supporting Radius of repulsion to confine a chain in nucleus a,g 800 nm u 1 Min. radius of attraction of telomere to envelope a,j,k 400 nm u 2 Max. radius of attraction of telomere to envelope a,j,k 800 nm u 3 Min. radius of repulsion of telomere from envelope a,j,k 950 nm c Strength of attraction of telomere from envelope a,h,j,k 3 a Values were chosen to reproduce the microscopic observation (11). b The same value as used by Zhang and Wolynes (2). c The value was chosen to reproduce the observed packing density of chromatins (8). d The same value as used by Rosa et al. (6). e k φ is defined in Eq.8 of Ref. (3) as an elastic constant of bending chain. Based on the assumption that chromatins are not the regular 30 nm fibers but are disordered and flexible, the value is set to k φ = 0 (i.e., the bending potential is not used) in the present model. f The value was chosen to allow enough sampling during the simulation. g Used the microscopic observation (12) as reference. h Used the microscopic observation (13) as reference. i Used the microscopic observation (14) as reference. j Used the microscopic observation (10) as reference.. k Used the microscopic observation (15) as reference.. Table 2: Averages and standard deviations of distributions of telomeretelomere distances. Distances between the left telomere of chr. III (3L) and the right telomere of chr. III (3R), the left telomere of chr. VI (6L) and the right telomere of chr. VI (6R), the left telomere of chr. V (5L) and the right telomere of chr. V (5R), the left telomere of chr. XIV (14L) and the right telomere of chr. XIV (14R), 6L and 14L, and 5L and 14R were examined in nm units. In simulation, 40 trajectories were calculated, and from each of them, 250 structures were sampled. d is the average of the simulated telomeretelomere distance in 250×40 = 10 4 structures, σ 0 is the average of the standard deviation for the distance in 250 structures that appeared in individual trajectories, and σ is the standard deviation for the distance in 10