{'fraction_non_alphanumeric': 0.027837259100642397, 'fraction_numerical': 0.03389830508474576, 'mean_word_length': 4.860059549128031, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, '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}}
Combining molecular dynamics simulations and scoring method to computationally model ubiquitylated linker histones in chromatosomes

The chromatin in eukaryotic cells plays a fundamental role in all processes during a cell’s life cycle. This nucleoprotein is normally tightly packed but needs to be unpacked for expression and division. The linker histones are critical for such packaging processes and while most experimental and simulation works recognize their crucial importance, the focus is nearly always set on the nucleosome as the basic chromatin building block. Linker histones can undergo several modifications, but only few studies on their ubiquitylation have been conducted. Mono-ubiquitylated linker histones (HUb), while poorly understood, are expected to influence DNA compaction. The size of ubiquitin and the globular domain of the linker histone are comparable and one would expect an increased disorder upon ubiquitylation of the linker histone. However, the formation of higher order chromatin is not hindered and ubiquitylation of the linker histone may even promote gene expression. Structural data on chromatosomes is rare and HUb has never been modeled in a chromatosome so far. Descriptions of the chromatin complex with HUb would greatly benefit from computational structural data. In this study we generate molecular dynamics simulation data for six differently linked HUb variants with the help of a sampling scheme tailored to drive the exploration of phase space. We identify conformational sub-states of the six HUb variants using the sketch-map algorithm for dimensionality reduction and iterative HDBSCAN for clustering on the excessively sampled, shallow free energy landscapes. We present a highly efficient geometric scoring method to identify sub-states of HUb that fit into the nucleosome. We predict HUb conformations inside a nucleosome using on-dyad and off-dyad chromatosome structures as reference and show that unbiased simulations of HUb produce significantly more fitting than non-fitting HUb conformations. A tetranucleosome array is used to show that ubiquitylation can even occur in chromatin without too much steric clashes. Author summary In eukaryotic cells the linker histones play a crucial role in the formation of higher order nucleoprotein complex of DNA, especially for the arrangement of the nucleosomes. Histones can undergo several modifications, but modification of a linker histone with a single udiquitin (mono-ubiquitylation) remains one of the least understood epigenetic modifications. One reason is the inaccessibility of homogeneously modified linker histones for experimental methods, which are crucial for distinct studies. We combine molecular dynamics simulations with machine learning-based approaches to study the influence of mono-ubiquitylation in linker histones on DNA interaction and their ability to form higher order chromatin structures. We were able to determine the probable states in six differently linked histone-ubiquitin complexes via accelerating classical molecular dynamics simulations and using advanced state characterization techniques. As it is computationally unfeasible to simulate the whole chromatosome with different modified histones we developed efficient geometric scoring technique to select biologically relevant structures of all six mono-ubiquitylated linker histone that can bound to nucleosome.

Chromatin is the nucleoprotein complex in cells in which DNA is packed. Besides the 2 DNA it is composed of various histone and non-histone proteins. [1] The basic building 3 block of chromatin is the nucleosome, a complex which is formed when ≈ 147-167 bp of 4 DNA wrap around an octameric protein complex of four pairs of core histones. [2] In 5 chromatin, the nucleosome is accompanied by a linker histone and this complex is then 6 called chromatosome. The linker histone interacts with the DNA linkers which connect 7 the chromatosomes and thus directly influences higher packing of the chromatin [3], 8 such as the formation and structure of the 30 nm chromatin fiber. The nature of linker 9 histone binding modes and chromatin packing is a topic of ongoing research. [4][5][6][7][8][9][10] 10 The binding affinity of the linker histone to the chromatin can be affected by the 11 histone's post-translational modifications, such as methylation, acetylation, 12 phosphorylation, etc. [11] One such post-translational modification is the formation of 13 an isopeptide bond between a substrate protein's lysine residue and the C-terminus of 14 the protein ubiquitin (Ub) (termed as "ubiquitylation" or "ubiquitination"). 15 Ubiquitylation is known to occur at the core histones, where it can regulate, among 16 others, gene silencing, [12] chromatosome dynamics, and ultimately DNA 17 accessibility [13]. Later ubiquitylation was discovered to also occur on the linker histone 18 (Fig 1). [14][15][16] Furthermore, polyubiquitylation of the linker histone has been 19 identified as an intermediate in the repair of double strand DNA break. [17, 18] The 20 chromatin system has also been studied in silico in various resolutions ranging from 21 highly abstracted models where the nucleosome comprises a rigid body to multiscale 22 and all-atom simulations of the nucleosome core structure. [19][20][21][22][23] However, most of the 23 studies exclude the linker histone from their models, despite acknowledging the 24 importance of this protein in the formation mechanism of higher order chromatin. The 25 position of the linker histone in the "pocket" of linker and core DNA, i.e. the binding 26 modes, has been subject of discussion but most probably varies between a symmetric 27 on-dyad mode and an asymmetric off-dyad mode. [10,[24][25][26][27] Recently, all-atom 28 simulations have been conducted to further investigate the nucleosome/chromatosome   [29] B). The numbers represent aligned lysine residues that have been found to undergo ubiquitylation. [14][15][16] The red numbers indicate lysines considered in this study that have not been found to be ubiquitylated in nature (K47 in 1GHC and K81 in human H1.2). Lysines within the orange outline are investigated in this study. Abbreviations: NTD (N-terminal domain), GD (globular domain), CTD (C-terminal domain).
ubiquitylated linker histones is still not available. 40 As many cellular processes depend on chromatin accessibility, static or dynamic 41 structural data of HUb and HUb inside a chromatosome would greatly benefit from 42 further research. Structural data on the chromatosome is sparse: the protein database 43 lists 8 such structures, as of this study. These chromatosomes mainly vary in the length 44 of the linker DNA. All of them exhibit the on-dyad binding mode. We decided to use 45 the three chromatosome structures obtained via X-ray crystallography as references 46 (first three entries in Table 1). All of these structures are synthetic biological constructs 47 reconstituted from core histones and the Widom 601 DNA sequence, which has been 48 specifically tailored to increase binding to the core histones and thus has come under 49 criticism. [34,35] No structures let alone dynamic data of ubiquitylated linker histone 50 have yet been published. Attempts to crystallize HUb in a chromatosome complex have 51 not been fruitful so far.
[36] In this study, we rely on the structural data summarized in 52 Table 1  The X-ray and HADDOCK methods provide valuable first steps for structural data 58 on the chromatosome. However, they strive to find a single structure of an inherently 59 dynamic system. In this paper we aim for the combination of molecular dynamics 60 simulations with machine learning techniques to explore the conformational landscape 61 of ubiquitylated linker histone. We offer possible conformations of HUb in 62 chromatosomes and chromatosome arrays to reveal a novel view into DNA regulation by 63 ubiquitylation. We generate an extensive library of conformations of differently 64 ubiquitylated linker histones using molecular dynamics (MD) simulations guided by our 65 own enhanced sampling algorithm [

77
We make general assessments of the conformational properties of the six HUb variants 78 and propose the identified clusters as structurally stable sub-states to be used as input 79 for the interpenetration and scoring algorithm presented in Materials and Methods. are proportional to their respective distance in sketch-map space, however, they are not 92 linearly related (Fig 3 (B) and (E)). Using this method we could identify large regions 93 of unique structures for every ubiquitylation variant and small intersections of these 94 regions, where the conformational space of two or more variants overlap (Fig 3 (A)).

95
The starting sidechain dihedral angle χ 3 is the second most prominent feature (Fig 3   96 (D)). These map regions can not be taken as statistical weights for certain 97 conformations as for that a notion of density is required (Fig 3 (C)). After the removal 98 of high density regions (S1 Text) the density can be compared with the similar 99 di-Ubiquitin protein (diUb). The density map is more shallow than similar diUbs. [41] 100 This is due to the histone's conformational space which is greater, when compared to 101 ubiquitin (i.e. the histone is more flexible) and the expansion scheme driving the 102 simulation away from meta-stable, local-minima conformations. We additionally 103 questioned the completeness of the conducted sampling and stability of obtained 104 meta-stable states, but assessing ergodicity is out of the scope of this work.

105
After creating the exhaustive structural library for all six HUb variants we can 106 address more interesting questions on: how many meta-stable states does HUb assume 107 and which of these states can bind to a nucleosome.
108 Figure 3. Sketch-map projections of all conducted simulations colored by various parameters and exemplary structures. A contour plot with transparent patches colored according to ubiquitylation site (A) shows that the most distinct feature of the six HUb variants is the ubiquitylation site. A density map of the combined projections (C) shows that the low-dimensional representation is mostly shallow (blue) with only small regions of higher density (red). Coloring the combined projection according to the starting χ 3 angles ((D) and (E)) shows a finer structure of the projection. Four of the basins in (E) have been chosen and their RMSD centroids are visualized with cartoon representation and colored according to secondary structure features (B). The blue transparent sphere is focused on the histone subunit, the grey one on the Ub subunit. Closeness in sketch-map space is related to structural similarity. In (A) two regions of high density and high RMSD variance were excluded, to not bias the density towards these regions. Subfigure (C) was rendered after removing the high-denisty, high-RMSD variance region (S1 Text). Considering size and properties of the given data there are some requirements to a 116 clustering method to be able to define meta-stable states in the 2D sketch-map 117 projection: 1) handle vast amount of points; 2) work without recurring a predefined 118 number of clusters; 3) distinguish between dense basins of any size and shape 119 (meta-stable configurations) and low density background (transitions states). 120 We found the hierarchical density-based clustering algorithm HDBSCAN with 121 iterative application particularly well suited for this task. HDBSCAN transforms the 122 space before building a single linkage tree thus pushing "noise" (sparse regions) further 123 away from the dense data. We applied this algorithm iteratively, first to exclude high 124 density regions with high RMSD variance (artifacts of the projection algorithm, see S1 125 Text), and then to gradually define clusters proceeding from large to small by As a first introduction to the ISA scores we will use K30Ub and place it into the 5NL0 147 nucleosome. The sketch-map projection can now be colored by the score, revealing 148 distinct regions of higher and lower scores (Fig 5 (A)). The lowest ISA score observed 149 for placing K30Ub into 5NL0 is 8.4 and the highest is 406.2. The previously identified 150 clusters represent (meta-)stable sub-states on the protein's folding pathway and as such 151 are good candidates to identify states of lower scores which are also likely to occur in 152 real-life systems. 153 We used 4 different chromatosome structures as receptors (rows 1-3 on-dyad H1 154 binding mode, row 4 off-dyad binding mode in Table 1). Thus, the resulting scores can 155 not only be split up into the different ubiquitylation sites (Fig 6 (C)) but also into the 156 different receptor nucleosomes (Fig 6 (E)). Although the number of atoms varies for The distribution of scores is skewed to lower scores. Although the HUb proteins were 168 simulated in water without any nucleic acid present, the overall majority of sampled 169 structures fits into the nucleosome (Fig 5 (B)). This coincides with studies observing 170 reconstituted HUb chromatosomes to still form orderly complexes and not break the 171 complex completely. [31] The chromatosome with off-dyad binding linker histone 172 exhibits greater deviations because its core histones contain the flexible tails, that the 173 other chromatosomes lack. The superposed conformations of the ubiquitylated linker 174 histones can overlap with these regions. However, due to the flexibility of these tails [42] 175 we expect the tails to be able to accommodate for the ubiquitin subunit and possibly 176 interact with it in a dynamic system. ISA score Figure 5. Results of the interpenetration and scoring algorithm scoring and example structures for scoring K30HUb into 5NL0. The sketch-map projection (A) colored according to ISA score exhibits distinct regions with fitting (green, low score) and non-fitting (red, high score) conformations. The histogram of scores has a skew to lower scores (B). In (C) and (D) exemplary HUb clusters are rendered within the nucleosome of 5NL0. Both (C) and (D) visualize the DNA as tubes and the core histones as a transparent blue surface. The HUb structure bundles are rendered using their secondary structure and colored accordingly. In (C) the K30HUb cluster with the lowest ISA score is shown. Using this cluster the HUb chromatosome complex could easily be used as input for further simulations. In (D) a cluster with an ISA score of 150 is shown. This cluster can be taken as an example of a high-score cluster, that does only slightly intersects with the DNA linkers. the majority of the obtained clusters did fit surprisingly well into the chromatosomes.

