# 3D polymer simulations of genome organization and transcription across
different chromosomes and cell types

Massimiliano Semeraro  Dipartimento di Fisica, Università degli Studi di Bari,
and INFN, Sezione di Bari, via Amendola 173, 70126 Bari, Italy
<EMAIL_ADDRESS>Giuseppe Negro  SUPA, School of Physics and
Astronomy, University of Edinburgh, Peter Guthrie Tait Road, Edinburgh EH9
3FD, UK
Dipartimento di Fisica, Università degli Studi di Bari, and INFN, Sezione di
Bari, via Amendola 173, 70126 Bari, Italy
<EMAIL_ADDRESS>Antonio Suma  Dipartimento di Fisica, Università
degli Studi di Bari, and INFN, Sezione di Bari, via Amendola 173, 70126 Bari,
Italy
<EMAIL_ADDRESS>Giuseppe Gonnella  Dipartimento di Fisica, Università
degli Studi di Bari, and INFN, Sezione di Bari, via Amendola 173, 70126 Bari,
Italy
<EMAIL_ADDRESS>Davide Marenduzzo  SUPA, School of Physics and
Astronomy, University of Edinburgh, Peter Guthrie Tait Road, Edinburgh EH9
3FD, UK
<EMAIL_ADDRESS>
###### Abstract

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.

##  1 Introduction

In every human cell about $2$ meters of DNA filaments are packed inside a
nucleus which is only $\sim 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 $\sim 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.

##  2 MODEL and METHODS

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

