{'fraction_non_alphanumeric': 0.025298632874596903, 'fraction_numerical': 0.014375255484398419, 'mean_word_length': 4.642619169656586, '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': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
3D polymer simulations of genome organization and transcription across different chromosomes and cell types

We employ the diffusing transcription factors model for numerical simulation of chromatin topology conformations and transcriptional processes of human chromatin. Simulations of a short chromatin filament reveal different possible pathways to regulate transcription: it is shown that the transcriptional activity profile can be regulated and controlled by either acting on the chain structural properties, or on external factors, such as the number of transcription factors. Additionally, comparisons between GRO-seq experimental data and large scale numerical simulation of entire chromosomes from the human umbilical vein endothelial cell and the B-lmphocyte GM12878 show that the model provides reliable and statistically significant predictions for transcription across different cell-lines.

reliable and statistically significant predictions for transcription across different cell-lines.

Introduction
In every human cell about 2 meters of DNA filaments are packed inside a nucleus which is only ∼ 10 micrometers in diameter. 1 In order to characterise the DNA chromatin degree of spatial organization in the nucleus and its connections with the vital transcription, replication and repair cellular functions, several advanced experimental techniques, such as, just to mention a few, "Hi-C" and super resolution microscopy, 2, 3, 4, 5 have been developed. Such novel experimental tools proved fundamental in revealing complex chromatin topological features as three-dimensional (3D) conformation at different length scales. At large scales, interactions between active regions were discovered to be the driver of the genome segregation into compartments; 6 at shorter scales, chromosomes were shown to be partitioned into contiguous regions termed domains, or topologically associating domains (TADs), 7 characterised by enriched self-interactions and chromatin loops. 8 In order to better frame these in vitro experimental data, in recent years much effort was put into the theoretical development of a number of polymer and statistical physics models for DNA and chromatin, 9,10,11,12,13,14,15,16,17,18,19 able to numerically reproduce in silico the stochastic organization of genomic loci or even entire chromosomes. 9,10,11,12 Two are the commonly adopted approaches: inverse, or fitting-based, models and forward, or mechanistic, models. The inverse modelling approach is top-down: one starts from experimental data and relies on sophisticated fitting procedures to infer most possible realistic polymer models. 20,12 On the opposite, the forward modelling approach is bottom-up: 12,21 one starts from few reasonable biological and biophysical assumptions to produce computer simulations outputs directly comparable to as many types of experimental data as possible.
In combination with experimental evidence, the use of forward models resulted instrumental in providing mechanistic models for genome organization. Three are the main polymer models offered by the literature: the block copolymer model, 22,23 the loop extrusion model 24,25,26,27 and the diffusing transcription factors model. 28,29,30,10,31 In the block copolymer model, each polymer bead coarse grains a chromosome region and is assigned a specific species. Beads of the same type interact attractively so as to model direct chromatin-chromatin interactions. The loop extrusion model instead assumes the existence of active factors, or extruders, able to bind the chromatin filament at a single point and generate a growing loop by tracking along the contour of the fibre.
In this article, we focus on the the diffusing transcription factors model, in which chromosomes are modelled as bead-and-spring polymers, each bead resolving ∼ 10 3 consecutive chromosome base pairs. The polymer beads modelling the promoter or enhancer regions interact attractively with transcription factors, i.e. free beads modelling proteins generically associated with the transcription process of DNA into messenger RNA. One of the main features of the diffusing transcription factors model is the so called bridging-induced attraction feedback mechanism driving the formation and stability of clusters. 31,30,32 The latter is a fundamental ingredient allowing the model to consistently reproduce Hi-C maps of Drosophila 33 and human chromosomes, 34,31 and to highlight the emergence of a small-world gene regulatory network driven by 3D structure. 30 Although both experimental Hi-C maps and transcriptional GRO-seq data could be compared to numerical data produced by the diffusing transcription factors model, extensive comparisons have been performed in the former case but not in the latter.
The aim of the present paper is to extensively test the predictive power of the diffusing transcription factors model. First, we provide a physical characterisation of the model by simulating a short fictitious chromatin filament to show how the structural properties of the polymer chain, i.e. the specific sequence of beads types, and external factors, such as the average number of active transcription factors, affect the transcriptional activity profile. Second, we perform extensive numerical simulations of three different chromosomes, HSA 14, 18 and 19, on two different cell-lines, the human umbilical vein endothelial cell and the B-lymphocyte GM12878, to outline how the model performs compared to GRO-seq experimental data.
The structure of the paper is as follows. In section 2 we detail the model, provide technical details on simulations and present the mapping between reduced and physical units. section 3 is completely devoted to the results of our investigation. In particular, in subsection 3.1 we consider the short chromatin filament, while in subsection 3.2 we report and comment the results concerning chromosome simulations. Finally in section 4 we summarise the conclusions of our investigation.

MODEL and METHODS
The aim of the present section is to introduce the diffusing transcription factors model and report simulations details.

Diffusion transcription factors model
Chromatin fibres are coarse grained into bead-and-spring polymer chains consisting of M monomers. Each monomer represents 3000 base-pairs (bp), and has a diameter of σ 30 nm. 28 The transcription factors (TFs) are modelled as N beads of the same diameter which can bind differently to three monomer types of the DNA chain. Along the polymer there are in fact n T U high affinity transcription units (TUs), n w low affinity (weakly-sticky) beads, and n nb nonbinding polymer beads, so that M = n nb + n w + n T U . TFs can be in either an active binding or an inactive non-binding state, and are able to switch state at fixed time intervals with switching off and on rates α of f and α on , respectively, so that on average there are N a active state TFs and N p passive ones, with N = N a + N p . Active state TFs interact attractively strongly to TUs, weakly to weakly-binding beads and repulsively with non-binding beads, while passive state TFs interact repulsively with all polymer beads. Adopting a discrete index i running from 1 to N + M for both TFs and polymer beads , r i = (r ix , r iy , r iz ) will denote the position of the i-th generic bead in the three-dimensional space, while r ij = α=x,y,z (r iα − r jα ) 2 the Euclidean distance between the i-th and the j-th beads.
Concerning the polymer beads, the generic i-th and j-th beads interact purely repulsively via the Weeks-Chandler-Anderson potential: where k B is the Boltzmann constant and T is the environment temperature. In order to enforce the bead connectivity, any two consecutive chain beads interact through the harmonic potential where k h = 100k B T /σ 2 is the harmonic constant andr is the equilibrium spring distance set atr = 1.1 σ. 30 Finally, the stiffness of the chromatin fibre is taken into account by making each three consecutive beads along the polymer interact via the Kratky-Porod potential where i and j are neighbouring beads, s i is the tangent vector connecting the i-th and the i + 1-th beads, θ is the angle formed by such tangents and and l P is the persistence length of the chain. Concerning polymer beads and TFs excluded volume interactions, any polymer bead, say the a-th, and TF, say the b-th, interact via the following truncated and shifted Lennard-Jones potential, where r c is the interaction cut-off set to r c = 2 1/6 σ for active and inactive TFs and non-binding polymer beads and for inactive TFs and weakly-sticky beads and TUs (purely repulsive interactions) or to r c = 1.8σ for all active TFs and weakly-sticky beads and TUs (attractive interaction to model binding). The U LJ potential is shifted to zero at the cut-off. We set ab = k B T for purely repulsive interactions, ab = 3k B T for attractive interactions between weaklysticky beads and active TFs, and ab = 8k B T for strong attractive interactions between active TFs and TUs. A TU bead is considered transcriptionally active when it is bound to almost one TF, that is when the centres of the TU and an active TF are closer than a fixed distance threshold. Finally the global time evolution of the system is described by the solution of the following system of 3(N + M ) Langevin equations where i is the bead index running from 1 to N + M , α = x, y, z is the dimension index, m i and γ i are the mass and friction coefficient associated to the i-th bead, U i is the total potential experienced by the i-th bead and η iα (t) are a set of independent zero-mean delta-correlated stochastic white noises where δ ij and δ αβ are two Kronecker deltas over the bead and three-dimensional indices, respectively, and δ(t − t ) is a Dirac delta. For the sake of simplicity we set the mass m i ≡ m and friction coefficient γ i ≡ γ, equal for all beads.

Simulation details
Coarse grained Brownian molecular dynamics (MD) simulations are performed using the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) software package. 35 We use the Velocity-Verlet algorithm with the addition of the Langevin thermostat modelling the stochastic thermal fluctuations and viscosity of an implicit solvent. For the sake of simplicity, hydrodynamic interactions are neglected and mass m, diameter σ, temperature T , and Boltzmann constant k B are all set to unity. The typical timescales are the Lennard-Jones time τ LJ = σ m/ , with = k B T the energy unit, and the brownian time τ B = σ 2 /D, with D = k B T /γ the diffusion coefficient of a single bead. Setting γ = 1, we have τ B = τ LJ = 1.
The system of equations Equation 3 is integrated with timestep ∆t = 0.01τ B . The simulation domain is a cubic box with periodic boundary conditions of side L chosen so that the system always results dilute. Protein switching is implemented by stochastically changing each TF state every n s timesteps as prescribed by the switching off and on rates α of f and α on . The actual box side length, TF number, TFs switching rates, n s and evolution time length will be specified when describing each case.
Simulations begin with the chromatin fibre monomers initialised with a random walk, and TFs placed randomly within the simulation domain. Bead-bead overlaps are relaxed by first evolving the system for a few timesteps with a soft repulsive potential between beads. The system is then thermalized for 10 5 τ B with only repulsive Weeks-Chandler-Anderson potential interactions between all beads. Finally the system is evolved including attractive interactions and the observables of interest are monitored. A TU is considered transcriptionally active if the distance between its centre and the centre of almost one active TF is less than 2.23σ. Averages are evaluated over 100 independent runs for each case. Cluster are identified using the python-implemented DBSCAN algorithm 36 with distance threshold equal to the threshold for transcriptional activity.