182
All clusters do exhibit a non-zero mean ISA score, with the best fitting cluster (Fig 5 183 (C)) having a mean ISA score of 17. Using the four chromatosomes from Table 1 we 184 found that for 4QLC and 5WCU the minimal ISA score is achieved by the same cluster 185 (Fig 6). Interestingly, the HUb structure bundles of all on-dyad chromatosomes exhibit 186 similar spatial positioning. Here, the Ub subunit of the lowest scoring cluster is situated 187 "below" the nucleosomal core. In contrast, the off-dyad nucleosome the DNA linkers exit 188 the nucleosomal core at other angles than in the on-dyad nucleosomes. Furthermore, the 189 linker histone in the off-dyad nucleosome is slightly shifted and rotated. This leads to To model 220 such a system more research needs to be conducted. The packaging (on which the scores 221 depend) of the nucleosomes is still a topic of ongoing discussion. However, our findings 222 strongly suggest that there are native HUb conformations that fit well into the limited 223 space available in a nucleosome array.

224
Additionally to the ISA scores for the HUb variants, we looked at the position of the 225 ubiquitin with respect to DNA. For this, all sampled conformations of the six variants 226 were placed into PDB ID 5NL0 and the mean per-residue distance to the DNA used to 227 color the crystal structure PDB 1UBQ (Fig 6 (I)). The α 1 helix of the Ub-subunit 228 exhibits shorter distances (white) to the DNA than most of the residues in the β -grasp 229 region (blue).    The DNA is rendered as light blue tubes, core histones are blue translucent. Secondary structure elements are colored as default in VMD. [43] real-life HUb protein might adapt. Compared to other structural methods like NMR 242 spectroscopy or X-ray crystallography we created dynamic data of a protein modeling 243 its time-resolved conformational changes. 244 We developed an interpenetration and scoring algorithm (ISA) that exploits the 245 geometric properties of the system and shows that complex and computationally  Besides the on-dyad and off-dyad chromatosomes we applied our algorithm to a 259 tetranucleosome array. We found that, even though the DNA-pockets of the 260 tetranucleosome array are more restrictive, multiple clusters could easily be placed into 261 it without significant clashes, independent whether a on-dyad or off-dyad chromatosome 262 was used to approximate the linker histone positions.

