{'fraction_non_alphanumeric': 0.03219076005961252, 'fraction_numerical': 0.017337307501241926, 'mean_word_length': 4.483791882320894, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The entropic lock and key of the histone code

The intricate pattern of chemical modifications on DNA and histones, the"histone code", is considered to be a key gene regulation factor. Multivalency is seen by many as an essential instrument to transmit the"encoded"information to the transcription machinery via multi-domain effector proteins and chromatin-associated complexes. However, as examples of multivalent histone engagement accumulate, an apparent contradiction is emerging. The isolated effector domains are notably weak binders, thus it is often asserted that the entropic cost of orienting multiple domains can be"prepaid"by a rigid tether. Meanwhile, evidence suggests that the tethers are largely disordered and offer little rigidity. Here we consider a mechanism to"prepay"the entropic costs of orienting the domains for binding, not through rigidity of the tether but through the careful spacing of the modifications on chromatin. An all-atom molecular dynamics study of the most fully characterized multivalent chromatin effector conforms to the conditions for an optimal free-energy payout, as predicted by the model discussed here.

In eukaryotic cells, DNA is wrapped around histone octamers to form nucleosomes and is then packaged into chromatin, a higher order structure. The tails of the histones extend out beyond the DNA and display numerous post-translational modifications (PTMs) [1]. Along with histones, nearly 30 additional proteins per nucleosome are associated with chromatin and there are tens of millions of nucleosomes per nucleus [2]. Some of these additional proteins -"effectors" -can bind to the histone PTMs and consequently regulate gene expression by impacting the shape and function of chromatin. The hypothesis that patterns of PTMs transmit biologically relevant signals through recruitment of specific binding partners is referred to as the histone code [3,4].
Multivalency was invoked within the histone code as a way to (i) overcome the promiscuous binding of different PTMs exhibited by isolated effector domains, (ii) to enhance binding relative to that of individual domains, and (iii) to maintain a mechanism for rapid response to a change in cellular signaling or environmental stress. [5,6] Since then, several studies demonstrated the possibility of multivalent chromatin-effector interactions in designed in vitro systems [7][8][9], but it was not until very recently that such interactions were revealed to exist in vivo and be responsible for cellular effects [10][11][12]. In these instances, multiple structured domains connected by flexible disordered tethers allow the effector to enjoy multivalent engagement with chromatin. [5,6,13,14] One example of a multi-domain effector is the E3 ubiquitin ligase UHRF1. The linked tandem tudor domain (TTD) and plant homeodomain (PHD) of UHRF1 are known to engage the H3 histone tail in a multivalent interaction. [15,16] The TTD-PHD cassette of UHRF1 is the most well-characterized example of multivalency in the histone code, complete with structural coordinates of the bound state [15] and an in vivo demonstration that the multivalent bound state is the biologically significant state [12].
While some initial effort was made to provide a qualitative depiction of the histone code thermodynamics, multiple questions still remain unanswered. For instance, multivalency was invoked in part to compensate for the weak binding of isolated effector domains. If multiple domains are connected by a flexible tether then some entropic penalty must be levied against the free energy of multivalent binding and this could be significant. To save the situation, one can suppose the tethers are rigid. [5,6] The problem is that the inter-domain tethers found in multidomain chromatin effector proteins are almost unanimously identified as disordered by a consensus of several reliable techniques [17].
In what follows we suggest an alternative to this conundrum of entropy compensation that does not rely on tether rigidity. Instead, we hypothesize that the spacings of PTMs on chromatin are "paired" to an effector so as to minimize the entropic penalty associated with orienting the effector protein domains -resulting in an entropic lock and key. Multivalent interactions have been studied extensively outside of the chromatin context [18][19][20][21][22] but the problem at hand is unique in that we explore the possibility that selectivity in the histone code is driven, in part, by the minimization of entropic losses upon chromatin binding via PTM spacing. In this view, effector proteins act as a "feeler gage" that selectively identify sets of PTMs spaced at particular intervals, a view which highlights the importance of the histone sequences which carry the PTMs as those sequences ultimately determine this spacing.
This alternative means of paying the entropic costs is discussed below through a simple model of the TTD-PHD cassette of UHRF1, wherein the TTD and PHD are non-interacting points connected through a third point representing the inter-domain tether. Given the weak binding enthalpies of the individual effector domains, we suggest that there are two possibilities for designing the chromatin-effector pair. One possibility is that the tether is designed to optimize the free energy of binding for a given chromatin site. The other possibility is that the binding site is such that the free energy of binding is optimal given an inter-domain tether -we refer to this situation as the entropic lock and key because, as we show below, the entropic penalty of orienting the two binding domains has been minimized. These two possibilities place different physical constraints on the geometry of the effector-chromatin pair and may imply different evolutionary pathways or functional contexts. For example, cis-and trans-histone binding may differentially exhibit these two possibilities.
Our model consists of two particles in a box with edge lengths l and volume l 3 and we limit ourselves to a special case of multivalency, a divalent binding event. The particles p 1 and p 2 represent two effector protein domains. We impose two binding sites in the box, located at points S 1 and S 2 , and we ascribe to each a radius of interaction R 1 and R 2 respectively. The binding sites are separated by a distance D. In this model, a binding event occurs if the effector p i is within R i of S i . The i-particle cannot bind the j-particle binding site.
In this system, the probability of finding p 1 bound to S 1 and p 2 bound to S 2 is straight forward, where X is the configuration space contained in the vol-ume l 3 , 1 [·] are indicator functions which are equal to unity when the argument is satisfied and zero otherwise. We use |p i − S i | to indicate the distance between the points p i and S i . We also use X dy to indicate integration over the three dimensional Cartesian space of particle y. P 0 evaluates to The probability of a "dual-binding" event is simply the product of the probabilities that p 1 is bound to S 1 and that p 2 is bound to S 2 . We take this as our reference for measuring the impacts of tethering the two particles together so that binding S 1 and S 2 becomes a divalent affair. Note that the free energy of binding is Suppose we introduce a third particle, p 12 , that serves as a "tether" and results in a three particle coarse graining of a dual-domain effector protein. We couple p 2 and p 1 to the tether with a potential energy function that is zero if |p i − p 12 | ≤ L/2 and infinity if |p i − p 12 | > L/2 for i = 1, 2. The two domains are thus always together, within a sphere of radius L/2 and we have an effective tether length of L. We refer to this sphere as the restraint sphere.
Formally, the configuration integral for this coupled system is dp 1 dp 2 dp 12 X X X 1 [|p1−p12|≤L/2] 1 [|p2−p12|≤L/2] dp 1 dp 2 dp 12 . (2) We make the following arguments under the approximations: i) collisions between the particles and the walls of the container may be ignored so long as the center of the restraint sphere, p 12 , is inside the container, i.e.,L << l, and ii) we consider only whole volumes of the binding sites even when the restraint sphere technically only covers a fraction of the volume.
Under these approximations, we can construct the probability of observing the dual bound state P (Equation (2)) as the product of three probabilities: The first factor is the probability that the restraint sphere containing p 1 and p 2 also contains S 1 and S 2 . The second factor is the probability that p 1 is in S 1 given that the restraint sphere contains the particles and the binding sites. The third factor is just like the second factor, except it is evaluated for p 2 in S 2 . These last two factors are expressed as the ratio of the volume of a given binding site to the volume of the restraint sphere.
The first factor of the probability can be deduced as follows. When p 1 is at S 1 , the possible positions of p 12 trace a sphere of radius L/2. When p 2 is at S 2 , p 12 traces out another sphere of radius L/2. If p 12 is found in the intersection of these two spheres, then divalent binding may be observed. The volume of this intersection can be computed as the sum of two spherical caps. The volume of a spherical cap is V c = πh 2 (3L/2 − h)/3 where h = (L − D)/2 is the height of the cap and L/2 is the radius of the sphere. By symmetry, the volume of the intersection is 2V c . Thus, the ratio of 2V c to the volume of the container is the probability that the restraint sphere contains both binding sites. the dual-bound state under the influence of this tether is An illustration of this model is given in figure 1. Notice that in this equation the absolute scale is irrelevant and only the ratios of the parameters are important. The approach of Diestler and Knapp [21] could be used to arrive at equation (3) but in this case it is easy enough to approximate the configuration integrals directly. Given the probabilities of the monovalent and divalent binding events, we can now estimate the "payout" of divalency We can relate the coupled and uncoupled systems with P = P 0 f and write the free energy of binding as In this case, by design, the enthalpies of binding events are additive but additivity will not hold in general [18] although it is a common approximation in the context of multivalency [5,6,20]. The second term in equation (5) replaces one entropic contribution of k B T ln(l 3 ) from ∆G 0 with k B T ln L 3 . Thus, the full entropy of the two independent domains is reduced, effectively leaving one domain to sample the full container while the other domain only samples the restraint sphere. When l 3 is large this contribution would be as important as the addition of binding enthalpies. This term "prepays" for the translational freedom of one of the domains that would be lost on binding and is responsible for transforming the probability of dual engagement from fleeting to expected. The average mammalian cell nucleus has a radius of about 3 × 10 4Å . If 80 to 99% of the nucleus volume is filled, the volume of the empty space is still on the order of 10 9 to 10 10 cubicÅ suggesting a payout of billions to trillions.
The third term in equation (5) compares the number of ways the restraint sphere can cover the binding sites (the volume of the spherical caps) to the number of ways the two domains can be arranged (the volume of the restraint sphere). This term is the entropic penalty for simultaneously organizing the two domains into their binding sites. As the number of ways to cover the two sites shrinks (as D approaches L) this term can completely abolish binding, reflecting a poor fit between the dual-domains and the binding site geometry.
Equation (3) can be maximized in two ways: i) the tether length L can be optimized for fixed D or ii) the binding site distance D can be optimized for fixed L. The second case is trivial, when D = 0 the entropic penalty vanishes. In practice D will be finite due to the excluded volumes of the domains. The first case can be treated by finding the maxima of (L + D 2 )(L − D) 2 /L 6 for fixed D. The optimal tether length is L = D/.618, which causes the domain-domain root-mean-squared separation (RMS) to coincide with D.
These two situations should be easily distinguished by comparing the domain-domain distance distribution with the binding site separation. If the binding site separation is close to the domain-domain RMS, then the tether has been optimized for the binding site. If on the other hand the binding site separation is close to the minimum domain-domain spacing, then the binding site geometry is optimal for binding but the tether has not necessarily been optimized for binding.
We were curious how effectors of the histone code might be paired with their chromatin targets, so we compared the dynamics of the multivalent effector UHRF1 and its target histone H3 with the model captured by equation (2). We have combined our previous molecular dynamics (MD) simulations of the TTD-PHD cassette of UHRF1 (nearly 4µs) [12] with a new simulation of the H3 tail and a Monte Carlo (MC) simulation of a 5-particle system obeying equation (2). In this example, p 1 represents the methyllysine binding pocket of the TTD and p 2 represents the N-terminal binding surface of the PHD. The H3 N-terminal nitrogen and the N ζ if lysine 9 represent the the binding sites S 1 and S 2 and were used to build the distance distribution p(D). We used the OH atom of Tyr188 and the backbone O atom of Pro355 to build the distance distribution p(|p 1 − p 2 |) for the TTD-PHD. These choices were based on proximity to the N ζ and N-terminus, respectively. Details of the TTD-PHD simulation are reported elsewhere [12] and details of the histone simulation are given in Supplemental Materials.
We computed the free energy along the coordinate D = |N − N ζ | with an adaptive biasing potential [23]. Error in the free energy along the N-N ζ distance was estimated as |A(ξ, t 1 ) − A(ξ, t 2 )|dξ for A(ξ, t) < 10kcal/mol, ξ = D, and t 1 = 96ns, t 2 = 128ns. The error in free energy was 1.8 kcal/mol and was confined to values of D < 20Å . The free energy along D leads to the distribution in D since p(D) ∝ exp[−βA(D)].
The distance distribution p(|p 1 −p 2 |) for the TTD-PHD (shown in blue in Figure 2) suggests a tether length of L UHRF1 = 65Å and a particle radius of 6Å. We performed Monte Carlo simulations of equation (2) that includes excluded volume effects where the ratio 6/65 is used to determine the length scale of the model. We used an excluded volume (or particle) radius of L × 6/65. The MC results shown in figure 2 (orange) were obtained with R 1 = R 2 , R/L = 0.117, L/l = 0.34 and various values of D. Optimal binding for the model system was observed when D/L = 0.313, which suggests an optimal PTM spacing of D = |N − N ζ | = 0.313 × 65 = 20.4Å for the TTD-PHD and H3 system. The distribution of D for H3 is shown in black in figure 2 and peaks near 22Å. This is in agreement with the peak in divalent binding probability which occurs at D = 20Å. This suggests that UHRF1 is paired with the N-terminus and methylated N ζ of lysine 9 of H3 because this spacing minimizes the entropic penalty associated with the divalent bound state. Independently of this spacing, binding can be impaired by chemical modification of other histone residues (methylation of Arg2 for example). Thus the entropic lock and key is only one component of the elaborate recognition mechanism that drives the histone code.
It was not possible to infer the parameters R and l from the existing TTD-PHD data. However, l only impacts the magnitude of the probability of observing the dualbound state and does not shift the optimal D : L ratio. We considered several R and found very little impact on the results for R/L = 0.06, 0.117, 0.17.
In 1894, Emil Fischer suggested an interaction model explaining the enzyme-substrate specified by the complementarity of their respective geometric shapes. [24] This type of interaction, generally referred to as "the lock and key", was mostly thought to be enthalpy-driven due to multiple contacts resulting from the perfect shape fit. Based on the results presented here, we suggest that an entropic "lock and key" might be a fundamental component of design that pairs effector proteins with their chromatin targets. Because the individual modules of multi-domain chromatin effectors are often weak binders to their respective sites and because they are connected by long disordered linkers, multivalent binding of such effectors to chromatin is likely controlled by entropy. The physical traits of multivalency resulting from the presented model may have genome-wide implications. For instance, the characteristics of tethers connecting the effector domains, such as length and plasticity, as well as the geometry and dynamics of the target PTMs are an as yet underappreciated component of the histone code. In the case of UHRF1, phosphorylation of the tether in the TTD-PHD cassette is known to weaken histone binding by about 30 fold [15] and correlates with enhanced production of TopoIIα [25]. Additionally, a mutation in the TTD-PHD linker (R295A/R296A) results in an increased tether length and drastically impacts histone binding [15] while causing a loss of function [12]. Hence, the function of UHRF1 may be altered by a single modification that changes the properties of the tether between TTD and PHD. Until now, the studies of the histone code effectors were mostly focused on the discovery of new effector domains followed by the elucidation of their structure and function. In light of the presented results, one might suggest that the inter-domain linkers could have been used as an evolutionary speed-dial; By mechanical alteration of a linker, nature could steer an effector away from the areas where its presence might be unwanted, which is less time-consuming and less risky than reengineering the effector domains.
The present work also provides a means of categorizing effector-chromatin pairs. The TTD-PHD of UHRF1, which binds cis-histone, appears to bind chromatin opportunistically by selecting chromatin sites with optimal D. Will other chromatin effectors demonstrate a similar relationship with their binding partners or will some of them display optimized tethers? These would be two fundamentally different classes of effectors and would represent a significant difference in effector design. It will be interesting to consider CHD4/5 [8,10] or BPTF [7], which are examples of trans-histone binders. Will the cis-/trans-histone preference correlate with the two modes of optimallity outlined here?
The research described here was supported by the US National Institute of General Medical Sciences; US National Institutes of Health (grant RC1GM090732 and R01GM100919); the Carolina Partnership and the University Cancer Research Fund; University of North Carolina at Chapel Hill.

Supplemental Materials
We used the biasing scheme from reference 23 with b = 0.8 and c = 0.1/∆t (∆t is the timestep) to compute the potential of mean force along the distance from the N-terminal nitrogen to N ζ of lysine 9 in an isolated H3 histone tail (residues 1-10) using a 128 ns trajectory. Simulation setup is the same as in reference 12 but here we use the GROMOS96ff force field.
To enhance sampling of backbone dihedrals of the histone we applied a bias potential along the dihedrals. This bias was derived from a converged free energy computation in the φ, ψ, and ω backbone dihedral angles of a tri-alanine peptide. The (φ, ψ, ω) computation was carried out with a recently introduced decomposition method [26]. In the end, this provides an angle-specific biasing potential that can enhance rotations along the histone backbone. The converged histograms in φ, ψ and ω from the tri-alanine simulation were used against the histone with b = 0.3, c = 0.1/∆t. The effects of this bias on the distribution of dihedral angles was removed through on-the-fly-reweighting [27,28] during the N-N ζ free energy computation.