RESULTS
The present section is devoted to the results of our investigation. We first report and comment the model phenomenology using a short chromatin filament, and then we focus on human chromosomes and comparisons with experimental GRO-seq data.

Short chromatin filament
As a first step, we consider a simple M = 1000 beads long polymer chain. We investigate and compare many possible configurations varying the number of TFs in the system, the non-TU bead types (either non-binding or weakly-sticky) and the TU sequence. The TU bead sequence can be chosen in two different ways: either placing a fixed number of TUs in random positions, or placing consecutive TUs at fixed bead distance. Concerning the simulation parameter choice, the simulation side box is set to L = 100σ, so that the the chain volume fraction is of order ∼ 10 −3 , while the number N of TFs is always set in such a way that their volume fraction is of order ∼ 10 −6 (the exact number will be specified case by case). The switching rates are set at α of f = α on = 10 −5 τ B , so that in the steady state N a = N p with n s = 100τ B , and the system is evolved for 8 10 5 τ B .
Reference case. The first case we considered is a chain with 39 TUs placed randomly, the remainder polymer beads set as weakly-sticky and N = 40 TFs (N a = N p = 20). This case will be used in the following as reference for comparisons and therefore will be referred to as the reference case. A cartoon of the reference case is shown in Figure 1a, along with a visual depiction of the actual polymer bead sequence. The polymer model under scrutiny is characterised by a rich topological phenomenology: 31 once a TU and a TF bind, it is more probable for other TUs to bind the same TF. The latter plays then the role of a bridge bewteen TUs and a chromatin loop forms. This makes in turn the local TU density increase, so that now it is more probable for other TFs to bind the same filament region, generating a cluster. This cycle repeats until the cluster reaches the maximum dimension allowed by the entropic cost for its formation. A visual depiction of the main stages of cluster formation is provided by Figure 1b. The phenomenon here described is a well known feedback mechanism called bridging-induced attraction. An actual typical fibre space conformation characterised by different clusters and loops is shown by Figure 1c. The insets report a close-up of the biggest cluster and underlines the involvement of a great number of TFs and TUs as well. bind strongly to TUs (red) and weakly to all other beads (blue). (b): bridging-induced attraction mechanism -cluster formation stages. (c): typical space conformation of the reference case bead polymer with active and inactive TFs (dark-red and orange, respectively) and enlargements of the boxed areas (inactive TFs not shown, weakly-binding beads depicted partially and completely transparent in the small and large inset, respectively). (d): kymograph for a single significant simulation run lasting 8 10 5 τ B . The TU index denotes the sequential number each TU appears on the chain. Black and red pixels denote a transciptionally passive or active state, respectively. The green boxes highlight the transcriptional state of TUs involved in three stable clusters. The polymer bead sequence on the right offer visual aid in locating such TUs on the entire polymer bead sequence. In our model, the vicinity of a TU and almost one active TF at a given time is used to characterise that TU as transcriptionally active at that time.
Significant observables to monitor are the transcriptional state over time of each TU (whether active or not), graphically depicted as a kymograph, and the average transcriptional activity of each TU, i.e. how long each TU has been active with respect to the total simulation time. Figure 1d shows the kymograph of a single typical run. Some of the TUs are essentially never transcribing, some others are always transcribing, and the remainder alternate between the two states. Interestingly, three continuous red stripes stand out from the black background. As highlighted by the polymer bead sequence and the green boxes, they are located in chain regions with higher local TU density, i.e. TUs are close to each other, which consistently sequester active TFs. As a consequence other TUs are prevented from binding to TFs, hence the transcriptionally passive and state-alternating TUs.
The cluster argument is supported by Figure 1e, reporting the average transcriptional activity for each TU as a function of the bead index. As highlighted by the green boxes, the activity profile shows in fact the appearance of three groups of high activity peaks, signalling the stable formation of the three aforementioned clusters in every run, alternated by low activity areas, corresponding to low transcribing TUs. Further support is provided by Figure 1f, reporting instead the transcriptional Pearson correlation matrix in which each ij , i, j = 1, . . . , n T U pixel is coloured according to the Pearson correlation between the i-th and the j-th TUs. Cluster TUs participate to the same action (being transcribed) for long time periods so, as underlined by the green squares, they result strongly positively correlated. Oppositely, as cluster TUs are involved in the sequestering dynamics, couples of TUs one inside and one outside clusters perform for long time different actions (being and not being transcribed), so they result negatively correlated.
The time series of the average number of clusters and of the average fraction of active TFs in clusters with respect to the total active ones N c a /N a reported in Figure 1g shows that the system reaches a stationary state in a relatively short time. At long times, the average cluster size slightly oscillates around 3. This value is in line with our previous identification of three highly active TUs regions involved in as many clusters. The small fluctuations are instead due to the formation and breaking of small temporary clusters. The average fraction of active TFs in clusters reaches instead a value ∼ 1. Combining the information that essentially all active TFs are occupied in transcriptional activities with the identification of three main clusters argued above, it appears clearer that TFs are actually segregated in such clusters and TUs undergo a non negligible competition to be expressed.
Variation of TFs number. As a first model perturbation, we investigated the effects produced by the variation of the TFs number with respect to the reference case with N = 40 TFs. Figure 2a, b and c report the transcriptional activity profile for the reference case with N = 10 (N a = N p = 5), N = 80 (N a = N p = 40) and N = 160 (N a = N p = 80) TFs, respectively. These panels immediately make it clear that as the number of TFs is increased, the activity profile shows the appearance of peaks in the same locations but with increasing value, up to saturation for N = 160 TFs. In order to quantify how much difference there is between the transcriptional profiles from the reference case with N = 40 TFs, and the ones with other values of N , a useful observable is the transcriptional difference d T . 30 This quantity is defined as the euclidean distance between the activity vector {a i }, i = 1, . . . , n T U in one configuration and the activity vector {a i }, i = 1, . . . , n T U in another one, in symbols Simple algebraic arguments show that the transcriptional difference is minimum with value d T = 0 when the activity of all TUs is the same, while it is maximum with value d T = √ n T U when the activity of TUs in one case are all 0 and the ones in the other case are all 1. As intuitively expected, Figure 2d shows the highest d T values to occur for the two cases with maximum and minimum number of TFs considered, while the smallest one for the case with the closest number of TFs to the reference case. The TFs number variation therefore appears as a viable pathway to alter the transcriptional behaviour of a chromatin fibre, with the most remarkable changes obtained when the number of TUs is either so low as to make the mean transcriptional activity of each TU low or negligible (see Figure 2a), or so high as to saturate it (see Figure 2c). More insight into the TF-mediated transcriptional regulation is provided by  Role of weakly-binding beads. We then investigated the role played by the weakly-binding beads. This is done by substituting the weakly-binding beads in the reference polymer chain with non-binding ones, leaving the TF switching rates, the TUs position and the TUs and TFs number unchanged. Figure 3a shows this to have an important impact on the transcriptional profile: although the highest peaks appear in the same locations as in the reference case, they are much lower, and all other previous bars are now either absent or negligible. In terms of transcriptional difference, the overall configuration is similar to Figure 2a. For the non-binding beads case we in fact have d T ∼ 1.76 10 3 and for the reference case with N = 20 TFs, where the dramatic transcription reduction was due to a 75% decrease in the TFs number, d T ∼ 1.36 10 3 . It is then natural to conjecture that the weakly-binding beads play an important role, along with the TFs, in cluster formation and stability. A first support to this hypothesis comes from the visual inspection of a number of nonbinding case configurations. We in fact observe the formation of smaller clusters, breaking and reforming during the system evolution. A more quantitative and conclusive support is provided by the average fraction of active TFs in clusters with respect to the total active ones N C a /N a and average number of clusters time series reported in Figure 3b. In contrast to the value ∼ 1 from the reference case, the average fraction of TFs in clusters reaches the constant value ∼ 0.5. In other words, on average only half of the active TFs are now involved in clusters, with a consequent much smaller average cluster dimension. The average cluster number is also strongly affected as it decreases from ∼ 3 in the reference case to ∼ 1.3. This value is larger than 1 beacuse it takes into account the cluster involving the TU with mean trancriptional activity ∼ 1 around bead index 900 plus some smaller temporary clusters quickly forming and breaking. From these data we infer that: (i) the dramatic role of weakly-binding beads lies in their ability to glue together TFs and TUs in clusters, allowing them to be larger in number and dimension and not break when some state-changing TF leaves the cluster; (ii) both TF number and chain bead type alterations result in similar and significant transcriptional difference values. Regular TU pattern. We also investigated the role of the sequence of TU beads in cluster formation and activity profile by simulating a chain whose consecutive TU beads are regularly spaced every n sp beads. The polymer length and TFs number are kept fixed at M = 1000 beads and N = 40 (N a = N p = 20) respectively. In order to simulate a case with a number of TUs close to the reference case, the TU bead distance is set at n sp = 25, resulting in n T U = 40 TUs. Figure 4a reports the average transcriptional activity for this case, which as expected now appears flat, around the value ∼ 0.5. TUs at the two ends of the filaments experience boundary effects: the left-most ones are shielded by weakly-binding beads, so it is a bit more difficult for them to access TFs, resulting in a slightly lower activity; the rightmost one is instead the very last bead of the filament, so TFs can bind it a bit more easily, hence its slightly higher activity. Interestingly, the flat transcriptional value can be modulated acting on the system parameters. Thus, setting α of f = 10 −5 τ −1 B and α on /α of f = 0.25 so that N a = 8, N p = 32 leads to ∼ 0.25, while keeping fixed the switching rates and setting N = 80 TFs leads to ∼ 0.75.
The shape of the activity profile is mirrored by the TU correlation. As shown by the Pearson correlation matrix in Figure 4b, consecutive or close TUs are positively correlated since they tend to be in the same cluster for a large fraction of time. Far apart TUs instead do not influence each other as much as in Figure 1f and result in very low reciprocal correlation values.
As the average cluster number time series in Figure 4c signals the formation of clusters, the only possibility to reconcile this observation with the flat transcriptional profile is that clusters can form with the same probability in every position along the polymer chain. This phenomenon is clearly driven by the regular TU bead sequence and supported by visual inspection of the activity profiles of single runs. In the reference case, where the only difference is the TU bead sequence, this phenomenon in fact does not occur. Here there are in fact groups of close consecutive TUs distinctly separated from other close groups of close TUs, each of them facilitating the formation of the same separate clusters in every run. We can therefore assert that the variation in the TU bead sequence is yet another pathway to regulate the topological features of the system regarding the formation of the same clusters in each run: as the TU spacing changes from uniform to random with groups of close TUs, the appearance of clusters becomes essentially deterministic, and consequently the transcripitonal profile transitions from flat to inhomogeneous, and characterised by groups of peaks.
Statistics of transcriptional activity. Finally, we make some comments about the transcriptional activity statistics. Figure 5 reports the transcriptional activity standard deviation (std) as a function of the average of the same observable for each TU for some significant cases. Data points roughly follow a parabolic trend. The reference case points span the largest std and average intervals and best depict the parabolic shape. Note that the maximum std is reached for average values ∼ 0.5, meaning that the TU resulting in such an average are those for which the std, hence transcriptional fluctuations, are the largest. Moving towards average transcriptional activity 0 and 1, the std reduces symmetrically. Together with the transcriptional profile results, this suggests that the TU bead sequence along the polymer chain induces topological conformations such that some TUs are saturated or non transcribed with a certain stability across different runs, thus having almost zero std.
The occurrence of high stds for transcriptionally fluctuating TUs is confirmed by looking at the regular TU pattern cases. As reported above, and as is clear from Figure 5, in these cases TUs on average are regularly transcribed with essentially the same parameter-dependent average transcriptional activity values. Their corresponding stds are comparatively high and consistent with the reference case TUs ones with similar average values. The characterisation for the extrema mean values 1 and 0 is instead better explained by looking at the case with N = 160 TFs and the case with non-binding beads instead of weakly-sticky. The former is in fact characterised by a transcriptional profile with almost all TUs saturated (see Figure 2c). As such TUs essentially perform very similarly in each run, their stds are clearly low, and as shown by the figure, much lower than data points with average ∼ 0.5. The hypothetical limiting case of a deterministic TU with transcriptional activity 1 in each run would lead to an std of 0, and thus complete certainty about its transcriptional outcome. In the latter case we have instead the opposite situation of most TUs non transcribed, and thus with average ∼ 0 and low fluctuations. As shown by Figure 3a, the remainder TUs are characterised by considerable average transcriptional activities and stds, which seems to hint to an approximately parabolic trend. Further simulations with non-binding beads and different number of TFs support the trend emerging from the reference case.