263
Additionally, we found that the ubiquitin subunit in our simulations has the α 1 -helix 264 (beside the isopeptide bound N-terminal region) closer to the DNA than the β -grasp of 265 Ub (Fig 7 (E)). Whether the 3 lysine residues in the α 1 -helix could favorably interact 266 with the DNA is certainly a worthwhile question for further examinations of this system. 267 Our insights into how ubiquitylated linker histone might fit into chromatosomes can 268 help in further developing models for the stabilization of chromatin.

285
To reduce bias of the starting rotation of the χ 3 angle and push the systems to fully 286 explore their conformational space with all-atom MD without adding biasing potentials 287 we chose to start more rotamers of the variants by additionally varying the χ 3 sidechain 288 dihedral angle in 20°steps. All in all 6 × 18 = 108 starting structures were created. As 289 the exploratory long simulations had collapsed quite rapidly, we decided on 20 ns long 290 initial simulations from which then an expansion scheme was started that is tailored to 291 push the system to explore new regions of conformational phase space instead of getting 292 stuck in few collapsed structures.

293
This expansion scheme has already been successfully used to study the 294 conformational space of intrinsically disordered peptides.
[40] It accelerates the sampling 295 by running many parallel simulations and assessing in regular intervals whether 296 simulations have been trapped in low-energy basins. Identification of these basins is 297 done by binning a low-dimensional projection of the conformational space of the HUb 298 variants. For this reduction in dimensionality, we first needed to obtain 299 high-dimensional descriptors of HUb that represent the protein's conformation. For this 300 system we chose to use the distances between the C α -atoms to provide these 301 high-dimensional descriptors, henceforth called collective variables (CVs). A 302 conformation of the HUb can then be represented as a collection of these distances.

