article: Chromatin organization can be probed by Chromosomal Capture (5C) data, from which the encounter probability (EP) between genomic sites is presented in a large matrix. This matrix is averaged over a large cell population, revealing diagonal blocks called Topological Associating Domains (TADs) that represent a sub-chromatin organization. To study the relation between chromatin organization and gene regulation, we introduce a computational procedure to construct a bead-spring polymer model based on the EP matrix. The model permits exploring transient properties constrained by the statistics of the 5C data. To construct the polymer model, we proceed in two steps: first, we introduce a minimal number of random connectors inside restricted regions to account for diagonal blocks. Second, we account for long-range frequent specific genomic interactions. Using the constructed polymer, we compute the first encounter time distribution and the conditional probability of three key genomic sites. By simulating single particle trajectories of loci located on the constructed polymers from 5C data, we found a large variability of the anomalous exponent, used to interpret live cell imaging trajectories. The present polymer construction provides a generic tool to study steady-state and transient properties of chromatin constrained by some physical properties embedded in 5C data.
  Chromatin is organized in heterogeneous sub-regions of various sizes, as recently revealed by Chromosome Capture (5C) data [1,2]. This multiscale organization is generated by short and long-range genomic interactions between DNA segments, observed in the statistics of a large number of cells. Mammalian chromatin at a resolution of 3kb, [3,4] contains an organization at 1Mbp scale, where several sub-structures are enriched with intra-connectivity, reflecting an increased encounter probability (EP) between genomic segments. This increased EP is described in the two-dimensional encounter frequency (EF) matrix, containing diagonal blocks called Topologically Associating Domains (TADs) [3,5]. TADs are associated with gene regulation [3], DNA replication timing [6], DNA entanglement or cross-linking by molecules such as cohesin, CTCF [7] and condensin [3]. Cross-linking between chromatin sites are precisely the events sampled by Chromosome Capture data (3C, 4C, 5C, HiC) [3,8], and single cell HiC confirms that positions of cross-links can vary between cell types and phases [9]. In that context, TADs represent average chromatin conformations, characterized by a higher numbers of binding molecules compared to non-TAD regions. Over the past ten years, polymer models have been used to analyze statistics hidden in Chromosome Capture (CC) data and to characterize the decay of the encounter probability with the genomic distance s. For example, the EP between two monomers A, B for a (linear) Rouse polymer decays with s−3/2, thus the exponent is 3/2 [10]. A range of decay exponents lower than 3/2 can be produced by other polymer models, where stable transient loops are formed between segments [11]. By varying the number of loops, a large range of chromatin configurations can be generated and the associated polymer characteristics are reflected in the EP decay exponent. By including transcriptional information, dynamic-loop model [11,12] can reproduce chromatin looping associated with transcriptional activity. Similar polymer models describe chromosomal territories [13], suggesting that inactive genes are located inside these territories. The strings-and-binders-switch polymer model consists of reversible binding between specific genomic segments located in close proximity, when an additional diffusing molecule is present at the binding site. In an extension of this model, genomic segments having similar epigenomic state can directly interact [14], revealing the phase diagram of polymer configurations. These models are used to interpret the contact probability decay in the Hi-C data [15–18]. In a new class of polymer model, each monomer interacts with any other through a potential well [5,19], where pairwise interaction between monomers is represented by either an attractive or repulsive potential. The parameters of the model are extracted from data by minimizing the chi-square norm between the EP empirical and Monte-Carlo simulation matrices. However, this highly connected model does not translate easily into molecular binding because the nature of these potential wells does not have a direct physical interpretation. When this model is applied to the region containing the X inactivation center of mammalian embryonic stem cells, it predicted novel ensemble of polymer configurations, representing TAD structures present in the 5C data. Other minimization procedures were used to inter-chromosomal distances at a large-scale resolution of 1Mbp [20]. Polymer models have also been used to interpret single particle trajectories (SPTs) of tagged DNA locus [21–27] revealing that chromatin is constantly remodeled. SPTs are characterized by their anomalous behavior, which can deviate significantly from classical diffusion and are usually quantified by the mean-square displacement (MSD). In that context, a model with a minimal number of parameters is still needed to account for both types of data: 1) the EPs decay rate and 2) dynamical parameters such as the anomalous exponent extracted from SPTs. In the absence of a systematic procedures to convert the EP into a polymer model with a similar EP decay rate as the one presented in the 5C data, the connection between 1) and 2) was left open. We present here a general computational and algorithmic procedure to estimate parameters of a randomly cross-linked polymer model following the 5C protocol. The procedure consists in constructing an ensemble of polymer models from the EP of 5C by randomly cross-linking monomers and in resolving the difficulty of assigning the minimal number of sparse interactions between monomer pairs. These interactions can be directly interpreted as binding molecules. The construction of the polymer model starts with the Rouse model [10], which consists of beads linearly connected by harmonic spring. We started with the coarse-grained Rouse polymer that describes accurately the statistics of the chromatin below a scale of few Mbp [28,29]. To further constrain monomer interactions, we determine monomer connectivity from the 5C data of mammalian X chromosomes. The construction procedure is divided into two steps: first, to account for heterogeneity in the 5C data, we added a minimal required number of connectors (cross-links) between genomic sites chosen at random that can reproduce TAD blocks. We show that the number of random connectors to be added is uniquely determined from data. However, this step is insufficient to recover the EP decay peaks contained in the 5C data. Thus, in the second step, we account for consistent long-range interaction present in the EP matrix within and between TADs. We calibrate our model by requiring that the EP matrix, constructed from simulations of the polymer, has the same decay exponent (for each monomer) as that of the empirical data. These calibrated polymer models allow us to study transient properties and to estimate the conditional encounter probability and the first encounter time between three specific genomic sites. Although the generalized Gaussian polymer model we are using here is well known, our reconstruction using minimal number of short and long-range connectors derived from empirical EP is new and permits generating novel statistics describing transient gene regulation. By exploring the statistics of simulated SPTs (by computing the anomalous exponent [21–24]), we further show that the large heterogeneity of the anomalous exponent present in live cell imaging can be explained by random binding locations on the chromatin that can vary from cell to cell.
  We previously described how we construct a polymer model from a symmetrized 5C matrix M (Fig 1A). By symmetrizing the EP matrix, we averaged-out asymmetrical fluctuations. The 5C data we used represent a sub-region of the X chromosome (≈ 92kbp), that was previously segmented into two regions called Topological Associating Domains (TADs) D and E [3]. The matrix M was further coarse-grained by binning the encounter frequencies into 307 monomers of 3kbp [5], where TAD D (resp. TAD E) is represented by the first 106 monomers (resp. 107–307), as shown in Fig 1A. We introduce a general polymer model (Fig 1B) with arbitrary configuration, the properties of which will be extracted from empirical 5C matrix M. The encounter probabilities between monomer m and monomer n are computed from the experimental 5C matrix M (see Material and Methods) by P e (| m - n | | n) = M n, n + | m - n | + M n, n - | m - n | ∑ m = 1 N M n, m, (1) which depends on the genomic distance |m − n| (Fig 1C). Although the average encounter probabilities decays with |m − n| for each n, they contain peaks that reflect consistent long-range interactions between monomers. To quantify the decay of the EP, we fitted its average value P e (| m - n |) = 1 N ∑ n = 1 N P e (| m - n | | n) (black dotted line in Fig 1C) with the function P ˜ (| m - n |) = C | m - n | β, (2) where C and β > 0 are two constants. For a Rouse polymer, the EP function P ˜ is characterized by a decay exponent β = 3/2 [10]. Fitting Eq 2 to data, revealed that β = 0. 77, from which we concluded that the polymer model should be modified to account for higher compaction than allowed by a Rouse polymer [30]. To better account for the heterogeneity in the EP of each monomer, we plotted the distributions of the exponent βn for n = 1‥307 along the polymer (Fig 1D blue dots). The exponents βn were extracted by fitting the function 2 to the empirical EPs Eq 1. The large variability in βn, n = 1‥307 reflects the local heterogeneity of the chromatin architecture at the current scale (a monomer represents 3kbp). The average value β for TAD D and E was found to be βD = 0. 81 and βE = 0. 74 respectively. The local minima of β deviated significantly from the mean values (Fig 1D red squares), where the deviation is 2–3 times the standard deviation computed from the β values in each TAD. The minima may represent chromatin features (Fig 1E) such as specific long-range interactions or boundary between chromatin sub-domains. Indeed, point m100 (around monomers 102–107 in Fig 1D red) is located at the boundary between TAD D and E, while m24 and m162 are characterized by strong long-range interactions (Fig 1E). To conclude, the distribution of β values extracted from the EP is quite heterogeneous, which can disclose chromatin subregions and long-range strong interactions. We shall account in the next two sections for these characteristic features and include in the polymer model both random and persistent long-range connections between monomers. To determine the level of connectivity of the generalized Rouse polymer which reproduce the EP-decay using a prescribed exponent β, we first studied the case of one TAD-like region using the simulation of a 307 monomer chain. For each realization, we added connectors between random non-nearest neighbor monomer-pairs in the subregion 103–203 (Fig 2A). The number of connectors, or the connectivity percentage ξ (fraction of the number of connected monomer-pairs to the maximum, described in Materials and Methods), was increased in the range 0–2%. By adding connectors, the EP between distant monomers has increased, as presented in the EP-matrix (Fig 2B). In contrast, outside the region 103–203 the EPs were similar to the case ξ = 0 (linear chain), showing that the connected region did not affect the EP in the non-connected ones. At this stage, we have shown that adding random connectors allows recovering the shape of TAD regions. To find the minimal numbers of connectors necessary to recover a given TAD, we aim of elucidating the relationship between the connectivity percentage ξ and the decay exponent β. For that purpose, we simulated an ensemble of polymers to their relaxation time (Materials and Methods), and used the equilibrium configuration to estimate the EP of each monomer for ξ ∈ [0,2]%. We calculated the exponent β by fitting the function 2 to the simulated encounter probability data: the values of βn for n ∈ [103,203] decreased with ξ. Indeed, for ξ ≈ 0. 2%, the coefficients βn decreases below the Rouse exponent (equals to βRouse = 1. 5), indicating compact polymer configurations. For ξ = 2%, the mean decay exponent βn and n = 103–203 was β ¯ = 0. 47, with a minimal value 0. 42 obtained for the boundary monomers 103 and 203 (Fig 2C). These results confirm that polymer condenses when ξ increases, as confirmed by computing the mean radius of gyration Rg [10], showing a decay from Rg = 11. 9 for ξ = 0 to Rg = 11 at ξ = 2% (measured for 307 monomers, Fig 2D). Values of β outside the TAD region (monomers 1 to 102 and 204–307) were mostly unchanged, fluctuating around β = 1. 5, confirming that statistical properties of a Rouse chain are unaltered when connectors are added to the middle region (n = 103 – 203). Finally, the average value of βn (computed over n = 103 – 203) versus ξ is shown in Fig 2D and, as we shall see, will serve to extract the connectivity percentage ξ from the empirical data. To reproduce the two TADs D and E of the X-chromosome, we started with a polymer of 307 monomers and added connectors randomly (green arrows) between monomers 1–106 and between monomers 107–307, as described in Fig 2E upper panel). This partition follows the empirical TAD segmentation described in [3,5] at the scale 3kbp of the polymer model. Three polymer realizations for ξ = 0,0. 2,1% are shown in Fig 2E bottom panel, showing polymer condensation into two distinct regions. The EF matrix shows that two TAD-like regions, named TAD1 and TAD2 (Fig 2F), emerge as the connectivity ξ increases from 0 to 2%. To extract the exponent β (Fig 2G, colored curves) we fitted the function 2 to the EP matrix for each monomer inside TAD1 and 2. In both cases, the exponent β decreased below βRouse = 1. 5 (ξ = 0% blue curve) and for ξ = 0. 2%, 11 and 39 random connectors were added for TAD1 and TAD2, respectively. The boundary between TADs is characterized by an abrupt decay of the β value, reflecting high long and short-range encounters. The average β exponent (averaged over each TAD), plotted with respect to the connectivity ξ (Fig 2H), was used to determine the number of connectors necessary to reconstruct the empirical data. Indeed, we extracted from the EP matrix (Fig 1A) that βD = 0. 81, βE = 0. 78 the associated connectivity percentages ξD = 0. 12, ξE = 0. 23, respectively (Fig 2H). To conclude, we obtain the minimal number of random connectors to be added on a generalized Rouse polymer such that the decay exponents of the reconstructed and empirical EP-matrix are as close as possible. A key feature present in the 5C EF-matrix (Fig 1A) is the ensemble of persistent long-range interactions between monomers (Fig 1C). To account for these interactions, we connected monomers corresponding to off-diagonal local maxima of the EF matrix for which their EP exceeds that of nearest neighboring monomers threshold (Materials and Methods). We found 24 long-range connections: 7 (resp. 13) within TAD D (resp. E) and 4 across the two (see SI for the list of monomers pairs) as shown in Fig 3A. We adjusted the spring constant κm, n between monomer m and n corresponding to consistent long-range interactions, based on the values of their EPs (Materials and Methods). The scaled coefficient κm, n are summarized in the SI, and were found to be 1. 1 to 3 times higher than the spring constant assigned to connectors of the linear backbone (κ = 0. 97). As revealed by simulations of a Rouse polymer with fixed long-range connectors, the polymer configuration are condensed, characterized by a radius of gyration of about Rg = 9. 1 (compared to Rg = 12 for the Rouse chain). Three realizations of the polymer are shown in Fig 3A. To quantify the effect of adding consistent long-range connections on the EP decay, we simulated a Rouse chain containing 307 monomers with the addition of fixed monomer connectivity extracted from the peaks of the empirical EP-matrix. We computed the decay exponents βn of each monomer by fitting the function 2 to the EPs from simulations, and compared it to the ones computed from the experimental data (Fig 3B). To estimate the quality of the reconstruction we use the L1-norm, defined for a function f by |f|1 = ∑k|f (k) |, and computed the difference between the experimental βExp and simulated βSim curves normalized by the norm ∣βExp∥1 and we find ∥ β E x p - β S i m ∥ 1 | β E x p ∥ 1 = 0. 312. (3) The experimental values βn (Fig 3B red) were generally lower than the ones obtained from simulations (blue), indicating that the reconstructed chromatin polymer is less condensed for both TADs. The mean β values for TAD D and E were quite similar with βD = 1. 02 and βE = 0. 99. Therefore, we conclude that long-range connectors are insufficient to reproduce the statistics of the 5C data. We previously evaluated separately the effect of adding random connectors and fixed long-range interactions on the EPs. We computed the decay exponent β and compared it with coarse-grained 5C data. We now combine these two constraints, such that specific and non-specific connectors are added to each realization of a generalized Rouse polymer (Fig 4). We first find the connectivity percentage matching that of the experimental data in each TAD ξD = 0. 23% and ξE = 0. 12%. These values summarize the contribution from the two types of connectors (Fig 2H). For long-range specific interactions (Fig 3), we previously obtained βD = 1. 02, βE = 0. 99, corresponding to ξD = 0. 12%, ξE = 0. 07% (Fig 2H) for TAD D and E, respectively. Therefore, we attributed the remaining percentages to the addition of random connectors, that is ξD = 0. 11, ξE = 0. 05. The number of added random connectors corresponding to ξD = 0. 11%, ξE = 0. 05% are 6 and 10 in TAD D and E, respectively. To reconstruct the polymer, we started with a Rouse chain (Fig 4A (gray) ) and added 24 connectors between monomer pairs corresponding to peaks of the EP matrix (red). Three polymer realizations, simulated with the two types of connectors, are shown in Fig 4B, characterized by a radius of Gyration Rg = 6. 4. We computed the EF-matrix (Fig 4C) that showed similarity with the experimental data (compare Fig 1A with Fig 4B), for which two TAD-like structures are visible. We find a satisfactory agreement between simulations and experimental data, measured by the Kolmogorov-Smirnov distance DMax = 0. 06 computed on the cumulative density function between the reconstructed and experimental data, as shown in S1 Fig). We further quantified the similarity between the two matrices by comparing the decay exponent for each monomer βn, (n = 1‥307) from numerical simulations to those of the experimental data. We used the function 2 to fit the EP of monomers 1–307 after long time polymer simulations (Materials and Methods). The fitted value for β shows an good agreement with the experimental β values (Fig 3C). Comparing the normalized difference between β signals we find ‖βExp − βSim‖1/‖βExp‖1 = 0. 2 (compared with 0. 312 Fig 3B). To conclude, accounting for long-range deterministic and short-range stochastic interactions leads to a more accurate polymer model for chromatin reconstruction. The quality of this approximation is measured by the decay norm of the β exponent between the data and the simulations, also confirmed by the Kolmogorov-Smirnov distance as shown in S1C Fig and to be compared with S1D Fig. We showed previously how to construct a polymer, the statistical properties of which match the ones extracted from 5C data. However, the 5C data cannot be used to study transient properties of the chromatin, as they represent a static genomic encounter interactions averaged over cell population. We shall now use the calibrated polymer described above, constrained by the steady-state properties of the 5C data (Fig 5A), to evaluate transient properties of the chromatin. We focus here on the 5C data harboring the Xist loci, which is the master regulator of X chromosome inactivation (XCI), and its antisense transcript, Tsix, which plays a key role in modulating Xist expression during mouse development [3,5]. Tsix is believed to play a key role in the choice of the Xist allele that will be expressed during X inactivation. Thus, we decide to estimate how chromatin conformation within TADs can contribute to this transcriptional variability using the present polymer model. We estimated the first encounter time distribution and the probability that monomer 26 (position of the Linx) meets monomer 87 (Xite) before monomer 62 (Chic1). These monomers represent three key sites on the X chromosome [3,5], located in TAD D. We show three realizations in Fig 5A and indicate the location of the three sites inside TAD D (yellow) and E (blue). We started the polymer simulations from the steady-state distribution and performed around 10,000 runs. As predicted by the narrow escape theory [31,32], the encounter time between two of the three monomers is Poissonian, and we confirmed this result by simulating the distributions (Fig 5B). The reciprocal of the mean encounter time is by definition [33] the encounter rate that we extract in Fig 5B–5D. We found also that the encounter probability, computed from the encounter rates (rate divided by the sum of the rates) as between Linx and Chic1 is P = 0. 55, while the mean encounter times between each pair (Linx-Chic1) and (Xist-Linx) are comparable of the order of 131s (table C in Fig 5). Finally, to check the impact of TAD E on the encounter time inside TAD D, we ran another set of stochastic simulations after removing TAD E (Fig 5D). Surprisingly, the encounter probability was inverted compared to the case of no deletion, while the mean time was increased by almost 50% to 195s (Linx to Chic1) and 205s (Linx to Xite). This result suggests that specific long-range interactions between TAD D and E (table in S1 Text) serve to modulate the probability and the encounter time between the three key genomic sites. This result further indicates that the search time inside a TAD can be influenced by neighboring chromatin configurations. To further study the statistical properties of loci trajectories, we simulated the three loci monomer 26 (Linx), monomer 87 (Xite) and 62 (Chic1), following classical single particle tracking experiments (Fig 6A). We used the calibrated polymer model reconstructed in section Encounter probabilities and distribution of search times of three genomic sites. Starting from an equilibrium configuration following a relaxation time, we estimated the mean-square displacement (MSD) and computed the anomalous exponent over all realizations. The MSD of a stochastic process X (t) is computed by averaging over trajectory Xi (t) realizations [33], ⟨ | X (t) - X (0) | 2 ⟩ = lim N p → ∞ 1 N p ∑ i = 1 N p | X i (t) - X i (0) | 2 (4) and for short time, we use the asymptotic behavior ⟨ | X (t) - X (0) | 2 ⟩ ∼ A t α, (5) where A is a constant and α is called the anomalous exponent. In practice, we computed the MSD from the estimator 1 N p ∑ i = 1 N p | X i (t) - X i (0) | 2, where the number of trajectories Np is of the order 1000. We fitted the function Atα to the MSD Eq 5, computed on the simulated trajectories and then extracted α. We explore now the consequence of computing the MSD when averaging over a cell population. To reproduce this situation, we consider a class of polymer configurations C6,10 obtained by fixing the number of random connectors 6, and 10 for TAD D and E respectively. However, their positions in each realization is allowed to vary, generating a variability. We shall now estimate the anomalous exponent αw of a locus for a given configuration w ∈ C6,10 and then compute the average <αw>w ∈ C6,10 over the ensemble C6,10. We shall emphasize that long-range specific connectors (table in S1 Text) are also accounted for in the polymer ensemble C6,10. For two specific configurations from C6,10, one with TADs E + D (Fig 6A Upper) and the other TAD D alone (Fig 6A Lower), we show SPTs for each monomer m = 26,62,87, projected in two dimensions. We then computed the MSDs curves by averaging over all realization (Fig 6B). Further, the distribution of the anomalous exponents is shown in Fig 6C. Interestingly, we found that each locus m = 26,62,87 had a different mean anomalous exponent α 26 D + E = 0. 39, α 62 D + E = 0. 4, α 87 D + E = 0. 31, revealing the intrinsic heterogeneity present in the chromatin. The α values we have computed for all monomers are indeed heterogeneous, as can be seen in histograms of Fig 6B (right column), 6C and S2A–S2B Fig left columns. Values are spreading in the range 0. 3 to 0. 5, however there is a peak in the histogram between values 0. 4 to 0. 45. To evaluate how the anomalous exponent is influenced by persistent long-range interactions present in the EP data, we repeated polymer simulations and MSD analysis by removing TAD E. The distribution of the anomalous exponent α shows that removing TAD E results in the loss of low values, which characterize high chromatin connectivity (Fig 6B, right). We found the following changes in the anomalous exponents: α 26 D = 0. 46, α 62 D = 0. 4, α 87 D = 0. 43. This result shows how specific long-range interactions affect the local chromatin dynamics and locus motion. Finally, the anomalous exponent averaged over all monomers is αD+E = 0. 425, whereas for TAD D and E we found αD = 0. 433, αE = 0. 42. However, when TAD E is removed, the average anomalous exponent for all monomers in TAD D is αD = 0. 458. To represent the contribution of a single cell experiments in a population, we simulated SPTs for a fixed polymer configuration (we chose an ensemble of connectors in C6,10). We computed the MSD of all loci n = 1‥307 that can be divided into 3 classes: low medium and high anomalous exponent (SI). The distribution of the anomalous exponents of all sites is quite uniform (S2 Fig). We then vary the connector positions and computed the spread of the anomalous exponents (Fig 6C). This result shows that changing the connectors position account for the variability of the anomalous exponents. This result clarifies how the local chromatin organization affects the MSD when computed across cell population [21]. Indeed, random connectors model molecular binding that can vary from cell-to-cell. To conclude, the construction of polymer models from 5C data can now be used to simulate SPTs of any loci of interest and thus to explore the inherent statistical variability found experimentally (Fig 6C). The present results can be used to interpret the variability of the MSD found in SPTs of live cell imaging, where the same locus in different cells can exhibit a different anomalous exponent, depending not only on the locus position, but also on the intrinsic variability due to the random arrangements of binding molecules between cells.
  We presented here a general method to construct a coarse-grained polymer model from the 5C encounter probability (EP) matrix. This construction preserves the statistical properties of the 5C data, such as the decay rate of the EP of each monomer. The present approach is not used to study the configuration space of chromatin geometry, because it is too large to be fully sampled by elementary polymer models. However, we used this approach to generate statistics of passage times and radius of gyration, which characterize more accurately chromatin dynamics and are not contained in the 5C data. We built here a coarse-grained polymer model of the chromatin based on Rouse and we disregarded the repulsion forces between monomers and possible cross over of bonds that certainly influences the dynamics of the chromatin and statistical results. Future models should clearly examine the effect of repulsion forces on the present reconstruction, especially to study refined spatial scales below few kbp. We added connectors between random monomer-pairs to characterize sub-configurations present in 5C data. Connectors are represented by springs between monomer-pairs and account for the chromatin architectures. By adding connectors between random monomers, we were able to recover TAD sub-regions. By randomizing connectors positions, we could reproduce the inherent variability in chromatin architecture of nuclei population captured in 5C experiments. Using our methodology, the characteristics of the reconstructed polymer are derived directly from the empirical data and do not require any minimization procedure [5]. One of the key result here is to determine the number of connectors directly from the experimental EP decay (Fig 2). Connectors can directly be interpreted as molecular interactions mediated by proteins such as cohesin, condensing and CTCF bindings. For example, cohesin could bind at random places scattered along the chromatin at 5C data acquisition time. These bounds could generate TADs, as shown here using simulations (Fig 2). 5C Contact maps represent steady-state distribution obtained from looping events in large ensemble of millions of cells, where TAD structures appears. The present approach differs from classical reconstruction methods, where 3D structures of a genome are inferred from 5C contact frequency data [34–37]. Previous models explored the effect of connectors between regions of the chromatin [18,38] and examined the consequences on the EP-decay rate, but the positions and the number of these connectors were not derived from data. Here, we use connectors to resolve a reverse engineering problem, which is to recover the degree of connectivity from the EP-decay rate (Fig 2). The relation between the mean number of connectors and the decay exponent β of the EP (Eq 2) is found using simulations in section Encounter probability in the random loop polymer model. The decay exponent of the EP characterizes the polymer scaling statistics [38] and for this reason, the present model extends the key switch-and-binders model developed in [38]. It was surprising to find that a TAD subregion could influence the encounter distribution between monomers located in different TADs (Fig 5). Indeed, the distribution of looping time in free space depends only on the distance between monomers, while in confined domain, the nuclear boundary has an effect [30,32,39]. We found here that this modulation of loop regulation is due to long-range inter-TAD interactions, that are present in the 5C data and accounted for here in the construction of the polymer model. Indeed, we reported (SI) significant inter-TAD interactions between three monomer pairs: 86–234,86–260,24–285, where spring constants, representing long-range interactions are at least twice larger than other interactions. The threshold procedure described in Material and Methods (Eq 11) disregards the peaks of the 5C data for which the EP falls below the nearest-neighbor threshold Tth. The peaks in the 5C below Tth were ignored because we interpret them as transient events, or a statistic that was not shared by the majority of the chromatin events or it could simply be due to some random fluctuations in the data. Long-range stable interactions have also been found in other mammalian chromosomes [4]. Finally, the present polymer reconstruction allows probing the dynamics of single particle trajectories (SPTs). We use our reconstructed polymer model to explore transient properties and the encounter time distribution between any two sites. Additionally, we simulated single particle trajectories from a reconstructed polymer model (Fig 4) and estimated the anomalous exponents, following the routine procedure in experimental SPTs. The present analysis reveals that the variability of the anomalous exponent of a given loci is due to the heterogeneity of the local polymer configurations constructed from the 5C data. This heterogeneity is simulated here as the random cross-links between monomer pairs. Indeed, for a Rouse polymer or any uniformly connected polymers, the anomalous exponent is constant [30,40,41]. We suggest that the inherent fluctuations of the chromatin interactions can be due to the random positions of binding molecules. Indeed, even for a fixed number of connectors, there still remains an intrinsic variability of monomer binding positions, leading to the same EP decay rate. This structural heterogeneity originating from connector positions certainly influences the anomalous exponents and can be used to re-interpret experimental SPTs [19,21,42,43]. The present method and algorithms are generic and can be used to reconstruct a polymer model at a given scale (number of monomers and number of bps coarse-grained in a monomer) in a limit of few Mbp. Analyzing chromatin condensation and its transient properties can now integrate chromosomal capture data and SPTs statistics.
  The Rouse model [10] describes a polymer as a collection of beads Rn (n = 1…N) linearly connected by harmonic springs and driven by Brownian motion. The energy of the polymer is given by [10] ϕ R o u s e (R 1, ‥, R N) = 1 2 ∑ j = 1 N - 1 κ (R j - R j + 1) 2, (6) where κ = 3 k B T γ b 2 and b is the standard deviation of the distance between adjacent monomers, γ is the friction coefficient, kB the Boltzmann coefficient, and T the temperature. To account for a sub-chromatin region C N, characterized by a higher EP than the rest, we will add connectors between monomer-pairs chosen randomly (with uniform distribution) inside this subregion such that an additional potential ϕ R a n d (R 1, ‥, R N) = 1 2 ∑ j, k ∈ C N κ (R j - R k) 2, (7) is added to ϕRouse, where C N is the ensemble of indices defining the sub-region. The number of connectors is a free parameter, that will be determined from experimental 5C data. In addition, to account for consistent long-range interactions, reflected by peaks in EP matrix (Fig 1C), we will fix connectors between monomer-pairs by adding a spring constant κm, n between monomer m and n, so that the associated energy related to the peak interactions is described by ϕ P e a k s (R 1, …, R N) = 1 2 ∑ n, m ∈ S M a x κ n, m (R n - R m) 2. (8) We will discuss below how the spring constant κm, n is computed from empirical data. In summary, the total energy of a polymer containing random connectors and prescribed peaks, is the sum of three energies Eqs 6,7 and 8: Φ (R 1, ‥, R N) = ϕ R a n d (R 1, ‥, R N) + ϕ P e a k s (R 1, …, R N) + ϕ R o u s e (R 1, ‥, R N) (9) and the stochastic equation of motion for n = 1, ‥, N is d R n d t = - ∇ R n Φ (R 1, ‥, R N) + 2 D d ω n d t, (10) where D = k B T γ is the diffusion constant, γ is the friction coefficient, and ωn are independent 3-dimensional Gaussian noise with mean 0 and standard deviation 1. To account for the 5C-data, comprised of a subsection of the X-chromosome from female mice embryonic stem cells reported in [3], showing TAD D and E as two diagonal blocks (Fig 1A), we use the coarse-graining procedure of [5], with a polymer of length N = 307. Each monomer represents a genomic segment of 3kbp and is connected to its 2 nearest neighbors by a harmonic spring (see subsection above). TAD D (resp. E) is represented by the range of monomers from 1 to ND = 106 (resp. 107 to 307), and NE = 201 for TAD E. To reproduce the empirical EP extracted from data (see formula 1), we connected non-nearest neighbor pairs of monomers chosen randomly with uniform probability 1 N 1 (resp. 1 N 2) ) where N1 = (ND − 2) (ND − 1) /2 (resp. N2 = (NE − 2) (NE − 1) /2). The number of connectors in each TAD is a percentage ξ of the total possible number of non nearest neighbor pairs, for TAD D (resp. E), we have C D = ξ D N 1 100 (resp. E C E = ξ E N 2 100) and ξD, ξE ∈ [0,100], which will be extracted from data. Random connectors were not added between monomers belonging to different TADs. The procedure of adding random loops to a Rouse polymer is implemented using the energy of random loops, as described in the previous subsection. Finally, 24 connectors were added (SI) in all polymers, corresponding to the selected peaks present in the EP matrix, obtained in the following procedure: we located the positions of peaks that form a subset SMax of the ensemble of the off-diagonal local maxima in the EP-matrix, such that their EP is higher than a threshold Tth. In practice, we assumed that the encounter probability between neighboring monomers is almost not affected by the global chromatin structure. Therefore, any EP value in the matrix Mi, j above the threshold Tth (equals to the the EP of the nearest neighbor monomers) is considered to be a stable loop. The threshold Tth is defined as follows: T t h = ∑ i M i i - 1 + M i i + 1 ∑ i, j M i j. (11) In a second step, we determine the spring constants κm, n between monomer m and n in the ensemble SMax from the empirical EP Pm, n at each peak position. We recall that for a Rouse chain the joint probability density function of beads Rm and Rn is given by [10, p. 15] Φ (R m, R n, Δ m, n) = 3 2 π b 2 Δ m, n 3 / 2 exp - 3 (R m - R n) 2) 2 b 2 Δ m, n (12) where Δm, n = |m − n|. We assume that the encounter probability between neighboring monomers is not affected by global polymer structure, thus for the nearest neighbors Δm, n = 1 the EP occurs at small distances such that the exponential term is almost 1, that is P m, n = 3 2 π b 2 3 / 2 ≈ κ m, m + 1 2 π 3 / 2, (13) We approximate the chromatin structure as a polymer chain with a variance b2 between adjacent monomers. Thus, the constant b ¯ is estimated as the mean EP over the sub- and super-diagonals: 3 2 π b ¯ 2 3 / 2 ≈ ∑ i (P i i - 1 + P i i + 1) = ∑ i M i i - 1 + M i i + 1 ∑ i, j M i j = T t h (14) To account for long-range interactions, we applied formula 13 to estimate the effective spring constant from the empirical EP P ˜, so that κ m, n = 2 π P ˜ m, n 2 / 3. (15) Finally, the energy related to peak interactions is described by ϕ P e a k s (R 1, …, R N) = 1 2 ∑ n, m ∈ S M a x κ n, m (R n - R m) 2. (16) Using the reconstruction method described in the two previous subsections, we generate polymer realizations, each differs in the position of random connectors inside a TAD. To generate statistics from the EP matrix, we started from an initial random walk configuration and simulated the polymer until its equilibrium after a relaxation time τR. The time τR is determined for each realization from the formula τR = 1/ (κminλ1), where κmin is the minimal positive spring constant, and λ1 is the smallest non-vanishing eigenvalue of the polymer’s connectivity matrix [44], which we computed numerically (in practice it is of the order of thousands of simulation steps with Δt = 0. 01s). For the numerical simulations, we divide Eq 10 by D and the spring constants are scaled by the friction coefficient such that κ = d k B T γ b 2 = d D b 2, in dimension d = 3. The encounter frequency matrix of the 307 monomers is computed at time τR, where two monomers are considered to have encountered if their distance is less than ϵ. The time step for all simulations is Δt = 10−2 s. The value of the parameter b is computed using formula 14. Using the b value in the relation κ = 3kBT/b2 [10], we computed the spring constant for random connectors and the linear backbone of the polymer. The spring constants for persistent long-range connectors were estimated using Eq 15 and are summarized in the table in S1 Text, other simulation parameters are summarized in Table 1.