Human chromosomes
We now turn to the simulation of human chromosomes. In this section our goal is to validate the predictive power and reliability of the diffusing trascription factors model by comparing its numerical outputs with GRO-seq experimental data across different chromosomes and cell lines. In particular, we study HSA 14, HSA 18 and HSA 19 in human umbilical vein endothelial cell (HUVEC) line; we also study HSA 14 in B-lmphocyte GM12878 cell line. Our choice falls on HSA 14 because it is an intermediate size human chromosome and it has been the object of a number of recent studies hence is well characterised by simulations which can be used as a benchmark. 28,30 The choice for the other two HUVEC chromosomes is instead driven by their gene densities. HSA 18 is in fact an example of gene-poor chromosome as, in polymer terms, only ∼ 9% of the beads are either weakly-sticky or TUs (∼ 3% weakly-sticky, ∼ 6% TUs, see below for workflow to determine weakly-sticky beads and TUs). HSA 19 is instead an example of gene-rich chromosome, as ∼ 17% of beads are binding (∼ 6% weakly-sticky, ∼ 11% TUs). HSA 14 is instead in an intermediate ranking, as ∼ 10% of its beads are binding (∼ 3% weakly-sticky, ∼ 7% TUs).
For each chromosome, the bead polymer sequence is built following the DHS model workflow 30 schematised in Figure 6a for the GM12878 HSA14. The starting point are the ChIP-seq for H3K27ac and DNase-hypersensitivity (DHS) data provided by ENCODE. 37 The former data are optimal for showing chromosome sites that significantly correlate with open chromatin, i.e. to locate weaklybinding beads, while the latter data can be used to locate promoters and enhancers, i.e. the TUs position. Experimental data are then coarse-grained into the chosen base-pair resolution. The type of each chromosome bead is determined in the following way: beads containing a peak in H3K27ac but not in DHS are marked as weakly-sticky, beads containing a peak in DHS are marked as TUs, and all other beads are marked as non-binding. The extracted bead polymer sequence schematised in Figure 6a is used to start the simulations. The box side length L, TFs number N and polymer bead length M will be specified case by case. However we anticipate that L is always chosen so that in all cases the chain volume fraction is of order ∼ 10 −3 and the TFs volume fraction is of order ∼ 10 −4 . Note that both chromatin and TF have higher density with respect to the short chromatin filament case from subsection 3.1: although the ratio between active TF and TU numbers remains similar, the higher concentration, as detailed below, leads to some difference in the qualitative behaviour. The switching rates are set at α of f = 10 −5 τ B and α on = α of f /4 with n s = 100τ B , so that in the steady state N a /N p = 0.25. This choice is inspired by the literature 30 and reproduces the experimentally observed N a /N p ratio. 38 The system is evolved for 2 10 5 τ B .