###  2.1 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 $\sigma\simeq 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_{TU}$ high affinity transcription units (TUs), $n_{w}$ low
affinity (weakly-sticky) beads, and $n_{nb}$ non-binding polymer beads, so
that $M=n_{nb}+n_{w}+n_{TU}$. 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 $\alpha_{off}$ and $\alpha_{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 , $\textbf{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}=\sqrt{\sum_{\alpha=x,y,z}(r_{i\alpha}-r_{j\alpha})^{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:

$U^{ij}_{WCA}=\begin{cases}4k_{B}T\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}+\frac{1}{4}\right]&\text{~{}if~{}}r_{ij}<2^{1/6}\sigma\\\
0&\text{otherwise}\end{cases}~{},$ (1)

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

$U^{ij}_{harmonic}=-k_{h}\left(r_{ij}-\bar{r}\right)^{2}~{},$ (2)

where $k_{h}=100k_{B}T/\sigma^{2}$ is the harmonic constant and $\bar{r}$ is
the equilibrium spring distance set at $\bar{r}=1.1~{}\sigma$ [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

$U^{ij}_{KP}=\frac{k_{B}Tl_{P}}{\sigma}\left[1-\frac{\vec{s_{i}}\cdot\vec{s_{j}}}{|\vec{s_{i}}||\vec{s_{j}}|}\right]=k_{kp}\left[1-\cos(\theta)\right]~{},$

where $i$ and $j$ are neighbouring beads, $\vec{s_{i}}$ is the tangent vector
connecting the $i$-th and the $i+1$-th beads, $\theta$ 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,

$U^{ab}_{LJ}=\begin{cases}4\epsilon_{ab}\left[\left(\frac{\sigma}{r_{ab}}\right)^{12}-\left(\frac{\sigma}{r_{ab}}\right)^{6}-\left(\frac{\sigma}{r_{c}}\right)^{12}+\left(\frac{\sigma}{r_{c}}\right)^{6}\right]&\text{~{}if~{}}r_{ab}<r_{c}\\\
0&\text{otherwise}\end{cases}~{},$

where $r_{c}$ is the interaction cut-off set to $r_{c}=2^{1/6}\sigma$ 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\sigma$ 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 $\epsilon_{ab}=k_{B}T$ for purely repulsive
interactions, $\epsilon_{ab}=3k_{B}T$ for attractive interactions between
weakly-sticky beads and active TFs, and $\epsilon_{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

$m_{i}\frac{d^{2}r_{i\alpha}}{dt^{2}}=-\nabla
U_{i}-\gamma_{i}\frac{dr_{i\alpha}}{dt}+\sqrt{2k_{B}T\gamma_{i}}~{}\eta_{i\alpha}(t)~{},$
(3)

where $i$ is the bead index running from $1$ to $N+M$, $\alpha=x,y,z$ is the
dimension index, $m_{i}$ and $\gamma_{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 $\eta_{i\alpha}(t)$ are a set of
independent zero-mean delta-correlated stochastic white noises

$\braket{\eta_{i\alpha}(t)}=0~{},\qquad\braket{\eta_{i\alpha}(t)\eta_{j\beta}(t^{\prime})})=\delta_{ij}\delta_{\alpha\beta}\delta(t-t^{\prime})~{},$

where $\delta_{ij}$ and $\delta_{\alpha\beta}$ are two Kronecker deltas over
the bead and three-dimensional indices, respectively, and
$\delta(t-t^{\prime})$ is a Dirac delta. For the sake of simplicity we set the
mass $m_{i}\equiv m$ and friction coefficient $\gamma_{i}\equiv\gamma$, equal
for all beads.

###  2.2 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 $\sigma$, temperature $T$,
and Boltzmann constant $k_{B}$ are all set to unity. The typical timescales
are the Lennard-Jones time $\tau_{LJ}=\sigma\sqrt{m/\epsilon}$, with
$\epsilon=k_{B}T$ the energy unit, and the brownian time
$\tau_{B}=\sigma^{2}/D$, with $D=k_{B}T/\gamma$ the diffusion coefficient of a
single bead. Setting $\gamma=1$, we have $\tau_{B}=\tau_{LJ}=1$.

The system of equations Equation 3 is integrated with timestep $\Delta
t=0.01\tau_{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
$\alpha_{off}$ and $\alpha_{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}~{}\tau_{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\sigma$. 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.

###  2.3 Mapping simulation units to physical units

Simulation units can be readily converted into physical units. Using the
realistic values of bead diameter $\sigma=3~{}10^{-8}~{}m$, room temperature
$T=300~{}K$ and nucleoplasm viscosity $\eta_{sol}=10-100~{}cP$, we find
$\epsilon\simeq 4.14~{}10^{-21}~{}J$ and
$\tau_{LJ}=\tau_{B}=\sigma^{2}/D=3\pi\eta_{sol}\sigma^{3}/\epsilon\simeq
0.6-6~{}10^{-3}~{}s$.

##  3 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.

###  3.1 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\sigma$, so that the the chain
volume fraction is of order $\sim 10^{-3}$, while the number $N$ of TFs is
always set in such a way that their volume fraction is of order $\sim 10^{-6}$
(the exact number will be specified case by case). The switching rates are set
at $\alpha_{off}=\alpha_{on}=10^{-5}~{}\tau_{B}$, so that in the steady state
$N_{a}=N_{p}$ with $n_{s}=100\tau_{B}$, and the system is evolved for
$8~{}10^{5}~{}\tau_{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.

Figure 1: Reference case. (a): graphical depiction of the reference case and
visual depiction of the polymer bead sequence. Active TFs (dark-red) switching
between active/inactive states at rate $\alpha=10^{-5}~{}\tau_{B}^{-1}$ 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}\tau_{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. (e) and (f): average transcriptional activity and Pearson
correlation matrix. The green boxes highlight the cluster transcriptional
activity and correlation values. (g): cluster characterisation. The gray
horizontal lines are a guide to the eye.

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,\ldots,n_{TU}$ 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 $\sim 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.

Figure 2: Variation of TFs number. (a), (b) and (c): average transcriptional
activity for the reference case with $N=10$, $N=80$ and $N=160$ TFs,
respectively. The green and blue boxes highlight the inside- and outside-
cluster TUs we focus on in panel (e). (d): transcriptional difference between
the random case with $N=40$ TFs and the cases with $N=10$, $N=20$, $N=80$ and
$N=160$ TFs. (e): mean average transcriptional activity for inside- and
outside-cluster TUs as a function of the TFs number. Bars report the standards
deviation. Points are coloured according to the boxes in panels (a), (b) and
(c).

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^{\prime}_{i}\\},~{}i=1,\ldots,n_{TU}$ in one configuration and
the activity vector $\\{a_{i}\\},~{}i=1,\ldots,n_{TU}$ in another one, in
symbols

$d_{T}=\sqrt{\sum_{i=1}^{N}(a^{\prime}_{i}-a_{i})^{2}}~{}.$

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}=\sqrt{n_{TU}}$ 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
Figure 2e, which reports the mean transcriptional activity for the leftmost
TUs respectively inside and outside clusters (blue and green boxes in Figure
2a, b and c) as a function of the TFs number. The figure shows that both
groups of TUs increase their mean transcriptional values. Due to cluster
formation and stability, inside-cluster TUs start from a value $\sim 0.6$ and
reach the saturation value $\sim 1.0$ already at $N=40$ TFs with a concave
trend. More repressed outside-cluster TUs start instead from a negligible
value and slowly increase their mean activity until reaching saturation only
for $N=160$ TFs thanks to the large availability of TFs.

Role of weakly-binding beads.

Figure 3: Role of weakly-binding beads. (a): average transcriptional activity
for the reference case with non-binding beads instead of weakly-binding ones.
(b): cluster characterisation. The eye-aiding gray horizontal lines are the
same as in Figure 1f.

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}\sim 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}\sim 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 non-binding
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_{a}^{C}/N_{a}$ and average
number of clusters time series reported in Figure 3b. In contrast to the value
$\sim 1$ from the reference case, the average fraction of TFs in clusters
reaches the constant value $\sim 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 $\sim 3$ in the reference case to $\sim 1.3$.
This value is larger than $1$ beacuse it takes into account the cluster
involving the TU with mean trancriptional activity $\sim 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.

Figure 4: Regular TU pattern. (a): average transcriptional activity for the
regular TU pattern case. (b): Pearson correlation matrix. (d): cluster
characterisation. The gray horizontal lines are a guide to the eye.

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_{TU}=40$ TUs.

Figure 4a reports the average transcriptional activity for this case, which as
expected now appears flat, around the value $\sim 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
$\alpha_{off}=10^{-5}~{}\tau_{B}^{-1}$ and $\alpha_{on}/\alpha_{off}=0.25$ so
that $N_{a}=8,N_{p}=32$ leads to $\sim 0.25$, while keeping fixed the
switching rates and setting $N=80$ TFs leads to $\sim 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 $\sim 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 $\sim 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 $\sim 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.

Figure 5: Statistics of transcriptional activity. Plot of the standard
deviation of the transcriptional activity as a function of the average of the
same observable for some significant cases. Each point represents the couple
transcription average-standard deviation for a single TU.

###  3.2 Human chromosomes

Figure 6: GM12878 HSA $14$ configuration. (a): DHS model workflow and scheme
of the GM12878 HSA $14$ polymer bead sequence. Bars denote the bead position
and are coloured according to the legend at the bottom right corner. (b):
typical space conformation of the GM12878 HSA $14$ polymer chain. Beads are
coloured according to their types as reported by the legend (inactive TFs not
shown). (c): close-up of a portion of the entire space configuration (inactive
TFs and non-binding beads not shown).

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 $\sim 9\%$ of the beads are either weakly-sticky or TUs
($\sim 3\%$ weakly-sticky, $\sim 6\%$ TUs, see below for workflow to determine
weakly-sticky beads and TUs). HSA $19$ is instead an example of gene-rich
chromosome, as $\sim 17\%$ of beads are binding ($\sim 6\%$ weakly-sticky,
$\sim 11\%$ TUs). HSA $14$ is instead in an intermediate ranking, as $\sim
10\%$ of its beads are binding ($\sim 3\%$ weakly-sticky, $\sim 7\%$ TUs).

For each chromosome, the bead polymer sequence is built following the DHS
model workflow[30] schematised in Figure 6a for the GM12878 HSA$14$. 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 weakly-binding 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 $\sim
10^{-3}$ and the TFs volume fraction is of order $\sim 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
$\alpha_{off}=10^{-5}~{}\tau_{B}$ and $\alpha_{on}=\alpha_{off}/4$ with
$n_{s}=100\tau_{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}~{}\tau_{B}$.

####  3.2.1 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_{TU}=2226$ TUs, the simulation side box side is set to $L=570\sigma$ 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: non-binding
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 7: GM12878 and HUVEC HSA $14$. (a) and (b): average transcriptional
activity for GM12878 and HUVEC HSA $14$, respectively. The bottom part of the
panels reports the chromosome polymer bead sequence generated by the DHS
model. (c), (d) and (e): time series for average cluster number,
$N_{a}^{C}/N_{a}$ and average cluster size, for both cell types, respectively.

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_{a}^{c}/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 $\sim 50$, while $N_{a}^{c}/N_{a}\sim 0.4$, corresponding to $\sim
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
$\sim 10$ shows a considerable amount of active TFs to be involved in each
cluster.

Figure 8: Chromosomes transcriptional overview. (a), (b), (c) and
$\textbf{(d)}:$ comparison between numerical transcriptional activity data and
experimental GRO-seq data ranked from $0-100\%$ for HUVEC HSA $14$, GM12878
HSA $14$ and HUVEC HSA $18$ and $19$, respectively. (e): comparison between
transcriptional activities from simulations and GRO-seq data trough Spearman’s
rank correlation coefficient evaluation.

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 $\sim
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.

####  3.2.2 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_{TU}=2226$ TUs, the
simulation side box side is set to $L=500\sigma$ 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_{TU}=2117$ TUs, the simulation side box is set to
$L=556\sigma$ 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_{a}^{c}/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$.

Figure 9: HUVEC HSA $18$ and HSA $19$. (a) and (b): average transcriptional
activity for HUVEC HSA $18$ and HSA $19$. The bottom part of the panels
reports the chromosome polymer bead generated by the DHS model. (c), (d) and
(e): time series for average cluster number, $N_{a}^{C}/N_{a}$ and average
cluster size, for both chromosomes, respectively.

Figure 10: Chromosome TU counts profiles. (a), (b) and (c): TU counts profile
along the polymer chain for HSA $14$, HSA $18$ and HSA $19$, respectively. In
(a) we cut the initial $6800$ beads as no TUs appear there (see Figure 7a).
Small ticks denote the $200$ beads bins width.

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 $\sim 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 $\sim 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 become 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.

##  4 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 weakly-
binding 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].

## Acknowledgements

This work was possible thank to the access to Bari ReCaS e-Infrastructure
funded by MIUR through PON Research and Competitiveness 2007-2013 Call 254
Action I. We acknowledge funding from MIUR Project No. PRIN 2020/PFCXPE.

## References

  * [1] M. Krieger, M. P. Scott, P. Matsudaira, H. Lodish, J. E. Darnell, L. Zipursky, C. Kaiser, and A. Berk, “Molecular cell biology,” _San Francisco_ , 2004\.
  * [2] J. Han, Z. Zhang, and K. Wang, “3c and 3c-based techniques: the powerful tools for spatial genome organization deciphering,” _Molecular Cytogenetics_ , vol. 11, no. 1, p. 21, 2018. [Online]. Available: https://doi.org/10.1186/s13039-018-0368-2
  * [3] H. Belaghzal, T. Borrman, A. D. Stephens, D. L. Lafontaine, S. V. Venev, Z. Weng, J. F. Marko, and J. Dekker, “Liquid chromatin hi-c characterizes compartment-dependent chromatin interaction dynamics,” _Nature Genetics_ , vol. 53, no. 3, pp. 367–378, 2021. [Online]. Available: https://doi.org/10.1038/s41588-021-00784-4
  * [4] D. L. Lafontaine, L. Yang, J. Dekker, and J. H. Gibcus, “Hi-c 3.0: Improved protocol for genome-wide chromosome conformation capture.” _Curr Protoc_ , vol. 1, no. 7, p. e198, Jul 2021.
  * [5] U. J. Birk, “Super-resolution microscopy of chromatin.” _Genes (Basel)_ , vol. 10, no. 7, Jun 2019.
  * [6] E. Lieberman-Aiden, N. L. van Berkum, L. Williams, M. Imakaev, T. Ragoczy, A. Telling, I. Amit, B. R. Lajoie, P. J. Sabo, M. O. Dorschner, R. Sandstrom, B. Bernstein, M. A. Bender, M. Groudine, A. Gnirke, J. Stamatoyannopoulos, L. A. Mirny, E. S. Lander, and J. Dekker, “Comprehensive mapping of long-range interactions reveals folding principles of the human genome,” _Science_ , vol. 326, no. 5950, pp. 289–293, 2009. [Online]. Available: https://www.science.org/doi/abs/10.1126/science.1181369
  * [7] J. R. Dixon, S. Selvaraj, F. Yue, A. Kim, Y. Li, Y. Shen, M. Hu, J. S. Liu, and B. Ren, “Topological domains in mammalian genomes identified by analysis of chromatin interactions,” _Nature_ , vol. 485, no. 7398, pp. 376–380, 2012\. [Online]. Available: https://doi.org/10.1038/nature11082
  * [8] S. S. P. Rao, M. H. Huntley, N. C. Durand, E. K. Stamenova, I. D. Bochkov, J. T. Robinson, A. L. Sanborn, I. Machol, A. D. Omer, E. S. Lander, and E. L. Aiden, “A 3d map of the human genome at kilobase resolution reveals principles of chromatin looping,” _Cell_ , vol. 159, no. 7, pp. 1665–1680, 2022/10/27 2014. [Online]. Available: https://doi.org/10.1016/j.cell.2014.11.021
  * [9] M. Nicodemi and A. Pombo, “Models of chromosome structure,” _Current Opinion in Cell Biology_ , vol. 28, pp. 90–95, 2014, cell nucleus. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0955067414000465
  * [10] A. Buckle, C. A. Brackley, S. Boyle, D. Marenduzzo, and N. Gilbert, “Polymer simulations of heteromorphic chromatin predict the 3d folding of complex genomic loci,” _Molecular cell_ , vol. 72, no. 4, pp. 786–797, 2018.
  * [11] C. A. Brackey, D. Marenduzzo, and N. Gilbert, “Mechanistic modeling of chromatin folding to understand function,” _Nature Methods_ , vol. 17, no. 8, pp. 767–775, 2020. [Online]. Available: https://doi.org/10.1038/s41592-020-0852-6
  * [12] S. Bianco, A. M. Chiariello, C. Annunziatella, A. Esposito, and M. Nicodemi, “Predicting chromatin architecture from models of polymer physics,” _Chromosome Research_ , vol. 25, no. 1, pp. 25–34, 2017. [Online]. Available: https://doi.org/10.1007/s10577-016-9545-5
  * [13] M. Conte, A. M. Chiariello, A. Abraham, S. Bianco, A. Esposito, M. Nicodemi, T. Matteuzzi, and F. Vercellone, “Polymer models of chromatin imaging data in single cells,” _Algorithms_ , vol. 15, no. 9, 2022. [Online]. Available: https://www.mdpi.com/1999-4893/15/9/330
  * [14] A. Esposito, S. Bianco, L. Fiorillo, M. Conte, A. Abraham, F. Musella, M. Nicodemi, A. Prisco, and A. M. Chiariello, “Polymer models are a versatile tool to study chromatin 3D organization,” _Biochemical Society Transactions_ , vol. 49, no. 4, pp. 1675–1684, 07 2021. [Online]. Available: https://doi.org/10.1042/BST20201004
  * [15] O. Henrich, Y. A. Gutiérrez Fosado, T. Curk, and T. E. Ouldridge, “Coarse-grained simulation of dna using lammps,” _The European Physical Journal E_ , vol. 41, no. 5, p. 57, 2018. [Online]. Available: https://doi.org/10.1140/epje/i2018-11669-8
  * [16] C. Bores, O. Henrich, C. Chen, and B. M. Pettitt, “A coarse-grained representation of dna immersed in an external protein force-field,” _Biophysical Journal_ , vol. 121, no. 3, p. 209a, 2022/11/21 2022. [Online]. Available: https://doi.org/10.1016/j.bpj.2021.11.1669
  * [17] R. Natesan, K. Gowrishankar, L. Kuttippurathu, P. B. S. Kumar, and M. Rao, “Active remodeling of chromatin and implications for in vivo folding,” _The Journal of Physical Chemistry B_ , vol. 126, no. 1, pp. 100–109, 2022, pMID: 34951303. [Online]. Available: https://doi.org/10.1021/acs.jpcb.1c08655
  * [18] S. Talwar, A. Kumar, M. Rao, G. Menon, and G. Shivashankar, “Correlated spatio-temporal fluctuations in chromatin compaction states characterize stem cells,” _Biophysical Journal_ , vol. 104, no. 3, pp. 553–564, 2013. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0006349512051491
  * [19] M. Foglino, E. Locatelli, C. A. Brackley, D. Michieletto, C. N. Likos, and D. Marenduzzo, “Non-equilibrium effects of molecular motors on polymers,” _Soft Matter_ , vol. 15, pp. 5995–6005, 2019. [Online]. Available: http://dx.doi.org/10.1039/C9SM00273A
  * [20] Y. Zhan, L. Giorgetti, and G. Tiana, “Modelling genome-wide topological associating domains in mouse embryonic stem cells,” _Chromosome Research_ , vol. 25, no. 1, pp. 5–14, 2017. [Online]. Available: https://doi.org/10.1007/s10577-016-9544-6
  * [21] N. Haddad, D. Jost, and C. Vaillant, “Perspectives: using polymer modeling to understand the formation and function of nuclear compartments,” _Chromosome Research_ , vol. 25, no. 1, pp. 35–50, 2017. [Online]. Available: https://doi.org/10.1007/s10577-016-9548-2
  * [22] D. Jost, P. Carrivain, G. Cavalli, and C. Vaillant, “Modeling epigenome folding: formation and dynamics of topologically associated chromatin domains,” _Nucleic Acids Research_ , vol. 42, no. 15, pp. 9553–9561, 08 2014. [Online]. Available: https://doi.org/10.1093/nar/gku698
  * [23] D. Michieletto, E. Orlandini, and D. Marenduzzo, “Polymer model with epigenetic recoloring reveals a pathway for the de novo establishment and 3d organization of chromatin domains,” _Phys. Rev. X_ , vol. 6, p. 041047, Dec 2016. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevX.6.041047
  * [24] K. Nasmyth, “Disseminating the genome: Joining, resolving, and separating sister chromatids during mitosis and meiosis,” _Annual Review of Genetics_ , vol. 35, no. 1, pp. 673–745, 2001, pMID: 11700297. [Online]. Available: https://doi.org/10.1146/annurev.genet.35.102401.091334
  * [25] E. Alipour and J. F. Marko, “Self-organization of domain structures by DNA-loop-extruding enzymes,” _Nucleic Acids Research_ , vol. 40, no. 22, pp. 11 202–11 212, 10 2012. [Online]. Available: https://doi.org/10.1093/nar/gks925
  * [26] G. Fudenberg, M. Imakaev, C. Lu, A. Goloborodko, N. Abdennur, and L. Mirny, “Formation of chromosomal domains by loop extrusion,” _Cell Reports_ , vol. 15, no. 9, pp. 2038–2049, 2016. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2211124716305307
  * [27] A. L. Sanborn, S. S. P. Rao, S.-C. Huang, N. C. Durand, M. H. Huntley, A. I. Jewett, I. D. Bochkov, D. Chinnappan, A. Cutkosky, J. Li, K. P. Geeting, A. Gnirke, A. Melnikov, D. McKenna, E. K. Stamenova, E. S. Lander, and E. L. Aiden, “Chromatin extrusion explains key features of loop and domain formation in wild-type and engineered genomes,” _Proceedings of the National Academy of Sciences_ , vol. 112, no. 47, pp. E6456–E6465, 2015. [Online]. Available: https://www.pnas.org/doi/abs/10.1073/pnas.1518552112
  * [28] C. A. Brackley, S. Taylor, A. Papantonis, P. R. Cook, and D. Marenduzzo, “Nonspecific bridging-induced attraction drives clustering of dna-binding proteins and genome organization,” _Proceedings of the National Academy of Sciences_ , vol. 110, no. 38, pp. E3605–E3611, 2013. [Online]. Available: https://www.pnas.org/doi/abs/10.1073/pnas.1302950110
  * [29] C. A. Brackey, D. Marenduzzo, and N. Gilbert, “Mechanistic modeling of chromatin folding to understand function,” _Nature Methods_ , vol. 17, no. 8, 2020. [Online]. Available: https://doi.org/10.1038/s41592-020-0852-6
  * [30] C. A. Brackley, N. Gilbert, D. Michieletto, A. Papantonis, M. C. F. Pereira, P. R. Cook, and D. Marenduzzo, “Complex small-world regulatory networks emerge from the 3d organisation of the human genome,” _Nature Communications_ , vol. 12, no. 1, 2021. [Online]. Available: https://www.nature.com/articles/s41467-021-25875-y#citeas
  * [31] C. A. Brackley, J. Johnson, S. Kelly, P. R. Cook, and D. Marenduzzo, “Simulated binding of transcription factors to active and inactive regions folds human chromosomes into loops, rosettes and topological domains,” _Nucleic Acids Research_ , vol. 44, no. 8, pp. 3503–3512, 04 2016. [Online]. Available: https://doi.org/10.1093/nar/gkw135
  * [32] C. A. Brackley, B. Liebchen, D. Michieletto, F. Mouvet, P. R. Cook, and D. Marenduzzo, “Ephemeral protein binding to dna shapes stable nuclear bodies and chromatin domains,” _Biophysical Journal_ , vol. 112, no. 6, pp. 1085–1093, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0006349517301467
  * [33] D. Michieletto, D. Marenduzzo, and A. H. Wani, “Chromosome-wide simulations uncover folding pathway and 3d organization of interphase chromosomes,” 2016\. [Online]. Available: https://arxiv.org/abs/1604.03041
  * [34] C. A. Brackley, S. Taylor, A. Papantonis, P. R. Cook, and D. Marenduzzo, “Nonspecific bridging-induced attraction drives clustering of dna-binding proteins and genome organization,” _Proceedings of the National Academy of Sciences_ , vol. 110, no. 38, pp. E3605–E3611, 2013. [Online]. Available: https://www.pnas.org/doi/abs/10.1073/pnas.1302950110
  * [35] L. M. D. Simulator. Knuth: Computers and typesetting. [Online]. Available: https://www.lammps.org/#gsc.tab=0
  * [36] Python dbscan algorithm. [Online]. Available: earn.cluster.DBSCAN.html
  * [37] ENCODE project portal. [Online]. Available: https://www.encodeproject.org
  * [38] P. R. Cook, _Principles of nuclear structure and function_. Wiley, 2001.
  * [39] J. L. Myers, A. D. Well, and R. F. L. Jr, _Research design and statistical analysis_ , 2nd ed. Mahwah, N.J. : Lawrence Erlbaum Associates, 2011.
  * [40] S. Bustin, “Molecular biology of the cell, sixth edition; isbn: 9780815344643; and molecular biology of the cell, sixth edition, the problems book; isbn 9780815344537,” _International Journal of Molecular Sciences_ , vol. 16, no. 12, pp. 28 123–28 125, 2015. [Online]. Available: https://www.mdpi.com/1422-0067/16/12/26074