{'fraction_non_alphanumeric': 0.04123905470294699, 'fraction_numerical': 0.012710667545428867, 'mean_word_length': 4.4901271580688515, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, '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}}
Ubiquitous nucleosome unwrapping in the yeast genome

Nucleosome core particle is a dynamic structure -- DNA may transiently peel off the histone octamer surface due to thermal fluctuations or the action of chromatin remodeling enzymes. Partial DNA unwrapping enables easier access of DNA-binding factors to their target sites and thus may provide a dominant pathway for effecting rapid and robust access to DNA packaged into chromatin. Indeed, a recent high-resolution map of distances between neighboring nucleosome dyads in \emph{S.cerevisiae} shows that at least 38.7\% of all nucleosomes are partially unwrapped. The extent of unwrapping follows a stereotypical pattern in the vicinity of genes, reminiscent of the canonical pattern of nucleosome occupancy in which nucleosomes are depleted over promoters and well-positioned over coding regions. To explain these observations, we developed a biophysical model which employs a 10-11 base pair periodic nucleosome energy profile. The profile, based on the pattern of histone-DNA contacts in nucleosome crystal structures and the idea of linker length discretization, accounts for both nucleosome unwrapping and higher-order chromatin structure. Our model reproduces the observed genome-wide distribution of inter-dyad distances, and accounts for patterns of nucleosome occupancy and unwrapping around coding regions. At the same time, our approach explains \emph{in vitro} measurements of accessibility of nucleosome-covered binding sites, and of nucleosome-induced cooperativity between DNA-binding factors. We are able to rule out several alternative scenarios of nucleosome unwrapping as inconsistent with the genomic data.

our framework allows us to predict both sequence-specific free energies of nucleosome formation and the unwrapping potential from paired-end high-throughput sequencing datasets.

Nucleosome unwrapping potential
We use a high-resolution in vivo map of nucleosome dyad positions based on chemical modification of engineered histones and DNA backbone cleavage by hydroxyl radicals [30]. With cleaved DNA segments sorted by molecular mass on the agarose gel and the ∼ 150 − 200 bp fraction sequenced using pairedend reads, the map provides a direct measurement of both dyad positions and distances between adjacent dyads. Although superior to methods based on micrococcal nuclease (MNase) digestion whose accuracy is affected by MNase sequences preferences and its tendency to over-or under-digest DNA [33,34,35], the map is biased by unknown hydroxyl radical cutting preferences for two alternate sites at each DNA strand [30].
A genome-wide histogram of inter-dyad distances shows that at least 38.7 % of all nucleosomes are partially unwrapped (blue line in Fig. 1C). In order to study the energetics of unwrapping, we introduce a simple model based on the 10-11 bp periodic pattern of histone contacts with the minor groove of the nucleosomal DNA [31,32] (Fig. 1B). As DNA is peeled off each contact patch, its free energy increases because hydrogen bonds and favorable electrostatic contacts between histone side chains and the DNA phosphate backbone are lost. However, once DNA breaks free from the contact patch, it may adopt multiple conformations, which allows it to increase its entropy and thus lower the total free energy. The favorable entropic term grows with the extent of unwrapping until the next contact patch is reached, completing one cycle in the oscillatory energy profile. The oscillations are superimposed on a straight line whose slope reflects the average free energy cost per bp of histone-DNA contact loss. Additional details of the potential construction can be found in SI Appendix, Model A. The histone-DNA potential constructed in this way has no sequence specificity, in contrast to the free energy cost of bending DNA into a nucleosomal superhelix [36]. To a good approximation, we expect sequence effects to average out in the genome-wide histogram of inter-dyad distances (this assumption is tested later).
Thus we aim to reproduce the observed distribution of inter-dyad distances with a model in which nucleosome energetics is sequence-independent but transient unwrapping is allowed. We compute the conditional probability P (c + d|c) of finding the nucleosome dyad at bp c + d, given that the previous dyad is located at bp c (Materials and Methods). Since inter-dyad distances cannot be used to distinguish between symmetric and asymmetric unwrapping, we assume the former for simplicity. The model is fit to the observed distribution of inter-dyad distances (SI Appendix). The free parameters of the model in-clude the amplitude of the oscillations, the slope of the free energy profile and a min(max) , the minimum (maximum) effective length of the nucleosome particle (SI Appendix, Model A). The maximum extent of nucleosome unwrapping is controlled by a min , while a max is allowed to exceed 147 bp in order to account for the effects of higher-order chromatin structure and linker histone deposition. We also fit the relative frequency of hydroxyl radical DNA cleavage at the −1 position with respect to the nucleosome dyad (SI Appendix) and the chemical potential of histone octamers. Our model reproduces both the overall shape and fine oscillatory structure of the observed inter-dyad distance distribution (Fig. 1C,D). In contrast, models without unwrapping are unable to capture even the overall shape of the observed inter-dyad distribution (gray line in Fig. 1C).

Higher-order chromatin structure and linker histone energetics
The effective length of the particle found in the fit, a max = 163 bp, is greater than 147 bp, the length of the DNA in the nucleosome core [32]. Indeed, a max = 147 bp is incompatible with the observed inter-dyad distribution (Fig. S3A,B). The model is less sensitive to the value of a min because extensively unwrapped particles are energetically unfavorable and therefore are not frequently seen in the data. The overall shape of the inter-dyad distribution is also sensitive to the slope of the energy profile in Fig. 1B, providing a robust estimate of the average nucleosome unwrapping energy (Fig. S3C). The fitted value of the slope yields 14.4 k B T for the average histone-DNA interaction energy in a fully wrapped nucleosome. Thus the energy profile in Fig. 1B describes both DNA interactions with the histone octamer core (up to 73 bp from the dyad) and the effects of higherorder chromatin structure, including, potentially, the attachment of Hho1p, the H1 linker histone of S. cerevisiae, to the DNA immediately outside of the nucleosome core [37,38,39]. Although Hho1p is less abundant in yeast than in higher eukaryotes, it is involved in higher-order chromatin organization, including chromatin compaction in stationary phase [40,39]. Relatively little is known about the molecular mechanism of H1 binding: there is no consensus yet whether the binding is symmetric or asymmetric, or even what the extent of the H1 footprint is [37,38]. H1 binding and other factors that mediate chromatin folding into higher-order structures cause linker lengths to be discretized [1,41]. Linker length discretization can be described by a periodic decaying two-body effective potential between neighboring nucleosomes, with the first minimum approximately 5 bp away from the nucleosome edge [1,42,43].
Based on these observations, we have constructed two models for the energy profile outside of the nucleosome core region. One model is a polynomial fit that extends the quasiperiodic profile of the unwrapping energy through another cycle ( Fig. 1B; SI Appendix, Model A). The depth and the position of the first minimum outside of the nucleosome core are additional free parameters. As can be seen in Fig. S4B, our fit robustly predicts the first minimum to be positioned 5-6 bp outside of the nucleosome core, in agreement with previous 0.5 1 Figure 1: Genome-wide distribution of inter-dyad distances. (A) A nucleosome is a dynamic structure. Transient nucleosome unwrapping followed by factor binding prevents subsequent rewrapping, mediating further binding events. (B) Nucleosome energy profile. The energy of a nucleosome that covers 2x + 1 bps is given by u SI nuc = 2u half (x) with symmetric unwrapping. The minima and maxima of the energy landscape are based on the crystal structure of the nucleosome core particle [31,32]. Dark gray bars show where the histone binding motifs interact with the DNA minor groove in the structure. Light gray bars show where the DNA major groove faces the histones. The energy profile was obtained by a polynomial fit as described in SI Appendix, Model A. (C) The inter-dyad distance distribution from a high-resolution nucleosome map [30] (blue), and from the model with (red) and without unwrapping (gray). In the model without unwrapping, a min = a max = 147 bp and the fitting parameters are E b , µ and f (SI Appendix, Model A). RMS -total root-mean-square deviation between the model and the data. (D) Oscillations in the observed (blue) and predicted (red) inter-dyad distance distributions, obtained by subtracting a smooth background from the data and the model with unwrapping in (C). The smooth background was found by applying a Savitzky-Golay filter of polynomial order 3 with 31 bp length (using the sgolayfilt function from the Signal Processing Toolbox of MATLAB). Correlation refers to r osc , the linear correlation coefficient between measured and predicted oscillations.
studies [1,42,43]. The depth of this minimum is comparable to the depth of the unwrapping minima (Fig. S4, SI Appendix, Model A).
The other model represents the energy profile outside of the nucleosome core by a linear function ( Fig. S5A; SI Appendix, Model B). The two free parameters are the slope and the length of the line, which are related to the H1-DNA interaction energy and the H1 footprint. This model is motivated by a picture of the H1 histone being gradually detached from its DNA site immediately outside of the nucleosome core. This alternative scenario, although likely oversimplified, can be used to check the sensitivity of our results toward a particular energy profile outside of the core region. We find that the linear profile fits the overall shape of the inter-dyad distribution somewhat less well than the oscillatory one (compare RMS values in Figs. 1B and S5B), although the 10-11 bp periodic fine structure is reproduced in both cases (Figs. 1C, S5C). The optimal linear profile is 7 bp long, yielding a symmetric H1 footprint with two 7 bp half-sites ( Fig. S5D) and the H1-DNA interaction energy of ≈ 5k B T (Fig. S5E).

Alternative models of nucleosome unwrapping
Next we have tested the sensitivity of our fits to the functional form of the unwrapping free energy profile. Although our primary model follows nucleosome crystal structures in creating a quasi-periodic energy profile with both 10 and 11 bp modes, strictly periodic 10 or 11 bp sinusoidal profiles yield nearly the same quality of fit ( Fig. S6; SI Appendix, Models C,D). Since the initial phase of the oscillations is not determined by the crystal structure anymore, it becomes another fitting parameter. The fitted initial phases in the 10 and 11-bp models make the periodic curves match the crystal structure further away from the dyad, where most of the observed unwrapping takes place (Fig. S6A). The phases diverge closer to the dyad, where they are not as strongly constrained by the data. RMS deviation is less sensitive to the initial phase than r osc , the linear correlation between predicted and observed oscillations in the inter-dyad histograms (Fig. S7). The primary peak in the dependence of r osc on the initial phase matches the crystal structure. There is also a secondary peak corresponding to the 5 bp shift in the unwrapping energy profile, which in turn leads to the 10 bp, in-phase shift in the distribution of inter-dyad oscillations (Fig. S7B).
Since the inter-dyad distance distribution has a distinct oscillatory component, it is not surprising that a purely linear model of unwrapping energy does not fit the data as well, although it does match its overall shape ( Fig. S8A; SI Appendix, Model E). Less trivially, it was suggested on the basis of singlenucleosome unzipping experiments that nucleosome unwrapping proceeds with 5-bp periodicity because histones interact with each DNA strand separately where the DNA minor groove faces the histone octamer surface, creating two distinct contact "subpatches" [44]. This single-molecule data was fit to a model with a step-wise unwrapping free energy profile [45]. Each step in the profile corresponds to breaking a point histone-DNA contact, and the steps occur every 5.25 bp on average. We do not find any evidence for 5 bp periodicity of nucleosome unwrapping in the genomic data. Indeed, both 5 bp step-wise and 5 bp periodic sinusoidal profiles fit the data poorly, about as well as the linear model (Fig. S8B,C). Even the 10-bp step-wise unwrapping profile, while clearly having the right periodicity, does not fit the data as well as the structure-based model (Fig. S8D). This observation suggests that the picture of gradual loss of favorable finite-range histone-DNA interactions followed by gain in DNA conformational entropy is closer to reality than abrupt disruption of short-range histone-DNA contacts. A direct comparison of single-molecule and genomewide energy profiles is unfortunately obscured by the fact that the reported single-nucleosome unzipping experiments are specific to the 601 nucleosomeforming sequence [46], in contrast to our methodology which provides the average, sequence-independent picture of unwrapping energetics.  [50,51]. Interestingly, promoter nucleosomes have shorter inter-dyad distances and are therefore more unwrapped (Fig. 2B). When averaged over all genes, the number of dyads at a given bp and the average inter-dyad distance at that bp are strongly correlated (compare blue and red lines in Fig. 2C). The profile of average inter-dyad distances is also correlated with the distribution of DNA fragment lengths in an MNase assay which mapped both nucleosomes and subnucleosome-size particles by pairedend sequencing (Fig. S9A, green line in Fig. 2C) [47]. The two profiles do not coincide completely because inter-dyad distances also depend on the distribution of linker lengths. The observed behavior is opposite of the naive expectation that unwrapping increases with occupancy due to nucleosome crowding, but can be reproduced in a simple sequence-independent model in which nucleosomes phase off a potential barrier placed in the promoter region ( Fig. S10) [52]. Partially unwrapped nucleosomes tend to have elevated histone turnover rates [49], both around TSS and genome-wide (Figs. 2C,S9B,S9D). We find that nucleosomes at loci enriched in PICs [48] are also more unwrapped (Figs. 2C,S9C). Finally, inter-dyad distances tend to increase with the fraction of A/T nucleotides, indicating that nucleosomes occupying A/T-rich sequences have longer footprints genome-wide (Fig. S9E). We note that it is misleading to equate inter-nucleosome distances with peak-to-peak distances in the average profile of nucleosome dyad counts (blue line in Fig. 2C). The peak-to-peak distances are 164-165 bp, while the average inter-dyad distance for the nucleosomes in the [TSS,TSS+1000] region is 149.6 bp. Thus nucleosome unwrapping is much more common than could be predicted by mapping single-nucleosome positions alone. Distribution of the average distance between neighboring dyads. For each bp, the distances between a dyad at that bp and all neighboring dyads were averaged. ORFs are sorted as in (A). In (A) and (B), values at bps without dyads were obtained by interpolation, and heatmaps were smoothed using a 2D Gaussian kernel with σ = 3 pixels. (C) Data in (A), (B) and Fig. S9A-C averaged over all genes. Blue: nucleosome dyad counts, red: average distance between neighboring dyads, green: average length of DNA-bound particles mapped by MNase digestion [47] (see Fig. S9A for details). Curve with light gray background: combined occupancy of 9 PICs (TBP, TFIIA, TFIIB, TFIID, TFIIE, TFIIF, TFIIH, TFIIK, PolII) [48], curve with light pink background: average histone turnover rate [49]. The peaks in the dyad count profile (blue) are marked with orange ovals representing nucleosomes, and peak-to-peak distances are shown.   Figure 3: Probability of binding site exposure within a nucleosome. The solid blue and dashed green lines represent model predictions with and without unwrapping, respectively. In the latter case, a min = a max = 147 bp and all other parameters are adopted from the model with unwrapping. The dyad is fixed at bp 74. (A) Restriction enzyme sites inserted into the 5S rRNA sequence at locations indicated by the centers of vertical red bars [10]. (B) Restriction enzyme sites inserted into the 601 sequence at locations indicated by the centers of the vertical red bars in the middle of each group [11]. Each group of three bars corresponds to independent measurements in which the 601 sequence was flanked by different DNA sequences. In (A) and (B), the height of each bar is the equilibrium constant for site exposure averaged over multiple experiments (error bars show standard deviation).

Accessibility of nucleosomal DNA to factor binding
Partial unwrapping of nucleosomal DNA results in differential accessibility of factor binding sites with respect to their position inside the nucleosome: sites on the edges are more accessible than those closer to the dyad. In contrast, all-ornone nucleosome formation should not be sensitive to the binding site position -a nucleosome, once unfolded, liberates its entire site. Polach and Widom [10] studied differential accessibility of 6 restriction enzymes to their target sites. The sites were placed at various positions throughout the 5S rRNA nucleosomal sequence (Fig. 3A). A later study used the 601 sequence and an extended set of 11 restriction enzymes (Fig. 3B) [11]. These studies measured equilibrium constants for site exposure K conf eq , which are related to the probability for a site to be accessible for binding: p open = K conf eq /(1 + K conf eq ) ≈ K conf eq [29]. We use our crystal structure-based unwrapping model ( Fig. 1B; SI Appendix, Model A) to fit the data on site accessibility [10,11]. Here the system consists of a single nucleosome and asymmetric unwrapping is allowed. We assume that a site becomes accessible for the enzyme only after an additional number of bps, d, have been unwrapped from the histone octamer surface [29]. We also assume that once the dyad is unwrapped from either end, the entire nucleosome is unfolded. Then the probability for a binding site to be accessible is given by  [45]. Interestingly, the 601 DNA has to unwrap more extensively past the binding site to allow access to restriction enzymes. Overall, our model reproduces the observed differential accessibility of restriction enzyme binding sites with respect to the nucleosome dyad (Fig. 3). The only outliers are StyI and BfaI binding sites in the 601 series which were not used in the fit and which, if unwrapping proceeds from the ends, cannot be more open than the Pml I site located further away from the dyad.

Nucleosome-induced cooperativity
If multiple biding sites reside within a single nucleosome, binding of one factor makes the other sites more accessible, in a phenomenon known as nucleosomeinduced cooperativity [16,18,25]. The cooperativity disappears in the absence of nucleosomes and reduces in extent with the distance between consecutive sites [16]. Moreover, the cooperativity is not observed if the two sites are on the opposite sides of the nucleosome dyad [13].
We can use our model of nucleosome unwrapping (SI Appendix, Model A with a min = 1 bp, a max = 147 bp) to capture all these aspects of nucleosomeinduced cooperativity (Fig. 4). Specifically, for sites located more than 40 bp away from the dyad site accessibility is strongly enhanced if DNA unwrapping is allowed (Fig. 4A). Interestingly, cooperativity between two TFs bound on the same side of the dyad is observed both with and without unwrapping (Fig. 4B). However, without unwrapping it is impossible to show that binding on the opposite sides of the dyad is not cooperative, as observed in experiments [13] (Fig. 4C). Furthermore, the decrease of cooperativity with distance [16] cannot be reproduced (Fig. 4D). Thus modeling transient nucleosome unwrapping is necessary for understanding how TFs and other DNA-binding proteins gain access to their nucleosome-covered sites.

Sequence-dependent nucleosome positioning and unwrapping
Here we test the assumption that sequence specificity of nucleosome formation may be neglected when inferring nucleosome unwrapping potential from genomic data. In general, the total free energy u nuc (k, l) of a nucleosome occupying bps k, . . . , l is a sum of a sequence-independent term u SI nuc (k, l) describing contacts between histone side chains and the DNA phosphate backbone (Fig. 1B), and a sequence-dependent term u SD nuc (k, l) describing DNA bending into the nucleosomal superhelix [36]. Our approach (Eq. (4), SI Appendix) can be used to infer both contributions from high-throughput maps of nucleosome positions. Since the number of different unwrapped species may be as high as several thousand depending on the maximum extent of unwrapping and the assumption of unwrapping symmetry, we expect robust inference to require levels of read coverage that are not currently available, especially in higher eukaryotes. Furthermore, resolution of MNase-based nucleosome maps is insufficient for capturing the 10-11 bp periodicity of the unwrapping potential. In the absence of high-resolution, high-coverage experimental data, we have tested our ability to predict nucleosome unwrapping energetics using a realistic model system.
Specifically, we assumed that u SD nuc depends only on the number of monoand dinucleotides in the nucleosomal DNA [54] (SI Appendix). The length of nucleosomal DNA ranges from a min = 74 bp to a max = 147 bp; all other parameters of u SI nuc are as in SI Appendix, Model A. Using Eq. (1) with µ = −13 k B T, we computed the exact nucleosome distribution n nuc 1 (k, l) for the S.cerevisiae chromosome I. We sampled paired-end nucleosomal reads [k, l] from n nuc 1 (k, l) until a desired level of read coverage (mean number of reads starting at a bp) was reached. From this finite sample, we constructed a histogram of nucleosome lengths P (∆) and used it to optimize the parameters of u SI nuc (in each step of the fitting procedure, u SI nuc alone was used to predict n nuc 1 , which in turn gave P (∆)). Next, we used Eq. (4) to predict u nuc (i, j), subtracted u SI nuc from it (assuming that the dyad is at the mid-point of each particle), and fit u SD nuc to the rest. Finally, u SI nuc + u SD nuc were used to computeñ nuc 1 , which was then compared with the exact profile n nuc 1 . We find that we are able to reproduce the unwrapping potential even if nucleosome formation is sequence-specific (Fig. S11A). However, the overall slope of the potential is overpredicted because the histogram of particle lengths is affected by well-positioned nucleosomes which have negative sequence-dependent energies. The average of these energies biases the slope. Nucleosome occupancy predicted using u SI nuc + u SD nuc yields a reasonable correlation with the exact profile even at 1 read per bp, although at least 10 reads per bp are required to reproduce dyad positions (Fig. S11B). In all panels, free energy of a fully wrapped nucleosome is − log(10 9 ) k B T; histone chemical potential is log(10 −6 ) k B T; TF binding energy is − log(10 10 ) k B T to cognate sites, and − log(10 6 ) k B T to all other sites; TF chemical potential is log(10 −9 ) k B T unless varied. Asymmetric unwrapping is allowed; in the model without unwrapping, a min = a max = 147 bp and all other parameters are the same as in the model with unwrapping.

Discussion and Conclusions
Since DNA wrapped in nucleosomes is sterically occluded, nucleosome formation prevents access of DNA-binding factors such as TFs or polymerase subunits to their cognate sites. Transient DNA unwrapping off the histone octamer surface, extensively studied in model systems [16,10,17,11,18,19,20,12,13,21,22,23], can facilitate such binding events but its importance on the genomic scale has been unclear. MNase-based maps of nucleosome positions employing pairedend sequencing have identified numerous subnucleosome-size particles [47,55].
However, it is impossible to deduce from these experiments which particles correspond to unwrapped nucleosomes. In addition, there is a concern that DNA fragments may have been under-or overdigested by MNase. These challenges have been overcome in a recent experiment in which chemically modified histones are used to map nucleosome dyads with single-bp resolution [30]. Used in conjunction with paired-end sequencing, the experiment provides information about both dyad positions and the distances between neighboring dyads. The latter is especially important since it probes the two-particle distribution in the same cell. The histogram of inter-dyad distances shows that at least 38.7% of all S.cerevisiae nucleosomes are partially unwrapped. Gene promoters are enriched in such nucleosomes, which also have higher histone turnover rates. Over coding regions, the distribution of wrapped DNA lengths oscillates in phase with the nucleosome occupancy profile, so that more stable nucleosomes are also more wrapped.
We have developed a statistical mechanics description of nucleosome arrays which allows DNA unwrapping while rigorously treating steric exclusion and sequence specificity of nucleosome formation. We have shown that prominent 10-11 bp periodicity in the distribution of inter-dyad distances [30] can be explained using an unwrapping energy profile based on the pattern of histone-DNA contacts in the nucleosome crystal structures. We were able to rule out alternative scenarios of step-wise unwrapping and 5 bp periodicity of the unwrapping process. Furthermore, fitting the distribution of inter-dyad lengths required accounting for linker length discretization, commonly thought to be imposed by chromatin fiber formation [41].
Our model yields estimates of nucleosome unwrapping energies consistent with previous biophysical experiments, and accounts for single-nucleosome observations which show that nucleosome-covered binding sites closer to the edge of the nucleosome are more easily accessible to DNA-binding factors and that binding of the first factor enhances binding of subsequent factors on the same side of the dyad. Finally, our approach reproduces the periodic distribution of wrapped DNA lengths in the vicinity of TSS if potential barriers are placed in the promoters. The barriers may be created in vivo by PICs [48], other DNA-bound factors, and chromatin remodelers. The extent of nucleosome unwrapping in the yeast genome suggests that its treatment should be included in all future models of nucleosome positioning and chromatin energetics.

Statistical mechanics of DNA-bound particles
Here we outline the exact theory for T molecular species simultaneously interacting with one-dimensional DNA (details and extensions are provided in SI Appendix). The DNA-bound particles are subject to steric exclusion and may also be partially unwrapped. Let u s (k, l) (s ∈ {1, 2, . . . , T }) denote the binding energy of a particle of type s attached to bps k, . . . , l. One can show that the one-body distribution of particles of type t is given by where Z is the grand canonical partition function, t, k|z|s, l = e β[µs−us(k,l)] δ t,s , β = 1/k B T is the inverse temperature, µ s is the chemical potential of particles of type s, and δ t,s is the Kronecker symbol. If the DNA length is L, vectors |t, k are T L-dimensional with a 1 at position (t − 1)L + k and 0 everywhere else. Similarly, the nearest-neighbor two-body distribution function is given by where Θ is the Heaviside step function, and Z − (k) and Z + (k) are partial partition functions for the DNA segments [1,k) and (k,L], respectively. Z − (k), Z + (k), and Z can be computed recursively, as shown in SI Appendix. Using Eq.
(1), one can find particle occupancy: where a t min (a t max ) is the minimum (maximum) length of a particle of type t. The inverse problem of predicting DNA binding energies from one-particle distributions can also be solved: where Z, Z − , Z + are again found recursively using only one-particle distributions n t 1 (i, j) as input (SI Appendix). Eq. (4) enables us to disentangle intrinsic sequence preferences and unwrapping energetics from steric effects.

Model of inter-dyad distances
To predict the distribution of inter-dyad lengths, we compute the conditional probability of having a nucleosome with the dyad at bp c + d, given that the adjacent upstream nucleosome has the dyad at bp c: where the probability distributions of the nucleosome centers can be computed using Eqs. (1) and (2) for a single particle type: Here, 2∆ 1,2 +1 are lengths of the particles centered at bp c and c+d, respectively. To estimate P (c + d|c), we use c = 5 kbp and a box of length L = 10 kbp, so that the boundaries of the box are far away. We ignore nucleosome sequence specificity, convolve P (c + d|c) with a kernel to account for site-specific chemical cleavage bias (SI Appendix), and fit model parameters to reproduce the observed distribution of inter-dyad distances. We consider a problem of mutually interacting particles (one-dimensional rods) that can be reversibly adsorbed to a one-dimensional lattice of L sites [here, DNA base pairs (bps)]. In order to model partial unwrapping of nucleosomal DNA off the surface of histone octamers, we allow the particles to cover a variable number of bps between a min and a max . We assume that the particles cannot overlap while they are attached to the lattice. This is implemented using hard-core interactions between adjacent particles. There are also hard walls at the ends of the lattice so that particles are prevented from running off it. In addition, we allow generic two-body interactions between nearest-neighbor particles. The attachment of a particle to the DNA modifies the total energy of the system in a sequence-specific manner. Physically, the binding energy may have contributions from DNA bending, electrostatic interactions, hydrogen bond formation, van der Waals contacts, etc. Thus a particle which covers bps k, k +1, . . . , l has a total one-body binding energy u(k, l). Note that for pairs of coordinates (k, l) such that l − k + 1 > a max or l − k + 1 < a min , u(k, l) = ∞, because all particles must have the length between a min and a max bps. The theory presented below is valid for arbitrary binding energies u(k, l).

Supplementary Information
Let Φ(j, k) be the two-body interaction between a pair of nearest-neighbor particles which cover base pairs . . . , j −1, j and k, k +1, . . . (k > j) (Fig. S1). In the case of nucleosomes, such interactions may be used to account for the effects of higher-order chromatin structure [52]. Although we do not focus on two-body interactions in this work, they are included below for the sake of generality. We impose where V (d) is an arbitrary interaction potential which depends only on the linear distance d between two neighboring particles.
For a fixed number of particles attached to the DNA, N , the canonical partition function is where β = 1/k B T is the inverse temperature (k B is the Boltzmann constant). Note that with our definitions of one-body energies, two-body interactions and Figure S1: Schematic illustration of one-body and two-body potentials in a multi-nucleosome system. Nucleosomes may be partially unwrapped, resulting in variable DNA footprints.
hard-wall boundary conditions, only legitimate configurations of non-overlapping particles will contribute to Eq. (S1).
In order to simplify the notation, we introduce two L × L matrices: Here k|M |l represents the element of matrix M in row k and column l. |l is a column vector of dimension L with 1 at position l and 0 everywhere else, and k| is a row vector with 1 at position k and 0 otherwise.
Let |J be a vector of dimension L with 1 at every position. Eq. (S1) gives If the particles are allowed to attach and detach from the lattice, the system has a variable number of particles, and the grand-canonical partition function is where µ is the chemical potential, N max is the maximum number of particles that can fit on L bp, I is the identity matrix, and k|z|l = e β[µ−u(k,l)] ≡ ζ(k, l).
Note that all particle configurations with N > N max do not contribute to Z, allowing us to extend the upper limit from N max to ∞. From the partition function, we can compute s-particle distribution functions (see the chapter by Stell in [56]): , and in general .
Using these relations, we obtain the one-particle distribution and the two-particle distribution (S4) In particular, the nearest-neighbor two-particle distribution is given by Eqs. (S3) and (S4) allow an obvious interpretation. To find the probability of starting a particle covering positions from k to l [Eq. (S3)], we have to add statistical weights of all the configurations that contain a particle at that position, and divide the resulting sum by the partition function. Similarly, in order to find the probability of a pair of particles, one covering positions i to j and the other covering positions k to l [Eq. (S4)], we need to sum statistical weights of all the configurations containing that pair of particles, and divide by the partition function.
With the one-particle distribution, we can define the occupancy of a bp i as the probability of finding that bp in contact with any particle. In other words, we need to sum the probabilities of all configurations in which particles cover bp i: Note that 1 − Occ(i) is the probability that bp i is not covered by any particles.

Φ(t,j; s,k)
u t (i, j) u s (k, l) i j k l Figure S2: Schematic illustration of one-body and two-body potentials in a system with multiple particle types. The model allows all particles to be in multiple stages of unwrapping. In practice, we allow nucleosomes to be partially unwrapped but the transcription factors (TFs) always have fixed DNA footprints.

Multiple particle types
The above formalism can be easily extended to the case in which T types of particles can attach to the one-dimensional lattice. Let the binding energy of a particle of type t ∈ {1, . . . , T } that covers bps i to j on the lattice be u t (i, j). The interaction between a particle of type t ending at position k, and the next particle of type s starting at position l will be denoted by Φ(t, k; s, l) (Fig. S2). Each particle of type t, when attached to the DNA, is in contact with a number of bps between a t min and a t max . Thus u t (i, j) = 0 if i and j do not satisfy the constraints a t min ≤ j − i + 1 ≤ a t max . Also, Φ(t, i; s, j) = ∞ for j ≤ i since the particles cannot overlap. With this notation, the grand-canonical partition function becomes where µ t is the chemical potential of the particles of type t. The sum is over all configurations, which can have variable numbers of particles of any type.
Using the matrix notation, t, k|z|s, l = e −β[ut(k,l)−µt] δ t,s ≡ ζ t (k, l)δ t,s , t, k|w|s, l = e −βΦ(t,k;s,l) , where δ t,s is the Kronecker delta symbol, the partition function can be written as Z = all states Each vector |t, i has dimension T L and only one non-zero element, set to 1 for normalization. For example, |1, i vectors have a 1 at position i, |2, i vectors have a 1 at position L + i, etc. As above, we denote by |J a vector in which all T L elements are equal to 1, to obtain Similarly to the previous case of a single particle type, we can compute the one-particle density We can also obtain the two-particle density Eqs. (S6) and (S7) give the joint probability that a particle of type t covers bps i to j, while a second particle of type s covers bps k to l. Using Eq. (S5) we can compute occupancy for each type of particles t and for each bp i: n t 1 (k, l).
In the following sections we will focus on the one-particle density function.

Direct problem: Recursive solution for hard-core interactions
A straightforward application of Eqs. (S5) and (S6) entails computationally intensive matrix manipulations. Fortunately, for particles that interact only through hard-core repulsion rather than long-range two-body interactions, the one-particle distribution can be computed recursively and therefore much more efficiently.

General case
With multiple particle types, Eq. (S5) can be rewritten as where Z − (i) and Z + (j) are the partition functions for the domains [1, i) and (j, L], respectively. Note that in the case of hard-core interactions alone, Z − (i) and Z + (j) do not depend on the type of the particle occupying positions i through j.
In the case of steric exclusion alone, these partial partition functions satisfy the following recursion relations: and Here each particle type s has two characteristic lengths, corresponding to its minimum and maximum DNA footprints, respectively: a s min and a s max . The boundary conditions are Z − (1) = 1 and Z + (L) = 1. The full partition function is given by Z = Z − (L + 1) = Z + (0). Note that all unphysical terms for which bound particles run off the lattice automatically vanish from Eqs. (S8) and (S9). To avoid numeric instabilities, the recursion should be done in log space. Let with the boundary conditions F (1) = R(L) = 0. Then the one-particle distribution function is Although we do not show it here, the two-particle distribution can be computed in a very similar way. The only new ingredient in Eq. (S6) is the partition function for the box with walls at two arbitrary positions, Z(t, j, s, k) ≡ t, j|w(I − zw) −1 |s, k . This partition function can be computed recursively, exactly as the partial partition functions Z ± discussed above.

Special case: No unwrapping
The special case in which all particles are fully attached to their DNA sites (i.e., there is no DNA unwrapping) can be easily obtained from our general formalism. Indeed, in this case we restrict a s min = a s max = a s in Eq. (S10), obtaining As before, the boundary conditions are F (1) = R(L) = 0. The one-particle distribution is given by

Inverse problem: Recursive solution for hard-core interactions
In the previous section, we have solved the direct problem: given the binding energies for all particle types, we compute s-particle distributions. However, typically it is particle distributions that are observed experimentally, and the energetics of particle-DNA interactions need to be inferred. This inverse problem can be solved recursively for the case of systems with multiple particle types, partial unwrapping (variable footprints), and steric exclusion. The recursive solution is efficient enough to be employed on the genome-wide scale. Here we focus on one-particle distributions and one-body energies; the exact matrix formulation of the inverse problem for a single particle type with the two-particle distribution and the two-body potential is available in Ref. [52].

General case
Using Eqs. (S5), (S8) and (S9), we obtain: where N R (i) = t i−a t max +1≤j≤i−a t min +1 n t 1 (j, i) represents the probability of finding a particle of any type with the right edge at bp i, and ξ(i) = Z − (i)Z + (i)/Z. Z + satisfies a similar recursive relation: where N L (i) is the probability of finding a particle of any type with the left edge at bp i: where the initial condition ξ(0) = 1 has been used. After we compute both Z − and Z + in this way, the total partition function is given by Z = Z − (L + 1) = Z + (0) as before, and the binding energy for any particle attached to the DNA is given by . (S14)

Special case: No unwrapping
In the case of the all-or-none binding, all matrix elements i|n t 1 |j vanish unless j = i + a t − 1, where a t is the length of the binding site for the particle of type t. Thus Using these expressions, we can employ Eqs. (S11), (S12) and (S13) to compute Z + and Z − in log space. Finally, Eq. (S14) can be used to compute the binding energies.
If all particles are of the same type, the quantity ξ can be simplified further: where Occ(i) is the probability that bp i is covered by a particle. Thus in this limit ξ(i) is simply the probability that bp i is not bound by any particles. The recursion relations for Z − and Z + become .
These expressions are equivalent to those previously obtained in [54].

Sequence-specific nucleosome formation energies
The simplest sequence-dependent model of nucleosome formation assumes that the nucleosome energy depends only on the mono-and dinucleotide counts in the nucleosomal site [54]. Thus for the sequence S 1 S 2 . . . S N , the nucleosome formation energy is SiSi+1 , where Si and SiSi+1 are the mono-and dinucleotide contributions, respectively.
Here we consider a symmetrized version of the model, with A = T , C = G and 10 unique dinucleotide energies: AA/T T , AC/GT , AG/CT , AT , CA/T G , Using Eq. (S14), we can estimate all 12 's for a model in which nucleosomes are the only particle type and unwrapping is allowed. For each sequence starting a bp i and ending at bp j, we have where m X and m XY are the counts of mono-and dinucleotides X and XY in the sequence, respectively.
Using all possible combinations of pairs (i, j) where a nucleosome can form, we obtain a large number P of equations of this type: Here, E − µ is a column vector of dimension P , where each row contains one u(i, j) − µ element from Eq. (S15). µ is the column vector from Eq. (S15), and M is a P × 13 matrix with mono-and dinucleotide counts and -1's in the last column. From this equation we can derive the energies and µ by a least squares fit. A problem that arises is that matrix M has a one-dimensional kernel, or null space [52]. Because in every DNA sequence the number of mononucleotides is equal to the length of the sequence, and the number of dinucleotides is equal to the length of the sequence minus 1, the columns of the matrix M are not linearly independent. Indeed, the column vector is the only linearly independent vector from the kernel of M : M |V = 0 (i.e., the kernel of M is spanned by |V ). Thus we have only 12 linearly independent equations, and cannot obtain a unique solution for the 13 parameters and µ.
We need to add one constraint, which we choose to be to make the solution unique.

Site-specific chemical cleavage bias
Hydroxyl radicals that cleave DNA near the nucleosome dyad have two preferred cutting sites, at positions -1 bp and +6 bp with respect to the dyad [30]. If DNA is cut at these positions with frequencies f and 1 − f , the distance between two consecutive cuts (one on the Watson and one on the Crick strand) is given by where d dyads is the distance between two neighboring dyads and the bias b is The cleavage bias has to be taken into account by convolving the predicted inter-dyad distance probability P (c + d|c) with a kernel corresponding to this bias: The convolution is then compared with the observed distribution of interdyad distances.

Parameter optimization
Parameter fitting for all models was carried out in a two-stage procedure using the genetic algorithm optimization function ga from the MATLAB Global optimization toolbox. First, the objective function to be minimized was set equal to the root-mean-square (RMS) deviation between predicted and observed inter-dyad distributions. Once RM S < 10 −3 had been achieved, the objective function was replaced by RMS−r osc −r osc , where r osc is the linear correlation between observed and predicted oscillations after the smooth background has been subtracted from inter-dyad distributions, as in Fig. 1D. We have found that the two-stage setup allows us to effectively fit both the overall shape and the fine oscillatory structure in the data. The best-fit parameters for all models are given in SI Results.

Supplementary Results
Model A: Crystal structure augmented with an additional well. The binding energy of a particle of length a = 1 + x 1 + x 2 (1 bp for the dyad, and x 1 and x 2 for the extra number of bps in contact with the histone octamer on each side of the dyad) is given by u = u half (x 1 ) + u half (x 2 ), with where E b is the binding energy of a fully wrapped particle in the absence of 10-11 bp oscillations. The oscillations are based on the the positions of the histone-DNA contacts in the crystal structure [57]. The MATLAB function interp1(...) was used to generate the oscillatory pattern by piecewise cubic Hermite interpolation between the following data points:   Table S1: Fitted parameters for Model A. a max and a min are the maximum and minimum lengths of the nucleosome particle, µ is the histone octamer chemical potential, A is the amplitude of the oscillations, and f is the hydroxyl radical cutting frequency. p and d are the position and the depth of the first minimum outside of the nucleosome core particle, respectively. E b is the binding energy of a fully wrapped particle in the absence of 10-11 bp oscillations. Model B: Crystal structure augmented with a linear function. Same as Model A for x ∈ [1, 73], followed by a linear function: where ∆E/∆X is the slope of the linear function (i.e., ∆E is the energy difference between the first and last points of the linear function and ∆X is the cardinality of the range of the linear function).   All residuals are defined as in Model A.

Parameter
Model C: 10-bp oscillations superimposed onto a linear function.
where A is the amplitude of the oscillations, x 0 determines the phase of the oscillations, and E b is the binding energy of a fully wrapped particle in the absence of the oscillations.   All residuals are defined as in Model A.
Model D: 11-bp oscillations superimposed onto a linear function.
where A, x 0 and E b have the same meaning as in Model C.  All residuals are defined as in Model A.

Parameter
Model E: Uniform unwrapping.
where E b is the binding energy of a fully wrapped nucleosome.  All residuals are defined as in Model A.

Parameter
Model F: 5-bp oscillations superimposed onto a linear function.
where A, x 0 and E b have the same meaning as in Model C.   All residuals are defined as in Model A.
where E step is the amount of energy lost in each step, and x 0 determines the phase of the stepwise profile.  All residuals are defined as in Model A.
where E step is the amount of energy lost in each step, and x 0 determines the phase of the stepwise profile.    Figure S4: Sensitivity of the predicted inter-dyad distribution to model parameters describing higher-order chromatin structure. The unwrapping potential is based on nucleosome crystal structures (SI Results, Model A). (A) Root-mean-square error (RMS) of the predicted inter-dyad distribution, as a function of the position and the depth of the first minimum outside of the nucleosome core (Fig. 1B). The depth of the first minimum is computed with respect to u half (x = 73 bp). (B) The linear correlation coefficient between oscillations in the predicted and observed inter-dyad distributions (r osc ), as a function of the position and the depth of the first minimum outside of the nucleosome core (Fig. 1B). The oscillations were obtained by subtracting the smooth background from inter-dyad distributions, as described in the Fig. 1 caption. In both panels, model parameters that were not varied were kept fixed at their best-fit values (SI Results, Model A).  Nucleosome unwrapping/higher-order structure potential energy profiles. Under the symmetric unwrapping assumption, the energy of a nucleosome that covers 2x+1 bps is given by 2u half (x). The minima and maxima of the energy landscape are either based on the crystal structures of the nucleosome core particle as in Figure  In each column, the upper panel shows the nucleosome unwrapping/higher-order structure potential energy profile (as in Fig. 1B), the middle panel shows the comparison of experimental and predicted inter-dyad distance distributions (as in Fig. 1C), and the lower panel shows observed and predicted oscillations in the inter-dyad distance distributions (as in Fig. 1D).  Figure S9: Genome-wide distribution of nucleosome lengths, histone turnover rates, and transcription pre-initiation complexes. (A) Distribution of average lengths of DNA-bound particles mapped by MNase digestion [47] in the vicinity of TSS. We considered particles with sizes between 80 and 200 bp and assigned particle lengths to the mid-point of each particle. Values for bps without dyads were obtained by interpolation. (B) Distribution of histone turnover rates [49] in the vicinity of TSS. (C) Distribution of the combined occupancy of 9 transcription pre-initiation complexes (PICs) [48] in the vicinity of TSS. PIC occupancies provided at 20 bp interval in Ref. [48] were interpolated. In panels (A)-(C), gene order is as in Figure 1B, and the heatmaps were smoothed using a 2D Gaussian kernel with σ = 3 pixels. (D) Correlation between inter-dyad distances and histone turnover rates averaged over 10 kbp windows tiling the yeast genome. (E) Correlation between average inter-dyad distances and the A/T ratio in 10 kbp windows tiling the yeast genome. A/T ratio is the fraction of A/T nucleotides in the window, divided by the genomewide A/T fraction. Correlation in (D) and (E) refers to the linear correlation coefficient. We align all yeast genes by their TSS as in Fig. 2 and for each bp compute the fraction of times that bp is found in an intergenic region, as opposed to the ORF of a neighboring gene (grey background curve). The intergenic fraction has an asymmetric shape with a maximum at about 50 bp upstream of the TSS. We use this shape as a guide for constructing an energy barrier for in vivo histone deposition (pink background curve). (x > c). The free parameters of the barrier are fit to maximize the sum of two correlations: between observed and predicted normalized dyad counts [30] (solid and dashed blue lines, respectively), and between observed and predicted average nucleosome DNA lengths [47] (solid and dashed green lines, respectively). Normalized dyad counts are computed as the total number of dyads at a given bp for all genes, divided by the average of this quantity in a [-1000,1000] bp window around the TSS. Average DNA lengths are computed for all nucleosomes with a midpoint at a given bp, for all genes (if the midpoint falls in between two bps, the one on the left is used). The fitted parameters are: H = 0.0545 k B T, D = 0.0243 k B T, c = x TSS − 32 bp, σ 1 = 162.7 bp, σ 2 = 28.0 bp, σ 3 = 2090.9 bp, where x TSS is the absolute position of the TSS in the box, c is the center of the 3 Gaussians, H is the height of the first Gaussian, D is the depth of the third Gaussian, and σ 1 , σ 2 , σ 3 are the standard deviations of the three Gaussian distributions. The simulations were done in a 15 kbp box with the barrier placed at its center to eliminate the boundary effects. Unwrapping was assumed to be symmetric and the nucleosome structure-based unwrapping potential (SI Results, Model A) was used. The total free energy u nuc (k, l) of a nucleosome occupying bps k, . . . , l is a sum of u SI nuc and u barrier nuc = l j=k j , where j is the value of the barrier at bp j. Note that the chemical mapping data underestimates the number of -1 and +1 nucleosomes due to gel selection bias [30].  Figure S11: Inference of the unwrapping potential in a model with sequence-specific nucleosome formation energies. (A) Exact unwrapping potential vs. unwrapping potentials predicted at three levels of sequence coverage. All calculations were done on S. cerevisiae chromosome I with L = 230 kbp. M L reads were randomly sampled from the exact distribution n nuc 1 (i, j) (see Materials and Methods), where M ∈ {1, 10, 100} is the desired level of read coverage per bp. Sampled reads were used to compile a chromosome-wide histogram of nucleosome DNA lengths, to which the unwrapping potential in SI Results, Model A was fit using a genetic algorithm optimization function ga from MATLAB, which minimizes the root-mean-square error of the predicted distribution of nucleosome lengths. (B) The relative errors between predicted and exact (SI Results, Model A) parameters of the unwrapping potential at three levels of sequence coverage (light blue background). P denotes any parameter on the horizontal axis. Linear correlation between predicted and exact distributions of dyad positions and nucleosome occupancy (light pink background). The height of each bar represents the mean relative error for the corresponding parameter or the mean correlation coefficient, obtained by averaging the results of 100 random sampling experiments. The uncertainty intervals represent standard deviations.