HSA 14 for GM12878 and HUVEC cell lines
To start, we focused on simulations of the same chromosome, HSA 14, from the two cell lines GM12878 and HUVEC. The polymer is M = 35784 beads long with n T U = 2226 TUs, the simulation side box side is set to L = 570σ and the TFs number to N = 5376 TFs (N a = 1075, N p = 4301). Figure 6b reports a typical space conformation of the GM12878 chromosome (similar comments apply to the HUVEC one) formed by all possible types of polymer beads: nonbinding beads, weakly-binding beads and TUs. The configuration reproduces features similar to the ones observed for the short chromatin filament of subsection 3.1, in particular we observe the formation of chromatin loops and clusters of various dimensions driven by the bridging-induced attraction feedback mechanism. These are clearly visible from the configuration close-up reported in Figure 6c, which at the same time better highlights the gluing role played by weakly-binding beads inside clusters discussed above. Figure 7a and b report the average transcriptional activity for GM12878 and HUVEC HSA 14, respectively, along with the polymer bead sequences produced by the DHS model. Both cases show regions with zero transcriptional activity, indicating the total absence of TUs in such areas, and regions with high activity peaks, indicating the formation of clusters in such other areas. The overall transcriptional profile is similar in both cases, with slightly lower transcriptional activity values for HUVEC in regions with low local TUs density. Concerning the cluster analysis, Figure 7c, d and e show the time series of the average cluster number, of the average fraction of active TFs in clusters with respect to the total active ones N c a /N a and of the average cluster size. These plots have the same trends and reach essentially the same stationary values in both cell lines. This is a first sign of the model reliability over different cell lines. The stationary average cluster number is ∼ 50, while N c a /N a ∼ 0.4, corresponding to ∼ 430 active TFs in clusters. The information that, as also shown by Figure 7a, cluster TUs are almost transcriptionally saturated and that a large fraction of active TFs is not bound at all to them suggests that the chosen N a is such that essentially all possible potential clusters have actually formed. Consequently in this case the sequestering dynamics discussed in subsection 3.1 plays a minor role and the transcriptional competition between clusters results is diminished. Interestingly, the stationary average cluster size ∼ 10 shows a considerable amount of active TFs to be involved in each cluster.
In order to monitor the model performance, we compare our numerical results with experimental GRO-seq data. A first qualitative indication pointing towards the validity and reliability of the model is provided by Figure 8a and b, reporting a comparison between numerical transcriptional activity data and experimental GRO-seq data ranked from 0 − 100% for HUVEC and GM12878 HSA 14, respectively. As a direct comparison is not feasible the idea here is to extract from both the numerical and the experimental transcriptional data distributions the quintile values, then rank each TU numerical and experimental transcriptional activity with respect to such quintile intervals (first quintile 0 − 20%, second quintile 20 − 40%, and so on up to the fifth quintile 80 − 100%), and finally counting the occurrences of all possible couple of rankings across all TUs. Strong values over the plot diagonals would denote a significant correspondence between ranked data, and this is exactly what the two panels show. A more quantitative proof is provided by the Spearman's rank correlation coefficient 39 between the same sets of numerical and experimental data, reported in Figure 8e and found to be on average ∼ 0.32 with an estimated uncertainty of order 10 −3 in both cases (mean and standard deviation obtained using bootstrap with 100 resamplings). The correlation values themselves denote a significant correlation between numerical and experimental data, while the occurrence of essentially the same values amongst different cell lines, is yet another sign of the fact that the model performs similarly in different biological conditions.