303
Instead of using the full set of 151 2 = 11 325 pairwise distances we chose to compare 304 two different subsets of CVs: the solvent accessible surface area (SASA)-CVs were 305 calculated according to Eq (1) by only using the C α atoms with a SASA greater than 1 306 in the respective crystal structures. The residue-wise minimal distance collective 307 variables (RMD-CVs) have been successfully employed to characterize di-Ubiquitin. [46] 308 These CVs were obtained by calculating the row-wise minimum of the pairwise distance 309 matrix between the C α atoms of the histone and ubiquitin subunit according to Eq (2). 310 K41Ub, K56Ub, and K63Ub with RMD-CVs) we were able to project the 18 × 20 ns 312 simulations into a two-dimensional projection (sometimes also called map) using the 313 sketch-map [47] method (S2 Fig, S3 Fig and S1 Text). In these projections every HDBSCAN was first used with a minimal cluster size of 1000 to exclude projection 333 artifacts (S1 Text). After the exclusion of these artifacts further clustering in the 334 2-dimensional space was carried out iteratively to obtain structurally consistent clusters. 335 HDBSCAN takes at least one hyperparameter for its clustering algorithm. This 336 parameter defines the minimum cluster size to differentiate between clusters and noise 337 points after HDBSCAN's space transformation. We reduced the hyperparameter 338 step-wise from large to small clusters as not every sub-state of the protein's phase space 339 exhibited a similar density. The minimal cluster size was set to 750, 500, 250, 125, 75, 340 50, and 25, respectively, in 8 clustering passes. After each pass the homogeneity of 341 structures in the obtained clusters was checked as follows. First, a cluster's RMSD 342 centroid was determined by using the python package MDTraj [49]. Second, the RMSD 343 between the centroid and the remaining cluster points was calculated. The same was 344 done for the distance in the distance in the two-dimensional projection. Third, a set of 345 criteria were checked: If for any given cluster at least one of the listed conditions was satisfied, this cluster 352 was not assigned at that clustering iteration and points belonging to it were put back 353 into the pool of all possible points to be resolved by subsequent (smaller cluster size) 354 clustering (an example of such a cluster is shown in Fig 8).