summary: Chromatin organization remains poorly understood and polymer models are used to reconstruct such organization, to reveal hidden structures and to quantify genomic interactions. We use a generalized Rouse model (a linear chain of beads connected by springs) with additional interacting molecules that allow stable loop formation. The polymer models are constructed using the minimal number of binding molecules, positioned according to the encounter probability matrix obtained from experimental chromosomal capture data. We determine the conditional encounter probability of 3 key loci regulating gene inactivation from our calibrated polymer model. Using polymer simulations, we generate single particle trajectories and explore their transient properties. The present results suggest that the heterogeneity of anomalous exponents measured in live cell imaging is due to the large combinatorics in reconstructing the chromatin organization from Chromosomal Capture data. The present method and algorithms are generic and can be used to reconstruct a polymer model at a given scale from any Chromosomal Capture data.
keywords: gene regulation
  simulation and modeling
  probability distribution
  mathematics
  materials science
  epigenetics
  macromolecules
  chromatin
  materials by structure
  research and analysis methods
  polymers
  polymer chemistry
  sex chromosomes
  genomics
  chromosome biology
  gene expression
  chemistry
  x chromosomes
  probability theory
  cell biology
  genetics
  biology and life sciences
  physical sciences
  computational biology
  genomics statistics
  chromosomes
title: Transient chromatin properties revealed by polymer models and stochastic simulations constructed from Chromosomal Capture data