HSA 18 and 19 from HUVEC cell line
To further test the model, we then fixed the cell line and considered different chromosomes. As aforementioned, we focused on HSA 18 and 19 from HUVEC.
For HSA 18, the polymer is M = 26026 beads long with n T U = 2226 TUs, the simulation side box side is set to L = 500σ and the TFs number to N = 1832 TFs (N a = 733, N p = 2931). For HSA 19, the polymer is instead M = 19710 beads long with n T U = 2117 TUs, the simulation side box is set to L = 556σ and the TFs number to N = 5052 (N a = 1010, N p = 4042). Figure 9a and b report the average transcriptional activity for HSA 18 and 19, respectively. As for HSA 14, both transcriptional profiles follow the TU disposition along the chain, with TUs in cluster almost saturated.
Concerning the cluster features, Figure 9 c, d and e show that also for HUVEC HSA 18 and 19 the average cluster number, N c a /N a and average cluster size time series have similar trends, in turn similar to the HSA 14 one from Figure 7. The novelty here is that while the stationary values of cluster number and average cluster size are of the same order of magnitude, there tend to be slightly more clusters but smaller in size with respect to the case of HSA 14.
Concerning the model performance, for HSA 18 the situation is essentially similar to the HSA 14 cases, even if some high off-diagonal values appear. Figure 8c shows in fact an overall good qualitative agreement between numerical and experimental data. However, the model results are validated by the Spearman's rank correlation coefficient, found to be ∼ 0.36 on average with an estimated uncertainty of order 10 −3 and graphically reported in Figure 8e. Such a measurement is completely consistent with the one concerning HSA 14, and is then a proof of the model reliability also across chromosomes. For HSA 19 the situation is instead a bit worse. We in fact notice from Figure 8d that the ranking comparison count plot is flatter with many of the off-diagonal counts similar to the diagonal ones. Moreover, the Spearman's rank correlation coefficient is now on average ∼ 0.15 with an estimated uncertainty of order 10 −3 . This result is yet statistically significant, but clearly much lower than in previous cases. A possible reason explaining this finding could be related to the DHS model mechanism. We recall that promoters directly bind proteins to initiate transcription, while enhancers can bind proteins to increase the probability that a specific gene is transcribed. In our machinery, DHS data, marking promoters and enhancers as well, are used by the DHS model to identify TU beads, whose sticky action should mimic the biological functioning of both promoters and enhancers. The discrepancy between numerical results and experimental data could therefore be traced back to the fact that our model essentially neglects the biochemical processes involved in the enhancers behaviour, treating the latter de facto as promoters. A logical consequence is that the model performance should worsen as i) the global chromosome gene density increases; ii) TUs be- come more evenly distributed. Point i) takes into consideration the mere fact that the more the TUs are, the more the ill-modelled enhancers. Point ii) takes instead into account that a close group of TUs, coarse graining promoters or enhancers or both, leads anyway to the formation of a cluster, while as seen in Figure 4, many almost evenly spaced TUs in the model would lead to the regular formation of clusters in every position of such region, neglecting the bias provided by the exact location of promoters and enhancers. For the sake of clarity, Figure 10 reports the TU counts profiles for the chromosomes of interest obtained by binning an interval of the same bead length of each polymer chromosome and counting the TU occurrences in each bin. Such profiles follow closely the average transcriptional profiles from Figure 7a and b and Figure 9a and b and show that different global TU densities are in this case associated to counts (transcriptional) profiles with different features: in the HSA 19 case, apart from the central region, the counts profile shows many consecutive high and similar counts values, while in the HSA 18 case the profile shows high counts peaks alternated by many low-counts regions. Coming back to the results, the hypothesised gene density and disposition correlation worsening seems to be confirmed by our data: the model in fact performs worst in the gene-rich and similar-counts HSA 19 case and best in the gene-poor peaked-counts HSA 18 case.