355
The resulting 710 clusters were then used as representative sub-states of HUb to be 356 evaluated further.

357
Fitting HUb into the chromatosome 358 To get a measure of how well a given HUb structure "fits" into a parent chromatosome 359 we developed an algorithm using the geometric properties of the HUb chromatosome histone subunit on the coordinates of the reference linker histone using RMSD 365 minimization. The Ub subunit will now possibly intersect with some parts of the 366 nucleosome because the nucleosome was not present during the simulations. The extent 367 of this clash (score) is now estimated by adding-up where the Ub subunit intersects with 368 the nucleosome. This is done by binning the atomic coordinates of the nucleosome, 369 creating a 3D histogram of atom counts. Every atom of the Ub subunit contributes to 370 the score with the value of the bin it intersects and the score is calculated as the sum of 371 the per-atom bin values.

372
An important factor to consider for this histogram (and histograms in general) is the 373 number of bins. The larger the number of bins, the more "empty space" will be present 374 in the 3D histogram. This creates a problem: The smaller the size of the bins, the more 375 zero-count bins are in the 3D histogram. However, a fine binning would be beneficial to 376 accurately model the surface of the nucleosome. To overcome this problem, we decided 377 to implement a cumulative histogram in the scoring process. Here, eight separate 378 cumulative histograms starting from the eight corners of the nucleosome's bounding box 379 are constructed. Following the diagonal of the bounding box each bin is assigned the 380 number of the contained atoms plus the values assigned to its preceding and 381 neighboring bins (see S1 Text: for examples in 1D and 2D spaces). These 8 cumulative 382 histograms are merged by assigning for every bin the respective lowest value of the 8 383 cumulative histograms. The resulting scoring histogram captures the surface of the 384 nucleosome via a fine binning, while still penalizing any Ub atoms that reach into the 385 nucleosomal core with a high score. The scoring histogram is then min-max normalized 386 onto the range [0, 1] to make values of scores more intuitive. At the end a single value is 387 obtained for a given HUb (ligand) and nucleosome (receptor) pose.

