Spatiotemporal organization of chromatin domains: role of interaction energy and polymer entropy
aa r X i v : . [ c ond - m a t . s o f t ] F e b Spatiotemporal organization of chromatin domains: role of interaction energy andpolymer entropy
Kiran Kumari,
1, 2, 3
J. Ravi Prakash, and Ranjith Padinhateeri IITB-Monash Research Academy, Indian Institute of Technology Bombay, Mumbai, Maharashtra - 400076, India Department of Biosciences and Bioengineering, Indian Institute of Technology Bombay, Mumbai 400076, India Department of Chemical Engineering, Monash University, Melbourne, VIC 3800, Australia
Chromatin is known to be organised into multiple domains of varying sizes and compaction. Whilethese domains are often imagined as static structures, they are highly dynamic and show cell-to-cellvariability. Since processes such as gene regulation and DNA replication occur in the context of thesedomains, it is important to understand their organization, fluctuation and dynamics. To simulatechromatin domains, one requires knowledge of interaction strengths among chromatin segments.Here, we derive interaction strength parameters from experimentally known contact maps, and useit to predict chromatin organization and dynamics. Taking α -globin domain as an example, weinvestigate its 3D organization, size/shape fluctuations, and dynamics of different segments withina domain, accounting for hydrodynamic effects. Perturbing the interaction strengths systematically,we quantify how epigenetic changes can alter the spatio-temporal nature of the domains. Computingdistance-distributions and relaxation times for different chromatin states, we show that weak andstrong interactions cooperatively determine the organization of the domains and how the solid-likeand liquid-like nature of chromatin gets altered as we vary epigenetic states. Quantifying dynamicsof chromatin segments within a domain, we show how the competition between polymer entropyand interaction energy influence the timescales of loop formation and maintenance of stable loops. PACS numbers:
Chromatin is a long polymer made of DNA and pro-teins; its four dimensional (spatial and temporal) orga-nization is crucial in regulating cellular processes liketranscription, replication and DNA repair [1–4]. Howthe chromatin self-organises spatially and temporally isa question of active research. Recent developments inchromosome conformation capture (Hi-C) [4–9] and mi-croscopy experiments [10–12] help us investigate contactsamong different segments and overall 3D organization ofthe chromatin polymer. Such experiments so far sug-gest that intra-chromatin interactions lead to the forma-tion of different types of domains within a single chromo-some [6, 7, 11–15].It has been hypothesized that loop extrusion and/orphase separation could be the mechanism for the forma-tion of these domains [14, 16, 17]. In the loop extrusionpicture, chromatin regions are actively brought togetherwith the help of proteins like cohesins/codensins and heldtogether by CCCTC-binding factor (CTCF) [13, 14, 18,19]. However, CTCF-dependent loops are found only ina fraction of the domains [13, 20]. Hence an alternativeproposal is that chromatin domains may also be passivelyformed via phase separation [17, 21, 22]. Recently, it hasbeen shown that in the absence of loop extruding fac-tors, chromatin does still form domains and execute thenecessary biological function[11, 17, 20, 23, 24] indicat-ing that micro phase-separation might be an importantmechanism. Phase separation is also known to bring to-gether certain enhancers and promoters segregating themfrom other regions [21, 22, 25]. In certain cases, as faras biological function is concerned, there is an ongoingdebate whether the actual contact is crucial or proximity — closeness in 3D without being physically in contact —would suffice [3, 17, 23, 26]. To probe this further, it isimportant to go beyond contacts and examine the wholeconfigurational space of chromatin polymer.Widely used experimental methods so far provide usonly static snapshots of chromatin contacts. Large Hi-Ccontact values are often assumed to represent static con-tacts holding together different segments of chromatin.However, note that experimentally inferred contact prob-ability values, for nearly all chromatin segments, are verysmall ( <<
1) [13]. This implies that the contacts will beoften broken and the chromatin polymer conformationcan be highly dynamic. Even in a steady-state, there arelikely to be large fluctuations, cell to cell and temporalvariabilities. While there have been several attempts tounderstand the static 3D configurations and their aver-age properties, the dynamic nature of chromatin domainsremains unclear.Complementing recent experimental efforts, there havebeen many computational studies investigating chro-matin organization [18, 19, 27–42]. However, to simulatethe formation and dynamics of domains, one requires op-timal interaction strength parameters such that experi-mentally known properties are recovered. We have re-cently developed an Inverse Brownian Dynamics (IBD)method to extract potential energy parameters that areconsistent with HiC-like experiments [29]. Using theintra-chromatin interactions computed from this method,in this work we go beyond the static picture and studychromatin fluctuations and dynamics perturbing epige-netic states. We use Brownian dynamics (BD) simula-tions with hydrodynamic interactions for this study. Go-ing beyond contact frequencies, we compute the full dis-tance distribution between all pairs of segments and ex-amine the cooperative nature of folding. We then studythe dynamics of the domain, compute relaxation times,loop formation times and contact times. We investigatehow the interaction strengths and polymer entropy influ-ence these measurable quantities. Finally, we discuss thesignificance of these findings.
Model and Methods
We consider chromatin as a bead spring chain hav-ing optimal intra-chromatin interactions derived from 5Cdata using an Inverse Brownian Dynamics (IBD) algo-rithm [29]. The total energy of the chromatin beadspring chain, made of N beads, is U = U S + U SDK where U S = P i H ( | r i − r i +1 | − r ) is the springpotential between the adjacent beads i and ( i + 1), r i is the position vector of bead i , r is the naturallength and H is the stiffness of the spring [43]. Tomimic protein-mediated interactions between bead-pairs,we use Soddemann-Duenweg-Kremer ( U SDK ) potentialwhose repulsive part is modelled by Weeks-Chandler-Anderson potential 4 (cid:20)(cid:16) σr ij (cid:17) − (cid:16) σr ij (cid:17) + (cid:21) − ǫ ij uptothe range 2 / σ . Here r ij = | r i − r j | is the distanceand ǫ ij is an independent parameter representing theattractive interaction strength between beads i and j .The repulsive part is unaffected by the choice of poten-tial parameter ǫ ij . The attractive part is modelled by ǫ ij (cid:2) cos ( αr ij + β ) − (cid:3) which, unlike the Lennard-Jonespotential, smoothly reaches zero at the cut off radius r c = 1 . σ [44, 45]. We simulate the chromatin polymerusing Brownian dynamics simulations, where the timeevolution of the bead positions are governed by an Itostochastic differential equation. This simulation accountsfor hydrodynamics interaction which is computed by theregularized Rotne-Prager-Yamakawa tensor. See supple-mentary information (SI) for details. For the simulation,all the length and time scales are non-dimensionalisedwith l H = p k B T /H and λ H = ζ/ H , respectively where T is the absolute temperature, k B is the Boltzmann con-stant, and ζ = 6 πη s a is the Stokes friction coefficient of aspherical bead of radius a with η s being the solvent vis-cosity. All the non-dimensional quantities are indicatedwith the asterik (*) symbol such as non-dimensional po-sition r ∗ i = r i /l H .For a chromatin, we do not know the interaction energystrengths `a priori . Given 5C/HiC data, we performed aninverse calculation and obtained the optimal interactionstrengths that are consistent with the experimental con-tact probability map [29]. The inverse algorithm worksas follows: we start with the initial guess values of in-teraction strengths, simulate the polymer following theconventional forward Brownian dynamics method and obtain the simulated contact probabilities in the steadystate. Based on a statistical method, the interactionstrengths are revised for the next iteration, dependingupon the difference between the simulated and the knownexperimental contact probabilities (see SI). We performseveral iterations of the loop (i.e., BD simulation, calcu-lation of contact probability and revision of interactionstrength) until the error between the simulated and ex-perimental contact probabilities is less than a predeter-mined tolerance value. Results and discussion
As a first step, to establish a connection between oursimulations and experiments, we computed the average3D spatial distance ( r ∗ ij = | r ∗ i − r ∗ j | ) between differentpairs of chromatin segments as a function of the corre-sponding 1D genomic separation ( s ij = | i − j | ). As shownin Fig. 1(a), h r ∗ ij i ∼ s νij with ν = 0 .
38 in the OFF state,suggesting near close packing within the chromatin do-main (red symbols). As a “control”, we also simulateda self-avoiding walk (SAW) polymer with no crosslinking( ǫ ij = 0) which results in ν = 0 . y -intercept ofthe experimental data gives us the size ( σ = l H ) of the10kb chromatin (a single bead in our simulation). Forthis experimental system (mouse chromosome 6, 1 . l H = 22nm. Even thoughwe do not have such extensive spatial distance data for α globin, we compared the available FISH data for α globinand deduced the l H = 36nm (see SI). Throughout thispaper, l H = 36nm and λ H = 0 .
1s are used to convert allnon-dimensional lengths and times into standard unitsand we will present quantities in both units. The reasonsfor the choice of both these specific values are discussedin greater detail in the SI.
Distance distributions and cooperative nature ofchromatin folding:
Even though the average dis-tance between two chromatin segments is often usedto represent chromatin organization, this may not de-scribe the accurate biological picture in a dynamic, het-erogeneous context. We compute the distribution of3D distance between different segments, p ( r ∗ ), as itcaptures the maximum information about variabilityand fluctuations of chromatin. As a control, we com-puted the p ( r ∗ ) for a SAW and it agrees well with theknown analytical expression of des Cloizeaux, p ( r ∗ ) = C ( r ∗ ) θ +2 e − ( Kr ∗ ) / (1 − ν ) [46] . Here ν is the Flory expo-nent, θ is a geometrical exponent and the coefficients C and K are functions of θ [46]. Recent work has ledto an accurate estimation for these constant which arediscussed in SI. Fig. 1(b) show the validation with inter-mediate beads of a SAW (with ν = 0 . , θ = 0 . , K =1 . , C = 2 . (a) (b) (c) (d)FIG. 1: Changes in distance probability distribution( p ( r ∗ ij ) ), cooperative nature of chromatin folding : (a)Average 3D distances between all bead-pairs as a function ofcorresponding genomic distances for the control simulations(SAW, black symbols), chromatin domain that we simulated(red symbols) and comparison with experimental data from[12] (blue symbols). The major axes (lower x and left y ) repre-sent quantities in dimensionless units (see methods) while theother axes (upper x and right y ) represent the same in stan-dard units. (b) Distance probability distribution p ( r ∗ ij ) fromsimulations compared with the analytical expression [46] fora pair of beads 9 and 34 of a SAW polymer. (c) p ( r ∗ ij ) forvarious bead pairs with different ǫ ij , but same s ij = 25 in theOFF state of α -globin gene. The interaction-driven peak ishighlighted in the inset. Vertical dashed lines represent h r ∗ ij i corresponding to each distribution. (d) Comparison of p ( r ∗ ij )between different “epigenetic states”. OFF: state with all WTinteractions in GM12878 (red), OFF GT1: when weak inter-actions are ignored; only with ǫ ij > k B T interactions (pink),OFF GT2: when only very high interactions ( ǫ ij > k B T ) areaccounted (blue), SAW: control simulation with no crosslink-ing (black). Vertical lines and inset are similar to (c). We then studied the p ( r ∗ ) for the α -globin gene lo-cus in GM12878 cell type. Examining various segments 250kb apart along the chain backbone (25 beads), we findthat all the distributions have a broad peak near their re-spective average distances (Fig. 1(c)). However, for beadpairs having high ǫ values, a sharp peak emerges near r ∗ ≈ r c — we call this an “attraction-driven peak” as itis within the attractive range of the potential. The heightof the peak is correlated with the strength of attrac-tion ( ǫ ). However, the average distances (vertical linesin Fig. 1(c)) appear independent of ǫ ij . This differenceis also reflected in the cumulative distribution functionas shown in SI. Together, these results imply that aver-age distances between bead-pairs may not represent thecomplete picture of chromatin organisation; understand-ing the whole distribution is necessary.Given that we have the optimal interaction strengthsthat satisfy the experimentally known contact probabil-ity constraints [29], we can answer the following impor-tant question: Are the measurable properties of a givenbead-pair (e.g. r , ) solely determined by the interac-tion between those two particular beads ( ǫ , ) or arethey influenced by the interactions among other bead-pairs as well? To answer this, we adopted the followingstrategy: we systematically switched off the attractive in-teraction among certain bead pairs and computed prob-ability distributions and other polymer properties. Wesimulated polymers for the following four cases: (i) allinteractions are considered – GM12878 (OFF), (ii) onlythose interactions above 1 k B T are considered – we call itOFF:GT1 – all weak interactions ( < k B T ) are switchedoff here, (iii) only strong interactions above 2 k B T areconsidered – OFF:GT2 – all weak and medium interac-tions ( < k B T ) are switched off (iv) all attractive inter-actions are switched off – the SAW polymer. These fourcases can be cosidered as four different epigenetic states –states having different interaction strengths due to under-lying epigenetic variations. Fig. 1(d) shows p ( r , ) for allthe four cases. When we switch off the weak interactionsbelow 1 k B T (OFF:GT1), compared to the OFF state, theheight of the interaction-driven peak of the distributiondecreases and overall the polymer swells resulting in theshift of the second peak (compare pink and red curvesin Fig. 1(d)). This implies that weak interactions hav-ing strengths comparable to thermal fluctuations can alsoinfluence the contact probability and polymer configura-tions. If we keep only the highly prominent interactionsand neglect all interactions below 2 k B T (OFF:GT2),the interaction-driven peak further diminishes and thedistribution function approaches the SAW distribution(compare blue with other curves in Fig. 1(d)). Notethat the interaction between beads 5 and 30 is present( ǫ , = 2 .
09) in all the cases except in the SAW case.These results suggest that the measurable properties fora given bead-pair (e.g. r , ) depends not only on theattraction strength of that particular bead pair but alsoon the interactions of the whole polymer chain. This re-sult implies that all bead-pairs collectively/cooperativilycontribute in determining the relative position for a par-ticular bead-pair. Epigenetic changes alter the volume and shape ofthe chromatin domains:
The above picture suggeststhat the chromatin folding is influenced by collective be-havior of all beads having different interaction strengths.To examine the nature of collective behavior, we probeda property of the whole polymer namely the radius of gy-ration ( R g ). To understand how folding is affected by dif-ferent epigenetic states, we did the following. We startedwith a polymer having no interactions (SAW), addedweak interactions (small ǫ ) that exist between beads inthe OFF state as the first step, equilibrated, computed R g and sequentially added stronger interactions betweenbeads step by step ( ǫ < . , ǫ < . , ..., ǫ < R g was computed (see Fig:2, top panel). From R g , wehave also computed the volume V = (4 / πR g of thechromatin domain as shown in the right side y -axis. Asseen from the figure, adding very weak interactions doesnot change the R g much. However, adding intermediateinteractions significantly reduces the R g and it saturatesas the interactions gets stronger, resulting in a sigmoidal-like curve showing signatures of cooperative/collectivebehavior. Since, we have equilibrated the polymer foreach set of ǫ values, the LHS of the curve can be in-terpreted in two ways: folding the polymer by addingstronger and stronger interaction starting with a com-pletely unfolded state or equivalently unfolding the poly-mer by removing the stronger interaction starting with acompletely folded OFF state. One can also ask how thepolymer would fold if one adds strong interactions (larger ǫ ) as the first step, starting with SAW, and then addweaker interaction sequentially step by step (denoted asGT1, GT2 etc). This is shown in the RHS of Fig. 2 (bluesymbols). The whole curve suggests that having promi-nent interactions alone or weaker interactions alone maynot take the system closer to its full equilibrium state.We also show typical snapshots of 3D chromatin config-urations corresponding to different epigenetic states. Asexpected, the OFF state is compact and the volume ofthe domain increases as we go towards the SAW state.The fold-change in volume we observe between the twoextreme states are roughly the same order as the densitychange observed experimentally [47].To quantify how the shape of the chromatin domainchanges with epigenetic states, we computed the as-phericity (B) and the acylindricity (C) parameters (see SIfor definition). Asphericity quantifies the extent of devi-ation from a spherical shape. If a polymer is coiled withthe average shape of a sphere, B = 0. Here a positive Bvalue suggests that even in the OFF state, the chromatindomain is not a perfect sphere. As we go from OFF toSAW, the asphericity increases by ≈
65% as shown in SI.
SAW LT1 LT2 OFF GT1 GT2 SAW202530354045
SAW LT1 LT2 OFF GT1 GT2 SAW4567
SAW LT1 LT2 OFF GT1 GT2 SAW678
FIG. 2:
Interaction strengths determining shapeproperties of chromatin domain:
Upper panel: Sizeof the chromatin domain ( R g ) as we perturb interactionstrengths. x -axis represents different interaction states withextreme ends representing the control (SAW) polymer and theOFF state (WT) in the middle. LT1 (LT x ) indicates that allinteractions below 1 k B T ( xk B T ) are present in the polymer.Similarly GT1 (GT x ) indicates that all interactions above1 k B T ( xk B T ) are present. The right y -axis indicates volumeof the chromatin domain in femtolitre. Snapshots from sim-ulations at various epigenetic states are shown around theperimeter of the graph. Bottom two panels represent thenormalized asphericity ( B/R g ) and acylindricity ( C/R g ), re-spectively. The x -axis is the same in all panels. However, the asphericity scaled with the polymer size(
B/R g ) changes by ≈ R g , we have shownthe GT (RHS, orange symbols) and LT (LHS, blue sym-bols) cases for the asphericity too. Even though bothsides are monotonically increasing, note that LT casesare not equivalent to the GT cases. We also compute theacylindricity parameter that quantifies the extent of thedeviation from a perfect cylinder. Here too, C >
C/R g ) is nearly a constant as shown in thelower panel of Fig. 2. Estimation of solid-like (stiffness) and liquid-like(drag) properties from domain relaxation timesand fluctuations:
Whether chromatin is liquid-like,solid-like or gel-like has been a matter of intense discus-sion in the recent literature [10, 24, 48, 49]. In the phaseseparation picture, chromatin segments are thought tobe “liquid-like”, dynamically exploring various configu-rations. Given that our model can study the stochasticnature of formation and breakage of bonds, and polymerdynamics, consistent with what is observed in Hi-C ex-periments, below we compute relaxation times and fluc-tuations of the chromatin domain and estimate effectiveelastic and drag properties. -1 OFF GT1 GT2 SAW5001000150020002500 (a) (b) -2 -1 OFF GT1 GT2 SAW02468101214 (c) (d)FIG. 3:
Stiffness and relaxation time affecting solid-like and liquid-like behavior: (a) Exponential decay ofend-to-end auto-correlation function with time for four epige-netic states computed with HI. (b) Relaxation times ( τ ) withand without HI reveal that HI helps in relaxing the poly-mer faster. (c) Effective stiffness between all the bead-pairs;OFF state is more stiff compared to less interacting states andSAW. Inset: Free energy as a function of bead pair distance r , . (d) Effective viscous drag felt by different chromatinstates. OFF state chromatin (with stronger interactions) ismore stiff and has a higher viscous drag. First, we computed the end-to-end vector autocorre-lation function h R ∗ E (0) · R ∗ E ( t ∗ ) i / (cid:10) R ∗ (0) (cid:11) where R ∗ E = | r ∗ − r ∗ | and extracted the longest relaxation time τ ∗ with and without HI. The autocorrelation decay com-puted with HI is shown in Fig. 3(a) and no-HI case isshown in SI. Fig. 3(b) shows that the relaxation times forall the epigenetic states are lower with HI, as observedpreviously for the protein folding simulations [50, 51]. Allresults presented in this paper are computed with HI, un-less stated otherwise. The chromatin in the OFF statehas a lower relaxation time compared to a crosslinking-free chromatin (SAW). This can be counter-intuitive asone would naively expect that the more the crosslinkingthe slower the chromatin will relax. A similar puzzle hasalso been observed in recent experiments [10, 52] wherethe repressed chromatin domain diffuses faster than the active one. To understand this apparent contradiction,we investigate the elastic and drag properties of the chro-matin domain. From the measurement of fluctuations ofeach bead-pair we can compute an effective stiffness de-fined as K ij = k B T / h| r i − r j | i (Fig. 3(c)). As expected,highly cross-linked OFF state is more stiff than the otherepigenetic states including SAW. This can also be under-stood from the free energy as a function of bead pairdistance F ∗ ( r ∗ ) = − ln ( p ( r ∗ ) / πr ∗ ) (see inset). Theabove behaviour is consistent with K ∗ ij ∼ ∂ F ∗ ij ∂r ij and stiff-ness ( K ∗ ij ) of different epigenetic states do show similarbehaviour. For the known stiffness and relaxation times,we can compute an effective drag coefficient defined as ζ ∗ eff = τ ∗ × K ∗ . Taking the effective stiffness of the endbeads ( K ∗ , ), we find that the drag for the OFF state ishigher than the other states suggesting that higher cross-linking reduces its ability to reorganize. Both the stiffnessand drag are greater for the OFF state than the SAW,but they combine to lead to a faster relaxation time forthe OFF state. Our findings agree with the recent ex-perimental report that crosslinked chromatin shows lessFRAP revealing gel-like nature of chromatin [24, 49]. Interplay between interaction energy and poly-mer entropy influences the dynamics of chromatindomain:
While we gained insights into steady state fluc-tuations and distance distributions, how the interactionswould affect chromatin dynamics can be further probed.We know that contacts between chromatin segments aredynamic; proteins that form contacts bind and dissoci-ate resulting in stochastic formation and breakage of con-tacts. This opens up interesting questions: How long dotwo beads remain in contact (looped)? When loops breakand beads diffuse away, how long does it take for the beadpairs to come back in contact? What are the factors (in-teraction strengths, polymer entropy etc.) dictating thephenomena of dynamic contacts?To study the temporal nature of chromatin, we defineloop formation time ( t ∗ L ) and contact time ( t ∗ C ) for allbead-pairs. t ∗ L is defined as the time taken for a pair ofbeads to meet ( r ∗ ij < r ∗ C ) for the first time, starting froma random equilibrium configuration. t ∗ C is defined as theduration that the bead-pairs remain looped/in contact.A schematic representation of a typical time trajectoryof 3D distance indicating t ∗ L and t ∗ C is shown in Fig. 4(a)and the actual data from our simulation, as an exam-ple, is shown in SI. Corresponding average quantities aredefined by h t ∗ L i and h t ∗ C i , respectively.Two possible factors that can influence these tempo-ral quantities are interaction strengths ( ǫ ) and polymerentropy. Since two beads having a larger segment lengthbetween them will have a higher entropy, it is expectedthat the time to come into contact is longer. In otherwords, the time of looping is expected to be dictated by r c t L t C r i j = | r i − r j | t (a) OFF GT1 GT2 SAW1.41.61.822.22.4 (b) (c)
10 20 30 40 500.40.60.811.21.41.61.82 (d) (e) (f) (g)FIG. 4:
Interaction energy and genomic seperationdictate the behavior of loop formation time h t ∗ L i andcontact time h t ∗ C i : (a) schematic representation of the dis-tance between two beads in a single trajectory showing t ∗ L and t ∗ c . (b) h t ∗ L i has a power law scaling with genomic length( h t ∗ L i ∼ s µ ) with exponent varying from 1.4 (OFF state) to2.3 (SAW) for different chromatin states. The exponents areshown in lower inset. The upper inset shows h t ∗ L i for various N values for a random walk polymer. (c) h t ∗ L i as a function ofinteraction strength with each point representing a bead pair.Note the huge spread in h t ∗ L i . Inset: h t ∗ L i binned and averagedover all bead pairs having same ǫ showing minimal influenceof ǫ . (d) h t ∗ C i as a function of s with each point representinga bead pair. Here too, note the spread. Inset: h t ∗ C i binnedand averaged over all bead pairs having the same s showingminimal dependence on the segment length. (e) h t ∗ C i increasesexponentially with the interaction strength. h t ∗ L i for all thebead pairs as a heatmap for (f) OFF and (g) SAW states. polymer entropy. To validate this hypothesis, we lookedat h t ∗ L i as a function of the genomic length with and with- out HI.As shown in Fig. 4(b) h t ∗ L i monotonically increaseswith s showing a power law behavior h t ∗ L i ∼ s µ . As acontrol, we matched our h t ∗ L i results with the previouslyknown exponents µ ≈ . µ = 2 . N = 10 , , ... ) we can infer that thedeviation from the power law for large s is due to finitechain effects (top inset). We have also computed h t ∗ L i forall the other epigenetic states revealing 2 . < µ ≤ . .
4. As we remove interactions from the sys-tem, µ gradually approaches the SAW limit. The scalingappears independent of HI as shown in SI. The changein power law may also be understood by looking at thefree energy plotted in Fig. 3(c) inset. One can see thatthe free energy has a higher tilt in the OFF state com-pared to the other states, implying that the bead-pairscan move along the landscape quicker in the OFF state.The results for h t ∗ L i suggests that even in the absence ofloop extrusion, the looping time is not too long (secondsto minutes). This also indicates that the micro phase-separation could be a viable mechanism for bringing to-gether chromatin segments and possibly explains the ex-perimentally observed fact that chromatin is functionaleven in the absence of loop extruding factors [11, 20, 23].We then examined how the interaction strength influ-ences h t ∗ L i , and found that there is a huge spread in the h t ∗ L i values, for a given ǫ (Fig. 4(c)), with the averageshowing a mild dependence on ǫ (inset).Interestingly the values of h t ∗ C i are nearly independentof genomic separation (Fig. 4(d)). Here too there is ahuge variability among different bead pairs with the insetshowing the behaviour when the segment length is aver-aged over all pairs having the same s . However, the in-teraction strength significantly alters the h t ∗ C i (Fig. 4(e))showing an exponential increase. This suggests that the h t ∗ L i is determined by the interplay between entropy (re-sulting from genomic separation) and energy (interactionstrength). Once bead-pairs come in contact h t ∗ C i is dom-inated by the interaction strength.For the OFF and SAW states, we also show h t ∗ L i be-tween all pairs of beads as a heatmap (see Fig. 4(f) &(g)). One can quickly note that the range of SAW timescales is much higher than that of the OFF state. Thisis the consequence of higher µ for the SAW compared tothe OFF state. In the SAW, one can observe that thetimes are similar for all points having the same distanceaway from the diagonal (a line parallel to diagonal axis),suggesting that what matters in this case is the interbead distance ( s ). In contrast, in the OFF state, thereis a heterogeneity and curvy color contours suggestingthat the time values are not just a function of segmentlength alone but also the identity (interaction strength)of the individual bead pairs. This once again points tothe interplay between entropy and energy. Nature of loop formation and contact time dis-tributions:
So far we have studied the average loopformation times and contact times; however, should oneassume that the average values describe these quantitiescompletely? To answer this, similar to p ( r ∗ ), here wehave investigated the nature of the distribution of thetemporal quantities. In Fig. 5(a) and (b), we presentthe probability distributions of contact ( p ( t ∗ C )) and loopformation ( p ( t ∗ L )) times, respectively. We observe that p ( t ∗ C ) ∼ exp ( − t ∗ C /τ c ) with the average time τ c that de-pends on the epigenetic state (SAW: τ c = 1 / .
6, OFF: τ c = 1 / . τ c is small for the SAW and it increasesas we add interactions to the system. However, interest-ingly, the probability of loop formation time ( t ∗ L ) has apower law decay ( p ( t ∗ L ) ∼ ( t ∗ L ) − γ ). This suggests thatthere is a huge diversity in loop formation times, and theaverage looping time alone may not be sufficient to de-scribe the loop formation phonomena. We find that theepigenetic states alter the slope of the distribution (SAW: γ = 0 .
4, OFF: γ = 1 .
0) keeping the overall nature thesame. Comparison of these two distributions reveals thatquantitatively the t ∗ L is much larger than τ ∗ c , indicatingthat chromatin segments take longer to come into contactbut stay in contact for a short time. -4 -3 -2 -1 0.1 0.2 0.3 0.4 -4 -3 (a) (b)FIG. 5: Distribution function for contact time ( t ∗ C ) and loopformation time ( t ∗ L ) for a specific bead-pair (bead 5 and bead30) in SAW and OFF state are shown in (a) and (b), respec-tively. Suggestions for experiments to test our predic-tions:
Since we take HiC-like data as input and predictpositional fluctuations and dynamics of chromatin seg-ments, microscopy is the ideal method to test our pre-dictions [10, 11, 52]. All p ( r ) predictions (Fig. 1) maybe tested either via live (without fixing) microscopy ex-periments or by collecting a large number of frozen snap-shots of segment-locations via FISH or equivalent meth-ods. Imaging experiments may also estimate the volumeoccupied by a domain (Fig. 2). From the positional fluc-tuation data, one can also obtain the effective stiffnessas described earlier in this paper (Fig. 3). To measurethe time-dependent quantities (Fig. 4, Fig. 5), apart from live microscopy experiments, one may also design appro-priate FRET pairs that can probe quantities like the con-tact time [54]. Obtaining all of these quantities for dif-ferent epigenetic states would facilitate comparison withour predictions. Conclusion
Even though there is a great improvement in our un-derstanding of static nature of chromatin organization,very little is known about the dynamics, which is a cru-cial aspect of in vivo chromatin. With the advancementof technology, it is now possible to experimentally probefluctuations and dynamics of chromatin polymer. How-ever, the main challenge to simulate dynamics of chro-matin is that we do not know the interaction strengthparameters among different segments. We have over-come this challenge by using an inverse technique andobtained optimal interaction strengths between all chro-matin segments and used it to investigate the dynamicsof a chromatin domain.We summarize our key findings: (i) The nature of 3Dspatial distances between chromatin segments in a do-main, obtained from simulations, are comparable to re-cent experimental findings. (ii) Going beyond the averageproperties, we computed the distance probability distri-bution and it shows two peaks – an interaction-drivenpeak and an entropy-dominated peak. (iii) Assumingthat interactions are arising from epigenetic states, weshow how perturbations in epigenetic states would alter p ( r ); the distance distribution between a given bead pairdepends on the interaction strength of all other pairs sug-gesting the cooperative nature of chromatin folding. (iv)Volume and the shape properties of the chromatin do-main depends on the epigenetic state. The OFF stateis highly collapsed/compact, more spherical comparedto the extended, less spherical SAW. (v) The relaxationtime of the domain is dependent on the epigenetic stateof the domain. Counter-intuitively, the relaxation timeof a highly crosslinked OFF state is much smaller thanthat of a crosslink-free SAW polymer. We explain thisphenomenon by computing effective stiffness of the do-main, from polymer fluctuations. We also show that theOFF state has a higher effective drag. (vi) We studydynamics accounting for crucial hydrodynamic interac-tions; we show that HI has a significant influence onthe relaxation time of the chromatin domain. With HI,the domain takes half the time to relax as compared tothe no-HI case. (vii) We compute the loop formationtime and the time for the looped bead pairs to remainin contact. We show that average looping time has dif-ferent scaling with genomic separation, depending on theepigenetic nature of the chromatin states. The loopingtimes show a power law distribution indicating multipletimescales that might be involved with looping. On theother hand, the contact time has an exponential distri-bution.Apart from understanding of the spatiotemporal na-ture of chromatin domains, quantities calculated herehave immense biological significance too. There is anongoing debate in the field about whether the gene regu-lation requires actual physical contact between two regu-latory segments or only the proximity would suffice. Cel-lular processes such as transport of proteins from oneregion to another (eg. enhancer-promoter), spreading ofhistone modifications in the 3D space etc would cruciallydepend on p ( r ). For example, given r , one can computethe time ( τ p ) for proteins/enzymes to diffuse from loca-tion r i to r j . The mean time would depend on the dis-tribution as h τ p i = R τ p p ( r ) dr . However, apart from thedistance among segments, the accessibility would dependon the local compactness and diffusivity too. That is,compactness of the domain (Fig. 2) and effective viscousdrag (Fig. 3) together with p ( r )(Fig. 1) would be crucialfor understanding how physics of chromatin would affectbiological function. Given that phase separation is ar-gued to be one of the important factors determining thedomain formation, our study also reveals how interplaybetween epigenetic states and polymer dynamics wouldaffect loop formation and contact times.This study can be further extended genomewide to ex-amine various gene loci and investigate fluctuations anddynamics of all domains in the genome. Such polymermodels are useful for examining aspects like spread ofhistone modifications and accessibility of the domains.We hope that this study would catalyse new experimen-tal and computational studies examining the interplaybetween epigenetics and polymer dynamics. [1] W. A. Bickmore, Annu. Rev. Genomics Hum. Genet. ,67 (2013).[2] B. D. Pope, T. Ryba, V. Dileep, F. Yue, W. Wu, O. De-nas, D. L. Vera, Y. Wang, R. S. Hansen, T. K. Canfield,et al., Nature , 402 (2014).[3] P.-C. Wei, C.-S. Lee, Z. Du, B. Schwer, Y. Zhang, J. Kao,J. Zurita, and F. W. Alt, PNAS , 1919 (2018).[4] A. Sanyal, B. R. Lajoie, G. Jain, and J. Dekker, Nature , 109 (2012).[5] E. Lieberman-Aiden, N. L. Van Berkum, L. Williams,M. Imakaev, T. Ragoczy, A. Telling, I. Amit, B. R. La-joie, P. J. Sabo, M. O. Dorschner, et al., Science ,289 (2009).[6] J. R. Dixon, S. Selvaraj, F. Yue, A. Kim, Y. Li, Y. Shen,M. Hu, J. S. Liu, and B. Ren, Nature , 376 (2012).[7] E. P. Nora, B. R. Lajoie, E. G. Schulz, L. Giorgetti,I. Okamoto, N. Servant, T. Piolot, N. L. van Berkum,J. Meisig, J. Sedat, et al., Nature , 381 (2012).[8] T. Nagano, Y. Lubling, C. V´arnai, C. Dudley, W. Le-ung, Y. Baran, N. M. Cohen, S. Wingett, P. Fraser, andA. Tanay, Nature , 61 (2017).[9] D. Ba`u, A. Sanyal, B. R. Lajoie, E. Capriotti, M. Byron, J. B. Lawrence, J. Dekker, and M. A. Marti-Renom, Nat.Struct. Mol. Biol. , 107 (2011).[10] T. Nozaki, R. Imai, M. Tanbo, R. Nagashima, S. Tamura,T. Tani, Y. Joti, M. Tomita, K. Hibino, M. T. Kanemaki,et al., Mol. Cell , 282 (2017).[11] B. Bintu, L. J. Mateo, J.-H. Su, N. A. Sinnott-Armstrong, M. Parker, S. Kinrot, K. Yamaya, A. N.Boettiger, and X. Zhuang, Science (2018).[12] Q. Szabo, A. Donjon, I. Jerkovi´c, G. L. Papadopou-los, T. Cheutin, B. Bonev, E. P. Nora, B. G. Bruneau,F. Bantignies, and G. Cavalli, Nat. Genet. , 1151(2020).[13] S. S. Rao, M. H. Huntley, N. C. Durand, E. K. Sta-menova, I. D. Bochkov, J. T. Robinson, A. L. Sanborn,I. Machol, A. D. Omer, E. S. Lander, et al., Cell ,1665 (2014).[14] M. J. Rowley and V. G. Corces, Nat. Rev. Genet. ,789 (2018).[15] H. Belaghzal, T. Borrman, A. D. Stephens, D. L. La-fontaine, S. V. Venev, Z. Weng, J. F. Marko, andJ. Dekker, Nat. Genet. pp. 1–12 (2021).[16] E. Alipour and J. F. Marko, Nucleic Acids Res. , 11202(2012).[17] M. Mir, W. Bickmore, E. E. Furlong, and G. Narlikar,Development , dev182766 (2019).[18] A. Goloborodko, J. F. Marko, and L. A. Mirny, Biophys.J. , 2162 (2016).[19] G. Fudenberg, M. Imakaev, C. Lu, A. Goloborodko,N. Abdennur, and L. A. Mirny, Cell Rep. , 2038 (2016).[20] A. Kaushal, G. Mohana, J. Dorier, I. ¨Ozdemir, A. Omer,P. Cousin, A. Semenova, M. Taschner, O. Dergai,F. Marzetta, et al., Nat. Commun. , 1 (2021).[21] S. J. Nair, L. Yang, D. Meluzzi, S. Oh, F. Yang, M. J.Friedman, S. Wang, T. Suter, I. Alshareedah, A. Gamliel,et al., Nat. Struct. Mol. Biol. , 193 (2019).[22] D. Hnisz, K. Shrinivas, R. A. Young, A. K. Chakraborty,and P. A. Sharp, Cell , 13 (2017).[23] N. S. Benabdallah, I. Williamson, R. S. Illingworth,L. Kane, S. Boyle, D. Sengupta, G. R. Grimes, P. Ther-izols, and W. A. Bickmore, Mol. Cell , 473 (2019).[24] B. A. Gibson, L. K. Doolittle, M. W. Schneider, L. E.Jensen, N. Gamarra, L. Henry, D. W. Gerlich, S. Red-ding, and M. K. Rosen, Cell , 470 (2019).[25] K. Shrinivas, B. R. Sabari, E. L. Coffey, I. A. Klein,A. Boija, A. V. Zamudio, J. Schuijers, N. M. Hannett,P. A. Sharp, R. A. Young, et al., Mol. Cell , 549 (2019).[26] Y. Zhang, R. P. McCord, Y.-J. Ho, B. R. Lajoie, D. G.Hildebrand, A. C. Simon, M. S. Becker, F. W. Alt, andJ. Dekker, Cell , 908 (2012).[27] D. Jost, P. Carrivain, G. Cavalli, and C. Vaillant, NucleicAcids Res. , 9553 (2014).[28] J. J. Parmar and R. Padinhateeri, Curr. Opin. Struct.Biol. , 111 (2020).[29] K. Kumari, B. Duenweg, R. Padinhateeri, and J. R.Prakash, Biophys. J. , 2193 (2020).[30] L. Giorgetti, R. Galupa, E. P. Nora, T. Piolot, F. Lam,J. Dekker, G. Tiana, and E. Heard, Cell , 950 (2014).[31] S. Bianco, D. G. Lupi´a˜nez, A. M. Chiariello, C. An-nunziatella, K. Kraft, R. Sch¨opflin, L. Wittler, G. An-drey, M. Vingron, A. Pombo, et al., Nat. Genet. , 662(2018).[32] G. D. Bascom, C. G. Myers, and T. Schlick, PNAS ,4955 (2019). [33] M. Di Pierro, B. Zhang, E. L. Aiden, P. G. Wolynes, andJ. N. Onuchic, PNAS , 12168 (2016).[34] A. Rosa and R. Everaers, PLoS Comput. Biol. ,e1000153 (2008).[35] M. Di Stefano, J. Paulsen, D. Jost, and M. A. Marti-Renom, Curr. Opin. Genet. Dev. , 25 (2021).[36] C. T. Clarkson, E. A. Deeks, R. Samarista, H. Ma-mayusupova, V. B. Zhurkin, and V. B. Teif, NucleicAcids Res. , 11181 (2019).[37] M. Conte, L. Fiorillo, S. Bianco, A. M. Chiariello, A. Es-posito, and M. Nicodemi, Nat. Commun. , 1 (2020).[38] Y. Qi and B. Zhang, PLoS Comput. Biol. , e1007024(2019).[39] G. Shi, L. Liu, C. Hyeon, and D. Thirumalai, Nat. Com-mun. , 1 (2018).[40] Q. MacPherson, B. Beltran, and A. J. Spakowitz, PNAS , 12739 (2018).[41] G. Bajpai and R. Padinhateeri, Biophys. J. , 207(2020).[42] C. A. Brackey, D. Marenduzzo, and N. Gilbert, NatureMethods , 767 (2020).[43] R. Bird, C. Curtiss, R. Armstrong, and O. Hassager,Dynamics of polymeric liquids, kinetic theory (volume 2)(1987). [44] T. Soddemann, B. D¨unweg, and K. Kremer, Eur. Phys.J. E , 409 (2001).[45] A. Santra, K. Kumari, R. Padinhateeri, B. D¨unweg, andJ. R. Prakash, Soft Matter , 7876 (2019).[46] J. Des Cloizeaux, J. Phys , 223 (1980).[47] R. Imai, T. Nozaki, T. Tani, K. Kaizu, K. Hibino, S. Ide,S. Tamura, K. Takahashi, M. Shribak, and K. Maeshima,Mol. Biol. Cell , 3349 (2017).[48] K. Maeshima, S. Tamura, J. C. Hansen, and Y. Itoh,Curr. Opin. Cell Biol. , 77 (2020).[49] H. Strickfaden, T. O. Tolsma, A. Sharma, D. A. Under-hill, J. C. Hansen, and M. J. Hendzel, Cell (2020).[50] T. T. Pham, M. Bajaj, and J. R. Prakash, Soft Matter , 1196 (2008).[51] T. T. Pham, B. Duenweg, and J. R. Prakash, Macro-molecules , 10084 (2010).[52] T. Germier, S. Kocanova, N. Walther, A. Bancaud, H. A.Shaban, H. Sellou, A. Z. Politi, J. Ellenberg, F. Gallardo,and K. Bystricky, Biophys. J. , 1383 (2017).[53] N. M. Toan, G. Morrison, C. Hyeon, and D. Thirumalai,J. Phys. Chem. B , 6094 (2008).[54] J. G. Yang, T. S. Madrid, E. Sevastopoulos, and G. J.Narlikar, Nat. Struct. Mol. Biol. , 1078 (2006). r X i v : . [ c ond - m a t . s o f t ] F e b Supplementary Information
Spatiotemporal organization of chromatin domains: role of interaction energy andpolymer entropy
Kiran Kumari,
1, 2, 3
J. Ravi Prakash, and Ranjith Padinhateeri IITB-Monash Research Academy, Indian Institute of Technology Bombay, Mumbai, Maharashtra - 400076, India Department of Biosciences and Bioengineering, Indian Institute of Technology Bombay, Mumbai 400076, India Department of Chemical Engineering, Monash University, Melbourne, VIC 3800, Australia
PACS numbers:
The main paper presents findings of studies on staticsand dynamics of a chromatin domain. Through chro-matin in vivo shows dynamicity and cell-to-cell variabil-ity, in literature, the nature of chromatin is mostly quan-tified through average static properties. In this work,we go beyond the average and quantify the whole phasespace by investigating probability distributions and sev-eral temporal quantities and we estimate the timescalesof loop formation and stable loop maintenance. The nec-essary supplementary information for results presentedin the main paper is provided here.This document is organised as follows : Sec. describesthe simulation technique, governing equations and In-verse Brownian dynamics (IBD) algorithm. The proce-dure to convert the non-dimensional quantities obtainedfrom simulation to its corresponding standard experi-mentally measured units is discussed in Sec. . We perturbthe interaction strengths systematically to model differ-ent epigenetic states as described in Sec. . Validation ofdistance probability distribution with SAW and cumula-tive distribution for α -globin gene is presented in Sec. .The shape analysis of α -globin domain is given in Sec. .Finally, Sec. contains the analysis of time-dependentquantities. SIMULATION DETAILS
We consider chromatin as a bead spring chain hav-ing optimal intra-chromatin interactions derived from 5Cdata using an inverse method described below. Takingexperimental contact probability data as input, we havecomputed the optimal interaction strength between dif-ferent segments of the chromatin using an Inverse Brown-ian Dynamics (IBD) method that we have developed [1].We perform Brownian dynamics simulations and com-pute various static and dynamic quantities as describedbelow. It is essential to include hydrodynamic interac-tions (HI) while computing dynamic quantities in orderto obtain even quantitatively accurate predictions [2–11].HI accounts for the propagation of velocity perturbationsdue to the motion of one segment of a polymer chain toother segments through the medium. Several studies thatattempt to probe chromatin dynamics have not consid- ered HI. In the present work, we simulate the dynamicnature of chromatin accounting for HI.
Governing equations for the bead-spring chainmodel
The total energy of the chromatin bead spring chain,made up of N beads, is U = U S + U SDK where U S is thespring potential between the adjacent beads i and ( i + 1),given by U S = X i H | r i − r i +1 | − r ) (1)where r i is the position vector of bead i , r is thenatural length and H is the stiffness of the spring.The Soddemann-Duenweg-Kremer potential ( U SDK ) isa Lennard-Jones-like potential used to mimic protein-mediated interactions [12, 13] which is given by: U SDK = "(cid:18) σr ij (cid:19) − (cid:18) σr ij (cid:19) + − ǫ ij r ij ≤ σ ǫ ij (cid:2) cos ( αr ij + β ) − (cid:3) σ ≤ r ij ≤ r c r ij ≥ r c . (2)Here r ij = | r i − r j | is the distance between beads i and j and ǫ ij is an independent parameter representing theattractive interaction strength between them. 2 / σ isthe distance at which U SDK is zero. The cut-off radius ofthe SDK potential is r c = 1 . σ . The SDK potential hasthe following advantages compared to the Lennard-Jones(LJ) potential:1. The repulsive part of the SDK potential ( r ij ≤ / σ ) is constant and remains unaffected by thechoice of the parameter ǫ ij .2. The SDK potential reaches zero at the cut off ra-dius r c = 1 . σ , unlike the LJ potential where theenergy goes to zero only at infinite distance.3. The SDK potential models the chromatin crosslinksmore efficiently [12, 14].For the simulation, all the length and time scales arenon-dimensionalised with l H = p k B T /H and λ H = ζ/ H , respectively where T is the absolute temperature, k B is the Boltzmann constant, and ζ = 6 πη s a is theStokes friction coefficient of a spherical bead of radius a ,with η s being the solvent viscosity. The evolution of beadpositions in BD simulations is governed by the followingIto stochastic differential equation, r ∗ i ( t ∗ + ∆ t ∗ ) = r ∗ i ( t ∗ )+ ∆ t ∗ D ij · ( F S ∗ j + F SDK ∗ j )+ 1 √ B ij · ∆ W j (3)Here t ∗ = t/λ H is the dimensionless time and r ∗ i = r i /l H is the dimensionless length. ∆ WWW j is a non-dimensionalWiener process with mean zero and variance ∆ t ∗ . Thebonded interactions between the beads are representedby a non-dimensional spring force, F S ∗ j , and the non-dimensional SDK force is F SDK* j . DDD ij is the diffusiontensor, defined as DDD ij = δ ij δδδ + ΩΩΩ ij , where δ ij is theKronecker delta, δδδ is the unit tensor, and ΩΩΩ ij is the hy-drodynamic interaction tensor. We use the regularizedRotne-Prager-Yamakawa (RPY) tensor to compute hy-drodynamic interactions (HI) [5], Ω ( r ∗ ) = Ω δ + Ω r ∗ r ∗ r ∗ (4)for r ∗ ≥ √ πh ∗ Ω = 3 √ π h ∗ r ∗ (cid:18) π h ∗ r ∗ (cid:19) ;Ω = 3 √ π h ∗ r ∗ (cid:18) − π h ∗ r ∗ (cid:19) (5)and for r ∗ ≤ √ πh ∗ Ω =1 − r ∗ h ∗ √ π ;Ω = 332 r ∗ h ∗ √ π (6)Here h ∗ is the dimensionless hydrodynamic parameter,defined as h ∗ = a/ ( l H √ π ). In the present work, we set h ∗ = 0 .
25 [10]. All the simulation parameters with theircorresponding values are indicated in Table S1.
Inverse Brownian Dynamics (IBD)
To simulate chromatin dynamics accurately, one re-quires the potential energy parameters that specifies theinteraction strengths ( ǫ ij ) between different regions ofchromatin. We have developed an IBD algorithm that FIG. S1: A flowchart representation of the inverse Browniandynamics algorithm. Each iteration consists of BD simula-tion, calculation of contact probability and revision of inter-action strength. Here, p (ref) represents the reference contactprobability matrix, and p (i) represents the contact probabilitymatrix from simulations at iteration i. This figure is repro-duced from Kumari et al. [1]. extracts the optimal interaction strengths between beadpairs, from an experimental contact probability map [1].The details of the algorithm is given below:The essence of the algorithm is to find the averagecontact probability and iteratively compare it to the ref-erence contact probability obtained from experiments.The average contact probability p m of the bead-pair m is given by p m = h ˆ p m i = 1 Z Z d Γ ˆ p m exp( − β H ) (7)Here Z = R d Γ exp( − β H ) is the partition function andˆ p m is an indicator function which indicates when contactoccurs between the bead pair represented by index m .As the interaction strengths are defined between a pair ofbeads, we use a single index ( m ) to define a specific beadpair ( i and j ) as m = [ i ( i − j − ( i − i variesfrom 2 to N , and j varies from 1 to ( i −
1) for a matrix ofsize N . ˆ p m is 1 if the distance between the beads is lessthan the cut-off distance of the indicator function, r ∗ c ,and 0 otherwise. We intend to target the experimentallyobtained contact probability p ref m by adjusting the well-depth of SDK attractive interactions ǫ m . The Taylorseries expansion of h ˆ p m i about the interaction strength ǫ m after neglecting higher order terms is h ˆ p m i ( ǫ m + ∆ ǫ m ) = h ˆ p m i ( ǫ m ) + X n χ mn ∆ ǫ n (8)where ∆ ǫ m is the change in the interaction strength, and TABLE S1: Numerical values of all the parameters used in the simulation
Simulation parameters symbols valuesno. of beads N r ∗ c . σ ∗ natural length of the spring r ∗ σ ∗ U SDK length scale σ ∗ h ∗ α β l H
36 nmtime scale λ H ǫ ij from IBDthe susceptibility matrix χ mn = ∂ h ˆ p m i ∂ǫ n = ∂∂ǫ n (cid:20) Z Z d Γˆ p m exp( − β H ) (cid:21) (9) χ mn = β (cid:20) Z Z ˆ p m b n exp( − β H ) d Γ − h ˆ p m i Z Z b n exp( − β H ) d Γ (cid:21) (10) χ mn = β [ h ˆ p m b n i − h ˆ p m ih b n i ]where b n = − ∂ H ∂ǫ n (11)Replacing the left hand side of Eq. 8 with the targetcontact probability p ref m obtained from experiment, weget p ref m − h ˆ p m i = X n χ mn ∆ ǫ n (12)Equation 12 can be solved for any particular iterationstep as ǫ ( i +1) n = ǫ ( i ) n + λ X m C ( i ) nm (cid:16) p ref − h ˆ p m i ( i ) (cid:17) (13)where the matrix C is the pseudo - inverse of the matrix χ [1], superscript i represents the iteration number, λ denotes the damping factor with 0 < λ <
1, and ǫ ( i +1) n isthe well-depth of the SDK attractive interaction for thenext iteration step.The flowchart of the IBD methodology is givenschematically in Fig. S1. We start with the initial guessvalues of interaction strengths, simulate the polymerfollowing the conventional forward Brownian dynamicsmethod and obtain the simulated contact probabilities in the steady state. The interaction strengths are revisedfor the next iteration, depending upon the difference be-tween the simulated and the known experimental contactprobabilities. We perform several iterations of the loop(i.e., BD simulation, calculation of contact probability,revision of interaction strength) until the error betweenthe simulated and experimental contact probabilities isless than a predetermined tolerance value. Using thisoptimal interaction strengths, we study dynamics of thechromatin domain. In the present work, we have appliedthe IBD technique to the α -globin gene locus and investi-gate the static and dynamics properties of the chromatindomain. CONVERSION OF NON-DIMENSIONALLENGTH AND TIME TO STANDARD UNITS
In our simulations, all quantities are computed in di-mensionless units as described earlier. To convert thesedimensionless numbers to standard units having appro-priate dimensions, we need to determine a lengthscaleand a timescale. By comparing our simulations with ap-propriate experimental observations, we deduce values ofcharacteristic length and time scales that can be used forthe unit conversion as follows:
Length scale:
Even though the 3D distances betweenthe genomic segments of α − globin are not available, the2D distance between the two probes located at 34 , ,
058 bp and 386 ,
139 - 425 ,
502 bp for GM12878(OFF state) was found to be 318 . ± . .
81) between the corresponding bead pair (bead 5and bead 40) from our simulation by averaging it over allthe three 2D planes ( xy , yz , zx ) as depicted in Fig. S2.By comparing 2D distance values obtained from sim-ulation and experiment, we estimate the characteristiclengthscale in our simulation as l H = 318 . / . ≈ l H to convert all non-dimensionallengths to standard units.
3D 2D "x-y" 2D "y-z" 2D ‘‘z-x"020406080100120
FIG. S2: Depiction of the calculation of 2D distance betweenbead 5 and bead 40 from simulation.
Time scale:
The timescale in our simulation is givenby: λ H = ζ H = 6 πη s a k B T (14)where H is the spring constant, T is the absolute tem-perature, k B is the Boltzmann constant, and ζ = 6 πη s a is the Stokes friction coefficient of a spherical bead ofradius a where η s is the solvent viscosity. For our prob-lem a = h ∗ l H √ π ≈ − Pa.s to 10 Pa.s. [16, 17]. Given thisdegree of variability, we decided to use a simple methodto estimate time, based on recent experimental reports ofchromatin dynamics. Chromatin segments under micro-scope seems to “diffuse” around in a region having thesize of the order of ≈ . µm ) within a timescale of ≈
50 seconds [18]. This leads to a diffusion coefficient ( D )of the order of 500 nm /s, and a timescale λ H = a D = (16 nm) ×
500 nm / s = 0 .
12 s (15)Since the calculation is to estimate the order of magni-tude number, throughout this work, we use λ H = 0 . EPIGENETIC STATES
The changes introduced into the pair-wise chromatininteractions in the simulations can be considered asequivalent to epigenetic changes. In this spirit, wehave perturbed the interaction strengths systematicallyto model different epigenetic states. The four major epi- (a) OFF state (GM12878) (b) OFF:GT1 (c) OFF:GT2
FIG. S3: (a) The interaction strength obtained form the IBDfor GM12878 α -globin gene. (b) A subset of (a) where onlythe interaction strengths greater than 1 k B T are considered.(c) A subset of (a) where the interaction strengths greaterthan 2 k B T are considered. genetic states considered in this work are as follows:1. OFF - all interactions obtained through the IBDfor α -globin gene locus in GM12878 (off/repressedstate) cell line are considered. These values areobtained in our prior work [1] and are displayed inFig. S3(a) with a heatmap representation.2. OFF:GT1 - only interactions above 1 k B T are con-sidered (see Fig. S3(b)). All other interactions( < k B T ) are taken as ǫ ij = 0 (i.e., there is onlysteric hindrance between these bead-pairs). This isa subset of the OFF state mentioned in 1.3. OFF:GT2 - only strong interactions above 2 k B T are considered (see Fig. S3(c)). All interactions be-low 2 k B T are taken as ǫ ij = 0.4. SAW - a polymer with only steric hindrance, as acontrol. In other words, all attractive interactionsare switched off.Interestingly, Fig. S3(a) and (b) do not look significantlydifferent to the eye, and even in terms of magnitude,only interactions of the order of thermal fluctuation areswitched off, yet as demonstrated in the main text, thiscauses a qualitative change to predictions. DISTANCE PROBABILITY DISTRIBUTIONS
The analytical expression given by des Cloizeaux [19]for the distance probability distribution for a self-avoiding walk (SAW) polymer is p ( r ∗ ) = C [ r ∗ ] θ +2 e − [ Kr ∗ ] − ν (16)Here ν is the Flory exponent, θ is a geometrical exponentand the coefficients C and K are given by K = Γ ([ θ + d + 2][1 − ν ]) Γ ([ θ + d ][1 − ν ]) C = 4 π [ Γ ([ θ + d + 2][1 − ν ])] θ + d [ Γ ([ θ + d ][1 − ν ])] θ + d +22 where d is the dimension. Since our simulations are in3D, d = 3. The geometrical exponent θ takes differentvalues in the following three cases1. Case 1: When both beads are the end beads of thepolymer ( θ = θ ),2. Case 2: When one of the beads is at the end andthe other bead is an intermediate bead within thechain ( θ = θ ),3. Case 3: When both the beads are intermediatebeads ( θ = θ )As the coefficients C and K depend on θ , they take dif-ferent values in each of the above cases. Following thefindings of des Cloizeaux [19], Witten and Prentis [20],Duplantier [21] and Hsu et al. [22], one can determinethat θ = 0 . θ = 0 .
461 and θ = 0 .
814 [13]. Sim-ulating a SAW polymer, we compared the probabilitydistribution for all the three cases with the correspond-ing analytical expressions using the appropriate values of θ . Fig. S4 show the validation for case 1 and 2, while thevalidation for case 3 has been shown in the main text.As can be seen, the simulations are in excellent agree-ment with the analytical expression. To the best of ourknowledge, this is the first comparison of exact numeri-cal results with the analytical expression proposed by desCloizeaux.We studied the distance probability distributions ofvarious bead-pairs revealing a distribution with twopeaks where one of the peaks is dominated by entropyof the polymer (genomic separation). The other peakemerges with an increase in the interaction strength.Fig. S5(a) depicts the Cumulative distance distribution C ( r ∗ ) for various bead-pairs at same genomic separation( s ij = 25) experiencing different interaction strengths.Differences in these plot can be easily noticed only at thesmall r ij while they looks similar overall (see inset). Thesame has been depicted for a specific bead-pair (5 ,
30) indifferent epigenetic states in Fig. S5(b). The differencein this case is not only observable for small r ij , but forthe whole regime of r ij as can be seen in the inset. FIG. S4: Comparison of distance probability distributions ob-tained from the simulation of a SAW chain with the analyticalexpression for case 1 - where both beads are end beads, andcase 2 - where one bead is at the chain end and the other beadis an intermediate bead.
SIZE AND SHAPE ANALYSIS
Here we discuss various quantities that are used tocharacterise a chromatin domain. The radius of gyra-tion of the chain, R g ≡ q h R i , where h R i is definedby h R i = h λ i + h λ i + h λ i (17)with, λ , λ , and λ being the eigenvalues of the gyrationtensor G (arranged in ascending order), with G = 12 N N X i =1 N X j =1 r ij r ij (18)Note that, G , λ , λ , and λ are calculated for eachtrajectory in the simulation before the ensemble averagesare evaluated [14, 23–28]. We examined a large ensembleof configurations and computed shape properties such asthe asphericity ( B ) and the acylindricity ( C ) B = h λ i − (cid:2) h λ i + h λ i (cid:3) (19)and C = h λ i − h λ i (20)Fig. S6(a) and (b) show the asphericity ( B ) and theacylindricity ( C ) for a chromatin domain for various epi-genetic states. Both B and C increases monotonically aswe go from the OFF state to a SAW. -4 -3 -2 -1 (a) -4 -3 -2 -1 (b) FIG. S5: (a) Cumulative distance distribution C ( r ∗ ij ) for var-ious bead-pair with the same genomic separation s ij = 25 inthe OFF state of α -globin gene in log-log scale. The sameis indicated in linear scale in the inset. (b) Comparison of C ( r ∗ ij ) for the chromatin domain under different “epigeneticstates” (see text). QUANTIFYING THE TEMPORAL NATURE OFCHROMATIN DOMAINS
To study the dynamics of chromatin domain, we quan-tified several temporal quantities such as relaxation time,loop formation time and contact time. Temporal varia-tion in properties of chromatin can be found by observ-ing spatial distance between two segments as a functionof time. Fig. S7 showcases the spatial distance betweena specific bead pair (bead 5 and bead 30) with time. Itcan be easily observed that within a single trajectory,bead-pairs come in contact several times.As a first step towards quantify dynamics, we extractedthe longest relaxation time from the end-to-end vec-tor autocorrelation function h R ∗ E (0) · R ∗ E ( t ∗ ) i / (cid:10) R ∗ (0) (cid:11) SAW LT1 LT2 OFF GT1 GT2 SAW2025303540 (a)
SAW LT1 LT2 OFF GT1 GT2 SAW4567 (b)
FIG. S6: (a) and (b) represent the asphericity and acylindric-ity, respectively. x -axis represents different interaction stateswith extreme ends representing the control (SAW) polymerand the OFF state (WT) in the middle. LT1 (LT x ) indicatethat all interactions below 1 k B T ( xk B T )are present in thepolymer. Similarly GT1 (GT x ) indicate that all interactionsabove 1 k B T ( xk B T ) are present. where R ∗ E = | r ∗ − r ∗ | . Fig. S8 displays the autocorrela-tion function decay with time not accounting for HI. Themain paper shows the same with HI and the extractedrelaxation time for both (with and without HI).We then looked at the loop formation time ( t ∗ L ) definedas the time taken for a pair of beads to meet ( r ∗ ij < r ∗ C )for the first time, starting from a random equilibriumconfiguration. We also looked at the parameters affectingaverage loop formation time ( h t ∗ L i ). Fig. S9(a) and (b)depict h t ∗ L i with genomic separation ( s ) and interactionstrength ( ǫ ), respectively. All the filled symbols in Fig. S9represent results from simulation including HI while theempty symbols represent no-HI case.The study shows that while the HI is crucial in deter-mining the relaxation time of the whole polymer, it doesnot affect temporal quantities such as h t ∗ L i . FIG. S7: Distance data from our simulation for a particularpair of beads for a three randomly chosen realisation.
FIG. S8: Exponential decay of end-to-end auto-correlationfunction with time. Here the HI is not included in the simu-lation.[1] K. Kumari, B. Duenweg, R. Padinhateeri, and J. R.Prakash, Biophys. J. , 2193 (2020).[2] T. T. Pham, M. Bajaj, and J. R. Prakash, Soft Matter , 1196 (2008).[3] T. T. Pham, B. Duenweg, and J. R. Prakash, Macro-molecules , 10084 (2010).[4] J. R. Prakash, Korea-Aust. Rheol. J. , 245 (2009).[5] J. R. Prakash, Curr. Opin. Colloid Interface Sci. , 63(2019).[6] R. Prabhakar, C. Sasmal, D. A. Nguyen, T. Sridhar, andJ. R. Prakash, Phys. Rev. Fluids , 011301 (2017).[7] R. Prabhakar and J. R. Prakash, J. Non-Newtonian FluidMech. , 163 (2004). -1 (a) (b) FIG. S9: (a) h t ∗ L i for different epinegetic states with and with-out HI. Filled symbols represent the simulation with HI andempty symbol represent the simulation without HI.[8] C. M. Schroeder, J. Rheol. , 371 (2018).[9] P. Sunthar and J. R. Prakash, Europhys. Lett. , 77(2006).[10] P. Sunthar and J. R. Prakash, Macromolecules , 617(2005).[11] C. M. Schroeder, E. S. Shaqfeh, and S. Chu, Macro-molecules , 9242 (2004).[12] T. Soddemann, B. D¨unweg, and K. Kremer, Eur. Phys.J. E , 409 (2001).[13] A. Santra, K. Kumari, R. Padinhateeri, B. D¨unweg, andJ. R. Prakash, Soft Matter , 7876 (2019).[14] M. O. Steinhauser, J. Chem. Phys. , 94901 (2005).[15] D. Ba`u, A. Sanyal, B. R. Lajoie, E. Capriotti, M. Byron,J. B. Lawrence, J. Dekker, and M. A. Marti-Renom, Nat.Struct. Mol. Biol. , 107 (2011).[16] F. Erdel, M. Baum, and K. Rippe, J. Phys.: Condens.Matter , 064115 (2015).[17] C. M. Caragine, S. C. Haley, and A. Zidovska, Phys. Rev.Lett. , 148101 (2018).[18] T. Germier, S. Kocanova, N. Walther, A. Bancaud, H. A. Shaban, H. Sellou, A. Z. Politi, J. Ellenberg, F. Gallardo,and K. Bystricky, Biophys. J. , 1383 (2017).[19] J. Des Cloizeaux, J. Phys , 223 (1980).[20] T. Witten Jr and J. Prentis, J. Chem. Phys. , 4247(1982).[21] B. Duplantier, J. Stat. Phys. , 581 (1989).[22] H.-P. Hsu, W. Nadler, and P. Grassberger, Macro-molecules , 4658 (2004).[23] W. Kuhn, Kolloid Z. , 2 (1934). [24] K. ˇSolc, J. Chem. Phys. , 335 (1971).[25] G. Zifferer, Macromol. Theory Simul. , 433 (1999).[26] C. Haber, S. A. Ruiz, and D. Wirtz, PNAS , 10792(2000).[27] D. N. Theodorou and U. W. Suter, Macromolecules ,1206 (1985).[28] M. Bishop and J. P. J. Michels, J. Chem. Phys.85