CONCLUSIONS
In this paper we have performed molecular dynamics simulations adopting the diffusing transcription factors model to study chromatin 3D conformations and transcriptional processes. The model is first employed on a short chromatin filament to highlight some generic emerging properties of the model. We then performed simulations of entire human chromosomes to prove the validity of the model predictions for transcriptional activity.
Concerning the toy short filament model, we considered various systems versions fixing the polymer bead length and varying either the number of TFs, or the polymer TUs, weakly-sticky and non-binding bead sequence. To start, we considered a reference case where the TUs are randomly positioned and the remaining polymer beads are weakly-sticky. We reproduced some known results, such as the cluster formation driven by the bridging-induced attraction feedback mechanism and a non-uniform transcriptional activity profile affected by the formation of such clusters. We then considered a number of perturbations and highlighted the main notable differences. First, we varied the TFs number in the system whilst keeping the polymer chain unchanged, and proved the transcription factor concentration is important for gene expression, as one can transition from a regime with negligible transcription to a saturation regime where all TUs are transcribed most of the time, as the TF count is increased. Second, we explored the role of the weakly-binding beads by mutating all such kinds of beads in the reference case into non-binding ones. In this way we were able to show that weakly-sticky beads play a major role in cluster formation, size and stability. Third, we showed the TUS sequence, i.e. local TU density, to have a significant impact on the overall model output. By simulating a short polymer chain with the same bead length as the reference case one and consecutive TUs regularly spaced, we obtained a flat average transcriptional profile. We traced back such a feature to the formation of clusters in random positions along the chain, with complete loss of TU-specificity in the predicted expression profile. Finally, we remarked a notable parabolic-shaped trend of the standard deviation transcriptional activity as a function of the average of the same observable. As a result of the chosen polymer bead sequence, switching rates and TFs number, some TUs have highly fluctuating transcription, while some others are forced into saturated or non-transcribing configurations characterised by very low standard deviation in transcription. Therefore, our results show that both chromatin topology and transcription can be regulated acting on both the transcription factor concentration and on the patterning of TUs and weaklybinding or non-binding beads.
Concerning the human chromosomes simulations, we considered HSA 14, 18 and 19 in HUVEC cells, and HSA 14 in GM12878 cells, so as to test the model performance on different chromosomes along the same cell-line and also along the same chromosome in different cell lines. We found that correlation between numerical and experimental data is statistically significant in all cases, in line with previous literature results. However, we observe the model to perform better as the chromosome gene density decreases. We hypothesise this to be a consequence of the chromosome chain building procedure, as it neglects some biochemical processes involved in enhancer and promoter functioning.
Future work will need to address these issues related to gene-rich chromosome simulation, and in general should be focused on further improving the correlation between numerical and experimental data. A first step could be the inclusion in the model of the difference between enhancer and promoters. Other potential avenues to explore include instead the generalisation to the case of multiple types of transcription factors and units, so as to take into account the action of different proteins on different genes, and the inclusion of dynamic modification of transcription units, for instance to take into account known repression pathways. 40