388
With this new algorithm we were able to score all 483 983 conformations sampled by 389 MD in record time (details are discussed in Results: Scoring HUb conformations into 390 nucleosome structures), so no pre-screening or "cherry-picking" had to be conducted.

391
Scores were calculated by superposing a given HUb conformation using the python 392 package MDTraj's superpose functionality and the Cα atoms of the crystal structure's 393 linker histones (Table 1). [49] If the number of Cα atoms of HUb (PDB 1GHC) did not 394 match the number of the crystal structure's linker histone a sequence alignment with 395 T-COFFEE [51] was carried out beforehand. Given the huge size of both the molecular 396 system and the dataset we designed the ISA to be particularly efficient by also 397 parallelizing it using python's built-in joblib [52] package.

398
Still, a fast algorithm needs to yield usable results, so we validate our novel and 399 purely geometric approach to scoring with tried and tested scoring functions from 400 ROSETTA [53, 54] (we used pyROSETTA). We were able to calculate ISA scores at a 401 speed of 140 poses per second. In contrast, pyROSETTA took ≈10 s per pose, which 402 makes ISA roughly 1000 times faster than pyROSETTA. Although ISA does not use 403 sophisticated potentials to model attractive and repulsive forces, a satisfying correlation 404 of 0.88 to pyROSETTA (see Fig 9).  Figure 9. Comparison of scores obtained by our interpenetration and scoring algorithm (ISA) and ROSETTA for a subset of structures. We chose a random subset of structures from the expansion scheme simulations to score using ROSETTA and ISA. Rosetta scores range from −250 to ≈ 18 000 Rosetta Energy Units (REU). ROSETTA's scoring functions with default weights support our claim that many structures exhibit lower scores, as the most structures can be found for lower scores where a high bin count can be observed (A). Our ISA algorithm correlates with the much more sophisticated ROSETTA algorithm with a Pearson correlation coefficient of 0.88 (B). However, our algorithm was 1000 times faster than pyROSETTA (both algorithms have been parallelized on a per-structure basis using the Python package joblib [52]).
Alignment of chromatosome structures into tetranucleosome. Chromatin 406 contains many chromatosomes stacked in an orderly fashion and not just a single 407 H1-nucleosome complex, so we also employed the ISA on a tetranucleosome array (PDB 408 1ZBB) published by Schalch et al. [4] This structure lacks the linker histones and thus 409 we used a sequence and RMSD alignment of 1ZBB's core histones with the core histones 410 of 5NL0 [30] and the off-dyad chromatosome structure. For this we aligned the 411 sequences of the 8 core histones with T-COFFEE and chose the longest uninterrupted 412 residue sequence for RMSD alignment using MDTraj. Some additional adjustments had 413 to be carried out to prevent the linker histone from clashing with the DNA of 1ZBB by 414 manually adjusting the rotation of the chromatosomes around the third principal axis 415 (S5 Table). Using these modelled linker histone positions we could directly use the ISA 416 to get scores for HUb poses in the tetranucleosome array.    (1)) as input CVs for sketch-map. The other 4 variants were 616 projected using the RMD-CVs (Eq (2)). Points are colored according to starting 617 χ 3 angle of the ubiquitylated lysine residue. A cyclic colormap was chosen to represent 618 the periodic nature of dihedral angles.