EXO-20-008

EXO-20-008

# Search for narrow resonances in the -tagged dijet mass spectrum in proton-
proton collisions at $\sqrt{s}=13\TeV$

###### Abstract

A search is performed for narrow resonances decaying to final states of two
jets, with at least one jet originating from a quark, in proton-proton
collisions at $\sqrt{s}=13\TeV$. The data set corresponds to an integrated
luminosity of 138collected with the CMS detector at the LHC. Jets originating
from energetic hadrons are identified through a tagging algorithm that
utilizes a deep neural network or the presence of a muon inside a jet. The
invariant mass spectrum of jet pairs is well described by a smooth
parametrization and no evidence for the production of new particles is
observed. Upper limits on the production cross section are set for excited
quarks and other resonances decaying to dijet final states containing quarks.
These limits exclude at 95% confidence level models of bosons with masses from
1.8 to 2.4and of excited quarks with masses from 1.8 to 4.0. This is the most
stringent exclusion of excited quarks to date.

##  0.1 Introduction

Searches for heavy resonances decaying into jet pairs provide a powerful tool
with the potential for discovering new physics at hadron colliders. These
resonances can be produced via parton-parton interactions, and then can decay
to two partons. The final-state partons hadronize and are observed as two
jets, referred to as dijets.

Such heavy resonances at the scale are a feature of several models that extend
the standard model (SM) to address some of its shortcomings, such as the
extreme fine tuning required in quantum corrections to accommodate a Higgs
boson observed at a mass of 125 [1, 2, 3, 4]. These models introduce heavy
resonances that couple to the SM bosons and fermions [5, 6]. A minimal
extension of the SM is represented by the sequential standard model (SSM) [7],
which introduces a spin-1 boson with the same couplings to fermions as the SM
boson, but with a much larger mass. The SSM, among many others, is generalized
in the heavy vector triplet (HVT) framework [8], which extends the SM by
introducing a triplet of heavy vector bosons, one neutral and two oppositely
charged , collectively referred to as $\PV^{\prime}$. The heavy vector bosons
couple to SM bosons and fermions with strengths
$g_{\text{H}}=g_{\text{V}}c_{\text{H}}$ and
$g_{\text{F}}=g^{2}c_{\text{F}}/g_{\text{V}}$, respectively, where
$g_{\text{V}}$ is the strength of the new interaction; $c_{\text{H}}$ is the
coupling between an HVT boson, the Higgs boson, and a longitudinally polarized
SM vector boson; $c_{\text{F}}$ is the coupling between an HVT boson and an SM
fermion; and $g$ is the electroweak coupling. In this search, we consider two
benchmarks, defined in Ref. [8]. In the Model A scenario of the HVT framework,
the coupling strengths of the heavy vector bosons to SM bosons and fermions
are of the same order, and the new particles decay primarily to fermions. In
the Model B scenario of the same framework, the couplings to fermions are
suppressed with respect to the couplings to bosons, resulting in a branching
fraction of the new particles to SM bosons that is close to unity.

Other models extending the SM postulate that leptons and quarks are composite
objects, composed of more fundamental constituents. At an energy beyond the
scale of constituent binding energies, the compositeness scale $\Lambda$, a
new interaction should emerge [9, 10, 11]. Such models foresee the existence
of excited states of quarks ($\PQq^{\ast}$) and leptons ($\ell^{\ast}$), which
could be produced in high-energy collisions and then detected through their
decays to SM particles. Excited states of composite quarks [12], which would
be produced through a gauge interaction or via contact interactions (CI) [13],
can result in large cross sections and could decay predominantly to a quark
and a gluon ($\PQq\Pg$). In the model considered in this search, the
compositeness scale $\Lambda$ is set equal to the resonance mass and the
couplings of excited quarks to other particles are assumed to be the same as
for non-excited fermions.

Since the early 1980s, searches for resonances in the dijet invariant mass
spectrum [14] have been common at hadron colliders. From proton-proton
($\Pp\Pp$) collision data at $\sqrt{s}=7$ and 8, the CMS Collaboration has
published searches for dijet resonances where at least one of the two jets in
the final state arises from a quark [15, 16, 17], based on a jet flavor
identification commonly referred to as tagging. The most recent search from
the CMS Collaboration at $\sqrt{s}=13\TeV$ [18] did not attempt to identify
the type of final-state parton that produced the jets, and considered
resonance models where the final-state partons could be gluons or any flavor
of SM quarks. Similar to a recent search from the ATLAS experiment [19], the
present search identifies heavy-flavored quarks among the two final-state
partons, and has increased sensitivity to models of dijet resonances that
decay to quarks. Examples are such models that predict an excited quark
($\PQb^{\ast}$), or a where the couplings to the third generation are enhanced
relative to the couplings to the first and second generations. Such enhanced
couplings are also favored in models created to accommodate the possible
anomalies observed in the low-energy heavy flavor sector of the SM [20, 21].

In this paper we report a search for decaying into a quark pair as well as
$\PQb^{\ast}$ decaying into quark and gluon, where two production modes are
considered for the latter: $\PQb\Pg\to\PQb^{\ast}$ and $\PQb^{\ast}$ via CI.
Both resonances are assumed to be narrow. The data used for the analysis were
collected by the CMS detector from $\Pp\Pp$ collisions at a center-of-mass
energy of 13during the LHC Run 2 (2016–2018), and correspond to a total
integrated luminosity of 138.

Tabulated results are provided in the HEPData record for this analysis [22].

##  0.2 The CMS detector

The central feature of the CMS apparatus is a superconducting solenoid of 6m
internal diameter, providing a magnetic field of 3.8T. Within the solenoid
volume are a silicon pixel and strip tracker, a lead tungstate crystal
electromagnetic calorimeter (ECAL), and a brass and scintillator hadron
calorimeter (HCAL), each composed of a barrel and two endcap sections. The
ECAL and HCAL provide pseudorapidity coverage up to $\abs{\eta}<3.0$, which is
further extended by forward calorimeters [23]. Muons are measured in gas-
ionization detectors embedded in the steel flux-return yoke outside the
solenoid. For non-isolated particles of transverse momentum $1<\pt<10\GeV$ and
$\abs{\eta}<1.4$, the reconstructed tracks have a resolution of 1.5%, a
transverse impact parameter resolution of 25–90, and a longitudinal impact
parameter resolution of 45–190 [24].

Events of interest are selected using a two-tiered trigger system. The first
level, composed of custom hardware processors, uses information from the
calorimeters and muon detectors to select events at a rate of around 100kHz
within a fixed latency of about 4 [25]. The second level, known as the high-
level trigger, consists of a farm of processors running a version of the full
event reconstruction software optimized for fast processing, and reduces the
event rate to around 1kHz before data storage [26].

A detailed description of the CMS detector, together with a definition of the
coordinate system used and the relevant kinematic variables, can be found in
Ref. [23].

##  0.3 Data and simulated samples

The data sample analyzed has been collected by the CMS experiment with trigger
algorithms that require large hadronic activity in the event. In particular,
two sets of trigger algorithms are used: single-jet triggers that require a
jet with $\pt>550\GeV$ to be present in the event and triggers requiring the
scalar sum for all jets with $\pt>30\GeV$ and $\abs{\eta}<3.0$ in the event,
to exceed 1050.

The $\qqbar\to\PZpr\to\PQb\PAQb$ signal Monte Carlo (MC) simulated events are
produced at leading order (LO) accuracy in quantum chromodynamics (QCD) with
the v2.4.2 matrix element generator [27]. Various mass hypotheses in the range
of 1.6 to 8are considered for these narrow resonance signals, which have a
natural width smaller than the experimental dijet mass resolution (2–3%). This
approximation is valid in a large fraction of the HVT parameter space, and is
fulfilled in the Model A and Model B benchmarks discussed earlier [8].

For the $\PQb^{\ast}$ production processes, the signal samples are produced
with 8.230 [28] for resonance masses up to 6. The characteristic dijet mass
shape corresponding to the $\PQb\Pg\to\PQb^{\ast}\to\PQb\Pg$ signal component
is reported in Fig. 1 for resonance masses from 1 to 6. The contribution of
the low-mass tail to the signal shape becomes more dominant at large resonance
masses. This is because the resonance natural mass shape coming from a Breit-
Wigner distribution is convolved with falling parton distribution functions
(PDFs), where for larger resonance masses the PDFs fall more steeply.

In all cases parton showering and hadronization processes are simulated by
interfacing the event generators with 8.230 with the CUETP8M1 [29, 30] tune.
The NNPDF2.3LO [31, 32] PDFs are used to model the momentum distribution of
the colliding partons inside the protons. The generated events include
additional $\Pp\Pp$ interactions (pileup) with a distribution that is chosen
to match that observed in the data. These events are processed through a full
detector simulation based on [33] and reconstructed with the same algorithms
as those used for data.

Figure 1: Simulated signal shapes of $\PQb^{\ast}$ from the process
$\PQb\Pg\to\PQb^{\ast}\to\PQb\Pg$. Shown are the wide jets (see text) used to
reconstruct the dijet mass spectra.

##  0.4 Event reconstruction and selection

The event reconstruction is performed using a particle-flow (PF) algorithm
[34], which uses an optimized combination of information from the various
elements of the CMS detector to reconstruct and identify individual particles
produced in each collision. The algorithm identifies each reconstructed
particle either as an electron, a muon, a photon, or a charged or neutral
hadron. The charged component of the pileup contribution is removed by
applying the charged hadron subtraction mitigation algorithm. The momenta of
neutral hadrons are rescaled according to their probability to originate from
the primary interaction vertex deduced from the distribution of associated
energy deposits in the calorimeters [35]. The PF candidates are then clustered
into jets using the anti-algorithm [36] implemented in the package [37, 38]
with a distance parameter $R=0.4$.

Jet energy corrections, extracted from simulation and data in multijet,
$\gamma+\text{jets}$, and $\PZ+\text{jets}$ events, are applied as a function
of and $\eta$ to correct the jet response and to account for residual
differences between data and simulation. Multijet events refer to SM events
composed uniquely of jets produced through the strong interaction. The energy
resolution for jets with $\pt\sim 1\TeV$ is approximately 5% [39]. Jets are
required to pass identification criteria, which have negligible impact on the
signal efficiency, in order to remove spurious jets arising from detector
noise [40].

Events are required to contain at least two jets, each with $\pt>30\GeV$ and
$\abs{\eta}<2.5$, and each separated by $\Delta
R=\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}>0.4$ from any electrons or muons.
Here, $\phi$ is the azimuthal angle in radians. Geometrically close jets are
combined into “wide jets”, and used to determine the dijet mass, as in our
previous dijet searches [18]. In forming wide jets, the two jets with the
largest are used as seeds, and the four-momenta of jets within $\Delta R<1.1$
of a seed jet are added to the seed jet’s four-momentum. The wide-jet
algorithm, designed for dijet resonance event reconstruction, reduces the
analysis sensitivity to gluon radiation from the final-state partons.

The pseudorapidity separation of the two wide jets is required to be
$\abs{\Delta\eta}<1.1$. This requirement suppresses $t$-channel multijet
production, and enhances the contribution to the signal from $s$-channel
production.

In order to ensure the full efficiency of the trigger, only events with a
dijet invariant mass $m_{\mathrm{jj}}$ larger than 1530are considered. The
trigger efficiency measurement is performed with a sample acquired with an
independent trigger algorithm that requires the presence of a muon with
$\pt>50\GeV$ at the high-level trigger.

Jets from quarks are identified using a tagging DeepJet discriminator that
relies on the application of convolutional neural networks on low-level
objects, such as PF candidates and secondary vertex information, and improves
the jet flavor identification capabilities of the CMS experiment, especially
at high jet momenta [41, 42]. The operation point of the discriminator was
chosen in order to maximize the sensitivity to the and $\PQb^{\ast}$ signal
searches. It corresponds to a misidentification probability of 1% for jets
from light (, , ) quarks and gluons with of about 80. For this mistagging
value, the efficiency of identifying genuine quarks is approximately 75, 45,
and 10% for a jet with of 0.2, 1.0, and 3.0, respectively. The simulated event
distributions passing the tagger are adjusted via scale factors such that they
agree with the observed data.

Events are divided into independent categories, depending on the signal model
considered. When searching for a $\PQb^{\ast}$ resonance decaying to a quark
and a gluon, only a single category with at least one -tagged jet is needed.
For a signal, three categories are defined: 2, 1and muon. Events are included
in the 2category if both jets fulfill the tagging requirements, and in the
1category if only one of the two jets is tagged. The muon category contains
the events in which neither jet passes the tag selection, but at least one jet
contains a muon, identified within both the inner tracker and the muon
detector. The presence of the muon naturally enriches the heavy-flavor
component of the jets, mitigating the loss of signal efficiency of the tagging
algorithm at very large . This category dominates the signal sensitivity for
mass larger than 5, a region where the muon presence request has a roughly 50%
larger signal efficiency than in the 2category. The product of acceptance and
efficiency of the selections for the signal were very similar in all three
data-taking years, and the average values are shown in Fig. 2. The dominant
source of inefficiency is the requirement $\abs{\Delta\eta}<1.1$, which
reduces the acceptance to approximately 41%.

Figure 2: The product of acceptance and efficiency of the event selection for
a $\PZpr\to\bbbar$ resonance as a function of the resonance mass.

##  0.5 Background estimation

Background estimates based on simulated data predict that it is largely
dominated by multijet production, which accounts for more than 95% of the
total background. The contribution of top quark pair production is
approximately 3–4% of the total background, depending on the tagging category.
Production of vector bosons in association with partons, and production of
boson pairs, contribute the remaining 1–2% of the total background.

The background is estimated directly from data assuming that it can be
described by a smooth, monotonically decreasing function. The validity of this
assumption is checked in simulation. The functions considered are power laws
of the variable $x=m_{\mathrm{jj}}/\sqrt{s}$, where $\sqrt{s}=13\TeV$.
Starting from the simplest functional form, an iterative procedure based on
the Fisher F-test [43] is used to check at 90% confidence level () if
additional parameters are needed to model the individual background
distributions. Depending on the dataset considered for the fit, a three-
parameter ($p_{0}(1-x)^{p_{1}}/x^{p_{2}}$), four-parameter
($p_{0}(1-x)^{p_{1}}/x^{p_{2}+p_{3}\log(x)}$) or five-parameter
($p_{0}(1-x)^{p_{1}}/x^{p_{2}+p_{3}\log(x)+p_{4}\log^{2}(x)}$) functional form
is necessary to describe the data. The observed dijet mass spectra are well
modeled by the background functional forms, as shown in Fig. 3, where the
widths of the dijet mass bins correspond to the dijet mass resolution.

Figure 3: The observed differential cross sections as a function of the dijet
mass, shown as fit with the background functions, for the four tagging
categories (rows) and the three data-taking periods (columns). The number of
parameters in the fit, and the goodness of fit “$\chi^{2}/\text{ndf}$”, are
listed where “ndf” is the number of degrees of freedom. The lower panel within
each row shows the pulls, $(\text{data}-\text{fit})/\text{uncertainty}$, in
units of the statistical uncertainty in data. The upper three rows are used to
search for models, the bottom row is used to search for the $\PQb^{\ast}$
model, and example shapes of these signal models are shown with the same
arbitrary normalization for three choices of resonance mass.

The expected shape of the reconstructed signal mass distribution is extracted
from the simulated signal samples. The $\PQb^{\ast}$ signal shapes are dijet
mass distributions obtained directly from the simulation, as shown in Fig. 1,
while the signal shapes are parametrizations that are fit separately to the
simulation samples for each category, using a double-sided Crystal Ball
function [44]. This functional form consists of a Gaussian distribution that
models the core of the shape, and two power law functions that model the upper
and lower tails, and requires a total of 6 parameters. The resolution of the
reconstructed is given by the width of the Gaussian core, and is found to be
constant at 2% of the resonance mass. The signal shape, for a of arbitrary
mass, is then obtained by a spline interpolation of the simulated signal shape
parameter values as a function of mass. Examples of the and $\PQb^{\ast}$
shapes are shown in Fig. 3, with arbitrary normalizations.

Statistical tests have been performed to check the robustness of the fit
method. Pseudo-experiments are generated after injecting a simulated signal
with a range of values for the signal mass and cross section. The dijet mass
distribution of the pseudo-data is then fitted with the signal distribution
combined with an alternate background function containing one more parameter
than the nominal function. The fitted signal yield is found to be compatible
within one third of the statistical uncertainty to the injected yield,
regardless of the injected signal strength and resonance mass. These tests
confirm that the background estimation is insensitive both to the choice of
the function used and to the possible presence of a signal in the data.

##  0.6 Systematic uncertainties

The background shape is estimated from the fit to the data in the considered
categories, and no assumption or constraint is applied to function parameters.
These parameters are considered as uncorrelated among categories and data-
taking periods. As described in Section 0.5, different parametrizations of the
fit function were studied and introduce no additional source of uncertainty.

The dominant uncertainties in the signal arise from the tagging and the jet
reconstruction uncertainties. The uncertainty in the tagging scale factors
[45] yields an uncertainty in the signal normalization of 2 and 4% in the 1and
${\geq}1\PQb$ categories, respectively, and 15% in the 2and muon categories.
Uncertainties in the reconstruction of the hadronic jets affect mainly the
shape of the reconstructed resonance mass. The four-momenta of the
reconstructed jets are scaled and smeared according to the uncertainties in
the jet and momentum resolution. These effects result in a 2% uncertainty in
the mean, and 10% in the width of the signal. The signal normalization is also
affected by the uncertainty in the selection of the muons inside the jets. The
efficiency of the identification of the muons inside jets is measured in a
statistically independent and almost pure sample of high- jets, originating
from the decay of pair-produced top quarks. This uncertainty is deduced from
the statistical uncertainty in the efficiency measurement, ranging up to 43%
for muons with $\pt<100\GeV$. An uncertainty of 100% is assumed for jets with
the of an associated muon beyond this 100threshold, because no data was
available in this region. Additional systematic uncertainties affecting the
signal normalization include the vetoes for high-leptons and missing
transverse momenta (accounting for 1% each), pileup contributions (0.1%), and
the integrated luminosity (1.2% in 2016 [46], 2.3% in 2017 [47], and 2.5% in
2018 [48]). The systematic uncertainty from the choice of PDFs [49] is
estimated to be 8–41% of the normalization of the signal cross section,
depending on the resonance mass. The factorization and renormalization scale
uncertainties are estimated by varying the scales up and down by a factor of
2, both simultaneously and independently, using the maximum obtained value.
The resulting effect is a variation of 6–14% of the normalization of the
signal cross section.

##  0.7 Results and interpretation

Exclusion limits are obtained by performing a background-only fit and a
combined signal-plus-background fit to the dijet mass distributions,
separately by category and data-taking year. In the fit, based on a profile
likelihood, the parameters and the normalization of the background in each
category are left free to float. The systematic uncertainties in the signal
are treated as nuisance parameters, with Gaussian constraints for the jet
energy scale and resolution, and log-normal constraints for the integrated
luminosity, and are profiled in the statistical interpretation [50]. The
uncertainties that affect the signal normalization (PDFs and factorization and
renormalization scales) are treated differently depending on how the exclusion
is presented. When deriving upper limits on the cross section, these
uncertainties are not varied in the fit, but are reported separately as the
uncertainty in the theoretical cross sections from the model. When placing
limits on the model parameters, these nuisance parameters are fixed at the
best-fit values, in the same manner as with the other systematic
uncertainties. A more detailed description of the statistical treatment and
the likelihood function is reported in Ref. [51]. Upper limits at 95% are set
using the modified frequentist method [52, 53], adopting the asymptotic
approximation [54].

A model-independent representation of the observed upper limit on the product
of cross section and branching fraction of a resonance, as well as the
expected limit and its relative 68 and 95% uncertainty bands, are shown in
Fig. 4. The acceptance, which is included in the product on the right side of
Fig. 4, is defined exclusively via the geometric requirements $\pt>30\GeV$,
$\abs{\eta}<2.5$, and $\abs{\Delta\eta}<1.1$. This analysis sets an observed
mass limit of 2.4, and an expected mass limit of 2.3, at 95% , on a narrow
resonance in both the SSM and the HVT Model A. No upper limit can be set on
Model B with suppressed fermionic coupling.

Figure 4: The observed 95% upper limits (solid curve) on the product of the
cross section and branching fraction (left), and multiplied by signal
acceptance (right), for a resonance decaying to . The corresponding expected
limits (dashed curve) and their variations at the 1- and 2-standard deviation
levels (shaded bands) are also shown. Limits are compared to predicted cross
sections for bosons from the sequential SM (SSM) and the heavy vector triplet
(HVT) models A and B. The latter two models follow the parameter choices of
$g_{\text{V}}=1$ and $g_{\text{V}}=3$ respectively.

The cross section limit shown in Fig. 4 (left), is reinterpreted in Fig. 5 as
limits on the coupling strengths of a heavy vector boson to SM bosons
($g_{\text{H}}$) and fermions ($g_{\text{F}}$). The subset of the parameter
space where the natural width of the resonance is larger than the typical
experimental resolution, and hence the narrow width approximation is invalid,
is indicated in Fig. 5 with a gray shaded area.

Figure 5: The coupling strengths to SM bosons ($g_{\text{H}}$) and fermions
($g_{\text{F}}$) of a boson with mass 2.0(blue) and 2.5(red) that are excluded
at 95% for the HVT model. The shading indicates the excluded side of the
contour. The benchmark scenarios corresponding to HVT models A and B are
represented by a purple cross and a red point, respectively. The gray shaded
area corresponds to the region where the resonance natural width is predicted
to be larger than the typical experimental resolution, and thus the narrow-
width approximation is not fulfilled.

For resonances that decay to a quark and a gluon, upper limits are set by
using the category where at least one of the two leading jets is tagged. Upper
limits on the joint product of cross section, branching fraction, and
acceptance are reported in Fig. 6. When considering the single $\PQb^{\ast}$
production process, $\PQb\Pg\to\PQb^{\ast}\to\PQb\Pg$, with the cross section
shown in Fig. 6, this analysis sets 95% mass limits of 2.5(observed) and
2.6(expected) for an excited quark. Additionally the $\PQb^{\ast}$ production
via contact interactions has been considered. Several processes contribute to
the $\PQb^{\ast}$ production via CI [13], the dominant one being
$\qqbar\to\PAQb\PQb^{\ast}(\text{or}\,\PQb\PAQb^{\ast})\to\PAQb\PQb\Pg$ where
three jets are present in the final state. In order to display the theory
cross section for this new process on the same exclusion limit plot, only its
resonant component has been considered. In particular events are selected when
the reconstructed dijet originates from the $\PQb^{\ast}$ decay, and does not
include the other quark in the event, which is the case for 50 (80)% of the
events within the acceptance at 2 (4). Including both processes in the total
signal cross section (single $\PQb^{\ast}$ production and the resonant
component of the $\PQb^{\ast}$ production via CI) gives a significantly more
stringent 95% mass limit of 4on excited quarks.

Figure 6: The observed 95% upper limits on the product of the cross section,
branching fraction, and acceptance for dijet resonances decaying to a quark
and a gluon (points). The corresponding expected limits (short dashed) and
their variations at the 1- and 2-standard deviation levels (shaded bands) are
also shown. Limits are compared to predictions for single $\PQb^{\ast}$
production (blue, dot dashed), the resonant component of the $\PQb^{\ast}$
production via contact interactions (magenta, long dashed), and the total
$\PQb^{\ast}$ signal from the sum of these two production modes (red, solid).

##  0.8 Summary

A search for heavy resonances decaying into quarks has been presented and no
excess has been found over the standard model (SM) expectations. The data were
collected by the CMS experiment at $\sqrt{s}=13\TeV$ during 2016–2018 and
correspond to an integrated luminosity of 138. Model-independent upper limits
are set on the product of the cross section of the resonance and its branching
fraction to quarks. Signals of bosons decaying to pairs of quarks are
considered, for both the previously explored sequential standard model (SSM),
and also for a new heavy vector triplet (HVT) model. The decays of bosons in
both the SSM and the HVT Model A are excluded at 95% confidence level for
masses from 1.8 to 2.4, and limits are set on the coupling strengths of the
HVT boson to SM bosons and fermions. Signals of an excited quark are
considered for a previously explored channel,
$\PQb\Pg\to\PQb^{\ast}\to\PQb\Pg$, and a production mode via contact
interactions. The excited quark is excluded at 95% confidence level for masses
from 1.8 to 4.0. This is the most stringent exclusion of excited quarks to
date.

## References

  * [1] CMS Collaboration, “Precise determination of the mass of the Higgs boson and tests of compatibility of its couplings with the standard model predictions using proton collisions at 7 and 8 TeV”,  Eur. Phys. J. C 75 (2015) 212, 10.1140/epjc/s10052-015-3351-7, arXiv:1412.8662.
  * [2] ATLAS Collaboration, “Measurement of the Higgs boson mass from the $H\rightarrow{}\gamma{}\gamma{}$ and $H\rightarrow{}Z{Z}^{*}\rightarrow{}4\ell{}$ channels in $pp$ collisions at center-of-mass energies of 7 and 8 TeV with the ATLAS detector”,  Phys. Rev. D 90 (2014) 052004, 10.1103/PhysRevD.90.052004, arXiv:1406.3827.
  * [3] ATLAS Collaboration, “Evidence for the spin-0 nature of the Higgs boson using ATLAS data”,  Phys. Lett. B 726 (2013) 120, 10.1016/j.physletb.2013.08.026, arXiv:1307.1432.
  * [4] CMS and ATLAS Collaborations, “Combined measurement of the Higgs boson mass in $pp$ collisions at $\sqrt{s}=7$ and 8 TeV with the ATLAS and CMS experiments”,  Phys. Rev. Lett. 114 (2015) 191803, 10.1103/PhysRevLett.114.191803, arXiv:1503.07589.
  * [5] V. D. Barger, W.-Y. Keung, and E. Ma, “A gauge model with light $W$ and $Z$ bosons”,  Phys. Rev. D 22 (1980) 727, 10.1103/PhysRevD.22.727.
  * [6] E. Salvioni, G. Villadoro, and F. Zwirner, “Minimal Z’ models: present bounds and early LHC reach”,  JHEP 09 (2009) 068, 10.1088/1126-6708/2009/11/068, arXiv:0909.1320.
  * [7] P. Langacker, “The physics of heavy Z’ gauge bosons”,  Rev. Mod. Phys. 81 (2009) 1199, 10.1103/revmodphys.81.1199, arXiv:0801.1345.
  * [8] D. Pappadopulo, A. Thamm, R. Torre, and A. Wulzer, “Heavy vector triplets: bridging theory and data”,  JHEP 09 (2014) 60, 10.1007/JHEP09(2014)060, arXiv:1402.4431.
  * [9] H. Terazawa, M. Yasuè, K. Akama, and M. Hayashi, “Observable effects of the possible sub-structure of leptons and quarks”,  Phys. Lett. B 112 (1982) 387, 10.1016/0370-2693(82)91075-9.
  * [10] E. Eichten and M. Peskin, “New tests for quark and lepton substructure”,  Phys. Rev. Lett. 50 (1983) 811, 10.1103/PhysRevLett.50.811.
  * [11] H. Harari, “Composite models for quarks and leptons”,  Phys. Rep. 104 (1984) 159, 10.1016/0370-1573(84)90207-2.
  * [12] U. Baur, I. Hinchliffe, and D. Zeppenfeld, “Excited Quark Production at Hadron Colliders”,  Int. J. Mod. Phys. A 2 (1987) 1285, 10.1142/S0217751X87000661.
  * [13] U. Baur, M. Spira, and P. M. Zerwas, “Excited quark and lepton production at hadron colliders”,  Phys. Rev. D 42 (1990) 815, 10.1103/PhysRevD.42.815.
  * [14] R. M. Harris and K. Kousouris, “Searches for Dijet Resonances at Hadron Colliders”,  Int. J. Mod. Phys. A 26 (2011) 5005, 10.1142/S0217751X11054905, arXiv:1110.5302.
  * [15] CMS Collaboration, “Search for narrow resonances and quantum black holes in inclusive and b-tagged dijet mass spectra from pp collisions at $\sqrt{s}$=7 TeV”,  JHEP 01 (2013) 013, 10.1007/jhep01(2013)013, arXiv:1210.2387.
  * [16] CMS Collaboration, “Search for resonances and quantum black holes using dijet mass spectra in proton-proton collisions at $\sqrt{s}=8\text{ }\text{ }\mathrm{TeV}$”,  Phys. Rev. D 91 (2015) 052009, 10.1103/PhysRevD.91.052009, arXiv:1501.04198.
  * [17] CMS Collaboration, “Search for narrow resonances in the b-tagged dijet mass spectrum in proton-proton collisions at $\sqrt{s}$=8 TeV”,  Phys. Rev. Lett. 120 (2018) 201801, 10.1103/PhysRevLett.120.201801, arXiv:1802.06149.
  * [18] CMS Collaboration, “Search for high mass dijet resonances with a new background prediction method in proton-proton collisions at $\sqrt{s}=$ 13 TeV”,  JHEP 05 (2020) 033, 10.1007/JHEP05(2020)033, arXiv:1911.03947.
  * [19] ATLAS Collaboration, “Search for new resonances in mass distributions of jet pairs using 139 fb-1 of $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector”,  JHEP 03 (2020) 145, 10.1007/JHEP03(2020)145, arXiv:1910.08447.
  * [20] G. Ciezarek et al., “A Challenge to Lepton Universality in B Meson Decays”,  Nature 546 (2017) 227, 10.1038/nature22346, arXiv:1703.01766.
  * [21] LHCb Collaboration, “Test of lepton universality in beauty-quark decays”,  Nature Phys. 18 (2022) 277, 10.1038/s41567-021-01478-8, arXiv:2103.11769.
  * [22] CMS Collaboration, “EXO-20-008”. HEPData record for this analysis, 2022.  https://doi.org/10.17182/hepdata.127768.
  * [23] CMS Collaboration, “The CMS experiment at the CERN LHC”,  JINST 3 (2008) S08004, 10.1088/1748-0221/3/08/S08004.
  * [24] CMS Collaboration, “Description and performance of track and primary-vertex reconstruction with the CMS tracker”,  JINST 9 (2014) P10009, 10.1088/1748-0221/9/10/P10009, arXiv:1405.6569.
  * [25] CMS Collaboration, “Performance of the CMS Level-1 trigger in proton-proton collisions at $\sqrt{s}=13$ TeV”,  JINST 15 (2020) P10017, 10.1088/1748-0221/15/10/P10017, arXiv:2006.10165.
  * [26] CMS Collaboration, “The CMS trigger system”,  JINST 12 (2017) P01020, 10.1088/1748-0221/12/01/P01020, arXiv:1609.02366.
  * [27] J. Alwall et al., “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”,  JHEP 07 (2014) 079, 10.1007/JHEP07(2014)079, arXiv:1405.0301.
  * [28] T. Sjöstrand et al., “An introduction to PYTHIA 8.2”,  Comput. Phys. Commun. 191 (2015) 159, 10.1016/j.cpc.2015.01.024, arXiv:1410.3012.
  * [29] P. Skands, S. Carrazza, and J. Rojo, “Tuning PYTHIA 8.1: the Monash 2013 Tune”,  Eur. Phys. J. C 74 (2014) 3024, 10.1140/epjc/s10052-014-3024-y, arXiv:1404.5630.
  * [30] CMS Collaboration, “Event generator tunes obtained from underlying event and multiparton scattering measurements”,  Eur. Phys. J. C 76 (2016) 155, 10.1140/epjc/s10052-016-3988-x, arXiv:1512.00815.
  * [31] R. D. Ball et al., “Parton distributions with QED corrections”,  Nucl. Phys. B 877 (2013) 290, 10.1016/j.nuclphysb.2013.10.010, arXiv:1308.0598.
  * [32] R. D. Ball et al., “Unbiased global determination of parton distributions and their uncertainties at NNLO and at LO”,  Nucl. Phys. B 855 (2012) 153, 10.1016/j.nuclphysb.2011.09.024, arXiv:1107.2652.
  * [33] GEANT4 Collaboration, “—a simulation toolkit”,  Nucl. Instrum. Meth. A 506 (2003) 250, 10.1016/S0168-9002(03)01368-8.
  * [34] CMS Collaboration, “Particle-flow reconstruction and global event description with the CMS detector”,  JINST 12 (2017) P10003, 10.1088/1748-0221/12/10/P10003, arXiv:1706.04965.
  * [35] CMS Collaboration, “Pileup mitigation at CMS in 13 TeV data”,  JINST 15 (2020) P09018, 10.1088/1748-0221/15/09/p09018, arXiv:2003.00503.
  * [36] M. Cacciari, G. P. Salam, and G. Soyez, “The anti-jet clustering algorithm”,  JHEP 04 (2008) 063, 10.1088/1126-6708/2008/04/063, arXiv:0802.1189.
  * [37] M. Cacciari, G. P. Salam, and G. Soyez, “FastJet user manual”,  Eur. Phys. J. C 72 (2012) 1896, 10.1140/epjc/s10052-012-1896-2, arXiv:1111.6097.
  * [38] M. Cacciari, G. P. Salam, and G. Soyez, “The catchment area of jets”,  JHEP 04 (2008) 005, 10.1088/1126-6708/2008/04/005, arXiv:0802.1188.
  * [39] CMS Collaboration, “Jet energy scale and resolution in the CMS experiment in pp collisions at 8 TeV”,  JINST 12 (2017) P02014, 10.1088/1748-0221/12/02/P02014, arXiv:1607.03663.
  * [40] CMS Collaboration, “Jet algorithms performance in 13 TeV data”, CMS Physics Analysis Summary CMS-PAS-JME-16-003, 2017.
  * [41] E. Bols et al., “Jet Flavour Classification Using DeepJet”,  JINST 15 (2020) P12012, 10.1088/1748-0221/15/12/P12012, arXiv:2008.10519.
  * [42] CMS Collaboration, “Performance of the DeepJet b tagging algorithm using 41.9/fb of data from proton-proton collisions at 13 TeV with Phase 1 CMS detector”, CMS Detector Performance Note CMS-DP-2018-058, 2018.
  * [43] M. H. Kutner, C. J. Nachtsheim, J. Neter, and W. Li, “Applied linear statistical models”.  Irwin, 1996.  Fifth edition. ISBN 9780256117363.
  * [44] M. J. Oreglia, “A study of the reactions $\psi^{\prime}\to\gamma\gamma\psi$”.  PhD thesis, Stanford University, 1980.  SLAC Report SLAC-R-236.
  * [45] CMS Collaboration, “Identification of heavy-flavour jets with the CMS detector in pp collisions at 13 TeV”,  JINST 13 (2018) P05011, 10.1088/1748-0221/13/05/P05011, arXiv:1712.07158.
  * [46] CMS Collaboration, “Precision luminosity measurement in proton-proton collisions at $\sqrt{s}=$ 13 TeV in 2015 and 2016 at CMS”,  Eur. Phys. J. C 81 (2021) 800, 10.1140/epjc/s10052-021-09538-2, arXiv:2104.01927.
  * [47] CMS Collaboration, “CMS luminosity measurement for the 2017 data-taking period at $\sqrt{s}=13\TeV$”, CMS Physics Analysis Summary CMS-PAS-LUM-17-004, 2018.
  * [48] CMS Collaboration, “CMS luminosity measurement for the 2018 data-taking period at $\sqrt{s}=13\TeV$”, CMS Physics Analysis Summary CMS-PAS-LUM-18-002, 2019.
  * [49] J. Butterworth et al., “PDF4LHC recommendations for LHC Run II”,  J. Phys. G 43 (2016) 23001, 10.1088/0954-3899/43/2/023001, arXiv:1510.03865.
  * [50] CMS and ATLAS Collaborations, “Procedure for the LHC Higgs boson search combination in Summer 2011”, CMS Note CMS-NOTE-2011-005, ATL-PHYS-PUB-2011-11, 2011.
  * [51] CMS Collaboration, “Combination of CMS searches for heavy resonances decaying to pairs of bosons or leptons”,  Phys. Lett. B 798 (2019) 134952, 10.1016/j.physletb.2019.134952, arXiv:1906.00057.
  * [52] T. Junk, “Confidence level computation for combining searches with small statistics”,  Nucl. Instrum. Meth. A 434 (1999) 435, 10.1016/S0168-9002(99)00498-2, arXiv:hep-ex/9902006.
  * [53] A. L. Read, “Presentation of search results: the CLs technique”,  J. Phys. G 28 (2002) 2693, 10.1088/0954-3899/28/10/313.
  * [54] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, “Asymptotic formulae for likelihood-based tests of new physics”,  Eur. Phys. J. C 71 (2011) 1554, 10.1140/epjc/s10052-011-1554-0, arXiv:1007.1727. [Erratum: 10.1140/epjc/s10052-013-2501-z].


††thanks<EMAIL_ADDRESS>
# Different relative scalings between transient forces and thermal
fluctuations
tune regimes of chromatin organization

Anna Colettia, Benjamin L. Walker b, Kerry Bloom c, and Katherine A. Newhalla,
aDepartment of Mathematics, University of North Carolina at Chapel Hill,
Chapel Hill, NC, 27510  bDepartment of Mathematics, University of California,
Irvine, Irvine, CA, 92697  cDepartment of Biology, University of North
Carolina at Chapel Hill, Chapel Hill, NC, 27510

###### Abstract

Within the nucleus, structural maintenance of chromosome protein complexes,
namely condensin and cohesin, create an architecture to facilitate the
organization and proper function of the genome. Condensin, in addition to
performing loop extrusion, creates localized clusters of chromatin in the
nucleolus through transient crosslinks. Large-scale simulations revealed three
different dynamic behaviors as a function of timescale: slow crosslinking
leads to no clusters, fast crosslinking produces rigid slowly changing
clusters, while intermediate timescales are optimal for producing flexible
clusters that mediate gene interaction. By mathematically analyzing different
relative scalings of the two sources of stochasticity, thermal fluctuations
and the force induced by the transient crosslinks, we predict these three
distinct regimes of cluster behavior. Standard time-averaging that takes the
fluctuations of the transient crosslink force to zero can predict the
existence of clusters, but not their timescale-dependent lifetimes. Accounting
for the interaction of both fluctuations from the crosslinks and thermal noise
with an effective energy landscape does capture the timescale-dependent
flexible cluster lifetimes. No clusters are predicted when the fluctuations of
the transient crosslink force are taken to be large relative to thermal
fluctuations. This mathematical perturbation analysis illuminates the
importance of accounting for stochasticity in local incoherent transient
forces to predict emergent complex biological behavior.

##  I Introduction

The genome has an intricate and hierarchical organization that allows cells to
fit an extraordinary amount of genetic material and perform important nuclear
functions within a nucleus mere micrometers in size. Within the nucleus, DNA
is complexed with a variety of proteins to form the fibrous material chromatin
that makes up chromosomes. Contributing to chromosomal architecture are
structural maintenance of chromosome (SMC) protein complexes, namely condensin
and cohesin, which facilitate chromosome segregation and gene regulation [18].
Such processes aid in the compaction of genetic material within the nucleus as
well as allow the cell to locate and access specific genes depending on cell
cycle, cell type, or environmental cue [20]. In particular, condensin helps
regulate the retrieval of genetic information by acting on chromatin to form
stochastic gene-gene crosslinks [6] and generate loops to induce gene mixing
[9, 11, 21, 10].

The emergent behavior induced by condensin crosslinks has been studied
computationally with polymer models of chromatin [8, 23, 24, 12, 26, 27] in
order to more directly observe the dynamics of DNA over experimental
techniques. By representing 5kbp of DNA as beads, linked through worm-like-
chain springs to form a long polymer chromatin chain, the timescale of the
random SMC crosslink binding between different beads was shown to control the
formation of clusters and gene interactions [12, 26]. These clusters, or gene
neighborhoods, reveal an underlying organizational framework that optimizes
gene interactions, retains a high level of genomic compaction, and yet is
highly flexible. A similar detailed level of information is not yet available
experimentally; high throughput population methods such as Chromosome
Conformation Capture (Hi-C) methods lack temporal evolution information, while
microscopic imaging lacks genome-level resolution [22]. The ability to access
both the spatial arrangement and the temporal re-arrangement of beads, which
is possible in simulation, is crucial for advancing our understanding of life.

Building and remodeling these high-level genetic neighborhoods has the
potential to reveal the underlying mechanics for the configuration of the
energy landscape within our nuclei. In this paper, we reveal the mathematical
mechanism that describes how fluctuations in the model lead to the formation
of bead clusters and the temporal mixing of beads between clusters. We vary
the timescale of crosslinking forces between beads to change the size of this
fluctuating force relative to thermal fluctuations. These various timescales
produce three different dynamic regimes. At fast bead crosslinking timescales,
rigid, unchanging clusters form. At slow crosslinking timescales, the beads
interact in an amorphic state with no clusters. At optimal mean crosslink
lifetimes, clusters can both form and exchange beads, indicative of the
required gene interaction for proper function. Tuning the timescales of
crosslinking provides a unifying mechanism of genome organization that gives
insight into how a large number of configurational states can be rapidly
remodeled as cells encounter biological challenges.

Configurations that can be rapidly remodeled were shown to take advantage of
the interaction between the crosslink binding force fluctuations and the
thermal noise fluctuations to promote faster bead exchange between clusters
[27]. We build upon this work by incorporating the timescale of the crosslink
binding into the mathematical analysis to show how different relative scalings
of the two sources of stochastic fluctuations, thermal noise and the force
induced by the transient bonds, predict the three distinct regimes of cluster
behavior mentioned above. Standard time-averaging that takes the fluctuations
of the transient binding force to zero (fast crosslinking timescale) predicts
the existence of long-lived rigid clusters. No clusters are predicted when the
fluctuations of the transient binding force are taken to be large (slow
crosslinking timescale) relative to the thermal fluctuations. Only when
accounting for the interaction of both fluctuations from the binding and
thermal noise with an effective energy landscape do we recover the flexible
clusters and their timescale-dependent mixing.

This work provides the mathematical mechanism underlying the timescale-
dependent effects of active agents in biological systems. The timescale
controls the production of fluctuations in the forces that the active agents
produce. The relative size of these fluctuations to the thermal fluctuations
produces different emergent temporal behavior. This mechanism is not
restricted to just the chromatin dynamics studied here but is potentially
applicable to a wide range of cell-level biological processes such as homology
searches for genetic recombination [4, 7, 2], defensive mucus barriers
emerging from polymer chains of mucin [28, 17, 25], and error correction
through configurational state changes [3, 13, 16]. Thus to deepen our
understanding of how such biological systems function, we have shown it is
necessary to properly analyze the interplay of different fluctuations and the
timescale on which they are generated.

##  II Results

###  II.1 Timescale of crosslink binding force drives cluster formation and
dynamics

We establish that it is the timescale of the stochastically-switching SMC
protein crosslinking force through random binding and unbinding that affects
the organization of the beads and the dynamics of the clusters. To show this,
we create an idealized 361-bead model that represents the nucleolus, with each
bead representing 5 kilobase pairs of DNA in the chromatin chain. The beads
bind and unbind based on a biologically feasible random model for the SMC
protein crosslinks that is more likely to choose pairs of nearby beads to
bind.

Fig. 1 shows the formation of clusters for fast-enough timescale, $\beta$, of
the random crosslink binding force. The repulsive excluded volume force and
thermal fluctuations acting on the beads are in competition with the
attractive stochastically changing pairwise binding forces that coalesce many
beads into a cluster. Furthermore, changing the timescale of the crosslink
binding force, $\beta>0$, from slow to fast results in three distinct
clustering behaviors: The system passes from an amorphic (dissolution of
clusters with many bead interactions) to a flexible (frequent bead exchanges
between clusters) and then rigid clustering (minimal bead interaction between
clusters) behavior.

The observed timescale-dependent clustering behavior mirrors that seen in [26]
which employed a larger polymer-like chromosome model and a different random
model for the crosslinking proteins in the nucleolus. This demonstrates the
universality of the clustering dynamics that are driven by the timescale of
the binding forces.

Figure 1: (A-C) Snapshot of 361 bead model for increasing binding timescale.
We observe three clustering regimes: (A) amorphic, (B) flexible, and (C)
rigid. (D) Mixing coefficient and average nearby number of beads for a range
of $\beta$-values for 361 bead model. Flexible clustering occurs for the range
of $\beta$-values within the grey region. In this regime, clusters emerge and
beads interact.

###  II.2 Time-averaged force predicts existence of clusters but not their
dynamics

To mathematically determine the effective attractive force generated by the
stochastic binding we work with a further reduced 3-bead model. We start by
showing that a simple time-averaging of the stochastic binding force can
predict the existence of clusters but not their timescale-dependent mixing
dynamics.

The chosen random model for the crosslinking forces has the added advantage of
fitting into a continuous-time Markov chain (CTMC) framework. For the
idealized 3-bead model, the four Markov states are easily enumerated: all
beads unbound ($s=1$), beads 1 and 2 bound ($s=2$), beads 1 and 3 bound
($s=3$), and beads 2 and 3 bound ($s=4$). The time evolution of the states
follows a general CTMC process with a transition rate matrix $S$. Included in
$S$ is the bead binding-rate that is proportional to the bead separation
distance given by affinity function $a(\cdot)$ and a constant bond breaking-
rate $c$.

Now that the binding stochasticity is formulated as a Markov chain, finding
the time-averaged force on each bead is straightforward [19]. The percent of
time spent in each Markov chain state is given by $\vec{r}(\vec{x})$, the
normalized null vector of the switching matrix $S(\vec{x})$, so that the time-
averaged force for each coordinate is

$\frac{dx_{k}}{dt}=\sum_{s=1}^{4}v_{k}(\vec{x};s)r_{s}(\vec{x})=\langle
v_{k}(\vec{x})\rangle\;\text{for}\;k=\\{1,\dots,6\\}$ (1)

where $v_{k}(\vec{x};s)$ is the force on coordinate $k$ when the system is in
Markov state $s$. (Note, we have concatenated the $x$ and $y$ coordinates of
each bead into a single vector $\vec{x}$.) This removes the stochasticity and
results in a deterministic effective force for the system,
$\langle\vec{v}\rangle$.

The stable fixed points of (1) are the observed clusters. One fixed point is a
3-bead cluster in which the beads are in a small triangle configuration; see
Fig. 2C. This is distinguished from the unbound state $s=1$, in which the
beads form a triangle but with larger pairwise distances; see Fig. 2A. The
other three fixed points are 2-bead clusters in which two beads are
superimposed and the third is farther away and unbound; see Fig. 2B. These are
distinguished from states $s=2,3,4$ as the clusters persist longer than any
single bond lifetime.

In addition to predicting the existence of clusters, we can test this time-
averaging procedure’s ability to predict the lifetime of a cluster. We do this
by considering the effects of small perturbations of noise about the average
by studying the equation,

$dx_{k}=\langle v_{k}(\vec{x})\rangle dt+\sqrt{2\epsilon}\;dB_{k}$ (2)

in the limit as $\epsilon\to 0$. Even though the vector of forces is the
gradient of a potential function for each state $s$, i.e.
$\vec{v}(\vec{x};s)=-\nabla_{x}U(\vec{x};s)$, it is not true that
$\langle\vec{v}(\vec{x})\rangle=-\nabla_{x}\langle U(\vec{x})\rangle$ since
the CTMC null vector $\vec{r}$ depends on the positions $\vec{x}$. However, we
numerically find an effective potential $U_{\text{eff}}$ such that
$-\nabla_{x}U_{\text{eff}}(\vec{x})\approx\langle\vec{v}(\vec{x})\rangle$
along the most probable transition path connecting a minimum of
$U_{\text{eff}}$ (recall these are the stable fixed points of (1)) to a saddle
point using the String Method [5]. We are therefore justified in
asymptotically approximating the escape times from a cluster using the well-
known Arrhenius law,

$\log(\mathbb{E}[\tau])\sim\frac{\Delta U_{\text{eff}}}{\epsilon}$ (3)

in the limit as $\epsilon\to 0$, where $\Delta U_{\text{eff}}$ is the change
in the effective potential generated by the time-averaged force along the most
probable transition path [29].

Figure 2: (A-C) plots of pairwise distances of three beads over time and a
snapshot of the three bead configuration: (A) all beads unbound: balance
between excluded volume and confinement forces keep beads in a large triangle
configuration, (B) 2-bead cluster state with beads 1-2 bound, (C) 3-bead
cluster state: all beads remain close, rapid switching between pairwise bonds
keeps beads in a small triangle configuration.

Despite predicting the existence of clusters, the naive time-averaged force
does not produce an effective potential $U_{\text{eff}}$ that accounts for the
observed differences in the cluster lifetime as the crosslink binding
timescale $\beta$ varies in the 361-bead model. The time-averaged force, and
thus its effective potential, is invariant to the overall timescale of the
transition rate matrix. Scaling $S$ by some constant $\alpha$ does not affect
the fraction of time spent in each state because $\vec{r}$ is also in the null
space of $\alpha S$. Thus, Eqs. (1) and (2) are unchanged when scaling the
transition rate matrix $S$ by $\alpha$ and therefore cannot predict subsequent
changes to the dynamic cluster lifetimes.

###  II.3 Relative scalings of fluctuations predicts three different
clustering regimes

The key to understanding clustering dynamics is accounting for fluctuations
from both the stochastic binding and the small perturbations of thermal noise.
To analyze the dynamics mathematically we perform an asymptotic expansion as
the size of the fluctuations goes to zero. We consider taking both
fluctuations to zero at the same rate and at different rates. The latter
allows one noise source to dominate. Rigid clusters arise when the binding
fluctuations are small and only thermal noise drives transitions between
clusters. Amorphic arrangements arise when the binding fluctuations are large
and take place on a longer timescale than the fast thermal fluctuations. This
binding drives transitions between clusters. Taking fluctuations to zero at
the same rate allows for interaction between the noises. We find that this
interaction is required to predict the flexible cluster dynamics that depend
on the timescale of the binding. Together, these three asymptotic regimes
explain the three regimes seen in [26] and Fig. 1.

We include the effects of fluctuations about the mean force so that the
position of each bead follows the overdamped Langevin-like equation given by

$dx_{k}=v_{k}(\vec{x};s)dt+\sqrt{2\epsilon}dB_{k}.$ (4)

where $s$ is the state of the CTMC. How we choose to scale the transition rate
matrix $S$ determines how the relative rates of the two different fluctuations
go to zero.

We start with the scaling $\frac{1}{\epsilon^{2}}S$ which takes the binding
fluctuations to zero faster than the perturbations from thermal noise. We seek
an effective potential $W$ to replace $U_{\text{eff}}$ in Eq. (3). Therefore
we assume a steady-state distribution

$p_{s}(\vec{x})=r_{s}(\vec{x})\exp\bigg{(}-\frac{1}{\epsilon}W(\vec{x})\bigg{)}$
(5)

as the asymptotic solution for the system of steady-state Fokker-Planck
equations,

$\frac{\partial p_{s}}{\partial t}=-\sum_{k=1}^{6}\frac{\partial}{\partial
x_{k}}[v_{k}(\vec{x};s)p_{s}]+\epsilon\sum_{k=1}^{6}\frac{\partial^{2}}{\partial
x_{k}^{2}}[p_{s}]+\frac{\alpha}{\epsilon^{2}}\sum_{j=1}^{4}S_{sj}p_{j}.$ (6)

The element of the CTMC transition matrix $S_{sj}$ gives the transition rates
from state $j$ into state $s$. Thus, the last term couples the process between
different states. Setting (6) to zero and plugging in (5) yields the leading
order equation

$S\vec{r}=0$ (7)

with next order equation

$r_{s}\vec{v}\cdot\nabla_{x}W+r_{s}\nabla_{x}W\cdot\nabla_{x}W=0.$ (8)

Thus, $\vec{r}$ is the steady-state distribution of $S$ and summing Eq. (8)
over $s$ reveals that $-\nabla_{x}W=\sum_{s}\vec{v}(\vec{x};s)r_{s}(\vec{x})$
is exactly the naive time-average from above that does not depend on $\alpha$.

If instead we take the scaling $\frac{1}{\epsilon}S$, both sources of noise
show up in the leading order equation given by

$[A(\vec{x},\nabla_{x}W)+D(\nabla_{x}W)+\alpha S(\vec{x})]\vec{r}=0$ (9)

for advection (drift) matrix $A$, diffusion matrix $D$, and switching rate
matrix $S$; see [27] for details. In this regime, both sources of noise are
important for flexible clustering. We see the effects of varying the binding
timescale parameter $\alpha$, which controls how quickly the simulation
switches between the Markov chain states. This role mirrors that of the
kinetic timescale parameter $\mu$ in the large-scale simulations of [12, 26]
in which $\mu$ is a “tuning knob” used to set the kinetic timescale on which
the crosslinks bind and unbind.

To further illustrate the changing dynamics with $\alpha$, we compare the
predicted lifetime of the cluster governed by Eq. (3) with Monte-Carlo
simulations. The asymptotic escape times from one 2-bead cluster to another
2-bead cluster are shown in Fig. 3. The slope value is given by the
quasipotential barrier height and well expresses the linear relationship of
the mean escape times. As $\alpha$ increases, the effective energy barrier
approaches the barrier predicted by the naive time average discussed
previously.

Figure 3: Transition: 2-bead cluster to 2-bead cluster. (a) Quasipotential
(solid) for different $\alpha$ values and for the time-averaged binding force
(dotted) along the transition path, normalized to unit length. (b) Comparison
of asymptotic escape times computed via Monte Carlo simulation to slope taken
from quasipotential barrier height. (c) Quasipotential energy barrier height
vs. switching timescale $\alpha$.

Finally, if we keep $S$ independent of $\epsilon$, retaining all the
fluctuations, the system equilibrates within each Markov chain state and there
are no clusters predicted. In this regime, the effects of the binding noise
are pushed into a higher order so that the leading order does not contain any
contributions from the CTMC. We observe the beads are often unbound, settling
into a “large triangle” configuration; see Fig. 2A. This long-observed state
of the system is a balance between the excluded volume and the confinement
forces. We do not consider this a cluster as it would dissolve in the absence
of the confinement forces. Fig. 4 shows the formation of clusters as $\alpha$
increases. The system is able to remain in a cluster state despite all the
beads being unbound as the system does not have time between binding events to
reach the large triangle equilibrium.

Figure 4: Average proportion of simulation spent in $s=1$ unbound state of the
Markov chain and average proportion of simulation not spent in a cluster state
vs. $\alpha$.

###  II.4 Effects of crosslink binding fluctuations on different transitions

Scaling $S$ by $\frac{1}{\epsilon}$ revealed the dependence of cluster
lifetimes on the crosslink binding timescale, $\alpha$. By accounting for both
sources of randomness in the system we computed an accurate quasipotential
that predicted a higher energy barrier height between minimum states as
$\alpha$ increased, and thus longer cluster lifetimes. However, this increase
is not uniform across all transitions as we show in Fig. 5.

Of the ‘uphill’ most probable transition paths the system takes from one
minimum state to the saddle point shown in Fig. 5A-D, only the 3-bead on its
way to a 2-bead cluster has noticeable dependence on $\alpha$ (Fig. 5C). This
transition also has the largest change in quasipotential energy barrier with
$\alpha$ (Fig. 5F) indicating removing a bead from a cluster is more sensitive
to the binding timescale that rearrangements within the cluster (pathway in
Fig. 5D). The transition shown in Fig. 5D to the “collinear” saddle point
configuration also has the lowest energy barrier height amongst the
transitions, with limited dependence on $\alpha$ for either the pathway or the
height. This indicates that the binding force is not important for this
transition, as the stochasticity of the binding does not help the system make
the transition.

Starting in the 2-bead cluster state, our system can transition into a
different 2-bead cluster (Fig. 5A) or into the “small triangle” 3-bead cluster
(Fig. 5B). These transitions have similar pathways up to each saddle point
that do not significantly depend on $\alpha$, and have similar quasipotential
energy barrier heights for the range of $\alpha$ values from 0 to 20. This
suggests that the system will take either transition with similar probability
and do so in approximately the same amount of time.

Figure 5: (A-D) Transition pathways to a saddle point: (A) 2-bead cluster on
the way to a different 2-bead cluster (B) 2-bead cluster on the way to a
3-bead cluster, (C) 3-bead cluster on the way to 2-bead cluster for
$\alpha=1$, $\alpha=10$, and time-averaged binding force ($\alpha\to\infty$),
(D) 3-bead cluster to the collinear saddle point configuration that continues
to a rearrangement of the 3-bead cluster, or to a different saddle point on
the way to a 2-bead cluster. (E,F) Energy Barrier heights vs. $\alpha$ for
these transition pathways

##  III Discussion

The underlying mathematical mechanism behind the observed chromatin clustering
behaviors is the competition between thermal fluctuations and the timescale of
the model condensin crosslinking force. This timescale controls the production
of random fluctuations in the crosslink binding force and in turn the
dominating terms in the asymptotic perturbation when searching for an
effective thermal equilibrium. This analysis emphasizes the importance of
accounting for stochasticity in local incoherent transient forces to predict
emergent complex biological behavior.

The key to understanding flexible clustering dynamics is accounting for
fluctuations from both the stochastic binding force and the thermal noise,
taking both fluctuations to zero at the same rate. To produce this interaction
mathematically, we scale the binding rate matrix $S$ by $\frac{1}{\epsilon}$.
Taking $\epsilon\to 0$ with this scaling, an effective binding force remains
that allows for stable clusters. This limit also accounts for the fact that
the fluctuations of the binding force aid in transitions between clusters, as
it creates finite periods of time when the thermal noise has to overcome
smaller forces to push the system into a new cluster state. Since our added
binding timescale parameter $\alpha$ controls the size of the fluctuations, it
too controls the size of the effective energy barrier for predicting
transition times between stable cluster states given by the Arrhenius law.

If we take the fluctuations to zero at different rates, one source of noise
“dominates” over the other in the limit as $\epsilon\to 0$. Choosing to scale
$S$ by $\frac{1}{\epsilon^{2}}$, the fluctuations of the binding force go to
zero faster than the thermal fluctuations. An effective binding force remains
that allows for stable clusters, but the fluctuations in the binding force are
now much smaller and do not interact with the thermal noise to aid in
transitions between clusters. Indeed, we recover the naive time average to
describe the effective energy barrier, which is also the limiting barrier as
$\alpha\to\infty$ of the above-mentioned fluctuation-interaction case.

If we flip which fluctuation source dominates in the limit as $\epsilon\to 0$
by keeping $S$ independent of $\epsilon$, switching between Markov chain
states occurs so infrequently that the system equilibriates while in this one
state. No effective force is generated by the superposition of different
states, thus no clusters are predicted. In the 3-bead model, two cluster-like
states appear but their generating mechanism would not produce clusters in the
361-bead model. The first is a 3-bead triangle configuration, but it is a
result of the balance between the repulsive excluded volume force and
confinement force. In the 361-bead model, this would correspond to all beads
in an amorphic arrangement. The second is a 2-bead-like cluster formed by the
fact that a crosslink binds two beads; this cluster returns to the 3-bead
triangle configuration once the bond is broken. In the 361-bead model, this
would correspond to another amorphic arrangement with many pair-wise bound
beads.

The arrangement of beads within clusters provides a physical model for the
dynamics of gene clusters within the nucleus. Gene clustering is a mechanism
that facilitates the co-regulation of un-linked genes or of long linear arrays
of repeated genes such as rDNA. The ability to reconfigure the network is
central to the cell’s ability to re-wire its transcriptional circuitry. Using
a relatively small number of beads to reveal root mechanisms, we account for
both state changes and dynamics between states with thermal fluctuations and
cross-linking forces. The energy barrier for rearrangements within a cluster
is lower than that required for exchanges between clusters (Fig. 5E, F). The
system is remarkably plastic, with small bursts of energy able to potentiate
new configurations depending on the biological demand.

The emergence of a wide range of clustering regimes (amorphic, flexible,
rigid)(Fig. 1) and exchange within and between clusters with a pared-down
3-bead model (Fig. 5) provides a framework for understanding the governing
principles of genome organization: The continuous chain of each chromosome
biases their individualization. The crosslinking forces provide a mechanism to
build gene (bead) clusters and generate informational circuits that are
agnostic to the position of a gene in a given chromosome. This work shows that
by controlling the timescale of the crosslinking forces, the plasticity of the
gene clusters can be tuned based on biological needs. The predictions of
cluster dynamics for larger systems will depend greatly on thermodynamics and
tuning the crosslinking forces. This model and analysis provide critical
insight into the diversity of network configurations with a minimal set of
parameters.

Computationally, our analysis is limited to relatively small numbers of beads
due to the need to enumerate all possible pairwise bound states for the Markov
chain. With the growing power of machine learning, it is likely possible to
build a physics-informed neural network [14, 15] to learn the effective
potential. This would allow predictions of cluster dynamics for larger systems
but would need to be re-applied for each timescale or set of model parameters.
The effective potential alone gives little insight into the mechanism
underlying the emergent behavior that our analysis of the 3-bead model has
provided. This new outlook on the importance of noise at the proper timescale
aids in deepening our understanding of life at the cellular level.

##  IV Methods

The idealized chromatin model is based on the polymer-like chain of beads
model in [12, 26]. The position of bead $i$ is given by the overdamped
Langevin equation

$d\mathbf{x}_{i}=(\mathbf{f}_{c}^{i}+\mathbf{f}_{\text{EV}}^{i}+\mathbf{f}_{\text{bond}}^{i})dt+\sqrt{2\epsilon}d\mathbf{B}.$
(10)

The parameter $\epsilon$ scales the vector of Brownian increments
$d\mathbf{B}$ preparing for asymptotic analysis to take the fluctuations to
zero. The deterministic confinement force,

$\mathbf{f}_{c}^{i}=-\eta\mathbf{x}_{i},$

replaces the hard-wall constraint of the nucleus membrane while the excluded
volume force,

$\mathbf{f}_{\text{EV}}^{i}=\sum_{j\neq
i}a_{ev}(\mathbf{x}_{i}-\mathbf{x}_{j})\exp\bigg{(}-\frac{|\mathbf{x}_{i}-\mathbf{x}_{j}|^{2}}{c_{ev}}\bigg{)},$

remains similar to its form in [12, 26]. Parameters $\eta$, $a_{ev}$, and
$c_{ev}$ are given in Table 1 for both the 361-bead
$(\mathbf{x}_{i}\in\mathbb{R}^{3},i=1\dots 361)$ and the 3-bead
$(\mathbf{x}_{i}\in\mathbb{R}^{2},i=1\dots 3)$ models. Note we have neglected
the spring force linking the beads together to form a chain, as we show it is
secondary to the clustering dynamics.

The stochastic binding force,

$\mathbf{f}_{\text{bond}}^{i}=\sum_{j\neq i}\kappa
b_{ij}(\mathbf{x}_{j}-\mathbf{x}_{i}),$

models the binding SMC proteins found in the biological system. This force
binds two beads, corresponding to $b_{ij}=1$ if beads $i$ and $j$ are bound
and 0 otherwise. Each bead can be bound to only one other bead; these
stochastic bonds form and break at exponentially distributed times with a
binding-rate proportional to the bead separation distance,
$a(\mathbf{x}_{i}-\mathbf{x}_{j})$, and constant breaking-rate $c$. The time
evolution of the states follows a general CTMC process; for the 3-bead model
the transition rate matrix is given by

$S=\begin{pmatrix}b&c&c&c\\\ a(\mathbf{x}_{1}-\mathbf{x}_{2})&-c&0&0\\\
a(\mathbf{x}_{1}-\mathbf{x}_{3})&0&-c&0\\\
a(\mathbf{x}_{2}-\mathbf{x}_{3})&0&0&-c\end{pmatrix}$ (11)

with
$b=-a(\mathbf{x}_{1}-\mathbf{x}_{2})-a(\mathbf{x}_{1}-\mathbf{x}_{3})-a(\mathbf{x}_{2}-\mathbf{x}_{3})$.

The affinity function $a(\mathbf{x})$ and parameter values are given in Table
1 for both the 361-bead and the 3-bead models. Both $a(\mathbf{x})$ and $c$
are scaled by $\beta$ for the 361-bead model and $\alpha$ for the 3-bead model
to explore the kinetic timescale on which the crosslinks bind and unbind,
mirroring the kinetic timescale parameter $\mu$ in the large-scale simulations
of [26]. Both $a(\mathbf{x})$ and $c$ are further scaled by different powers
of $\epsilon$ to perform the asymptotic analysis of the fluctuations.

Table 1: Model parameters and formula. | 361-Bead Model | 3-Bead Model
---|---|---
$a(\mathbf{x})$ | $\displaystyle\frac{2}{1+\exp({20(|\mathbf{x}|-75)})}$ | $\displaystyle\frac{2}{1+\exp({20(|\mathbf{x}|-0.75))}}$
$\eta$ | 0.002 | 1
$a_{ev}$ | 0.0332 | 2
$c_{ev}$ | 30600 | 0.5
$\kappa$ | 10.9 | 5
$c$ | 0.01 | 0.5

### Monte Carlo Simulations

We perform Monte Carlo simulations of the 3-bead system to compare the scaling
of the escape times with $\epsilon$ to the effective energy barriers predicted
by the theoretical calculations. The simulations are started in either the
2-bead or 3-bead cluster, and continued until a stopping condition is met,
indicating that the system has left the basin of attraction of the initial
state. From the 2-bead cluster, we look for either a different 2-bead cluster
or the 3-bead triangle cluster. From the 3-bead triangle cluster, we look for
either a 2-bead cluster or a collinear configuration that is the saddle point
between different arrangements of the 3-bead triangle cluster. The mean escape
time, $\tau$, is computed by the maximum likelihood estimate that divides the
sum of all escape times by the number executing the desired transition.

Code for reproducing results is available on GitHub [1].

###### Acknowledgements.

AC and KN were partially supported by National Science Foundation (NSF) grants
DMS-1816394 and DMS-2307297 and AC was also partially supported by NSF grant
DMS-1929298. KB was partially supported by the National Institutes of Health
(NIH) grant R01 GM32238.

## References

  * [1] bwalker1/quasi-string-reprod: Quasipotential Strings v1.
  * [2] Jason C. Bell and Stephen C. Kowalczykowski.  RecA: Regulation and Mechanism of a Molecular Search Engine.  Trends in Biochemical Sciences, 41(6):491–507, June 2016.
  * [3] Kerry Bloom and Elaine Yeh.  Tension Management in the Kinetochore.  Current Biology, 20(23):R1040–R1048, December 2010.
  * [4] Alyson J. Conover, Claudia Danilowicz, Ruwan Gunaratne, Vincent W. Coljee, Nancy Kleckner, and Mara Prentiss.  Changes in the tension in dsDNA alter the conformation of RecA bound to dsDNA–RecA filaments.  Nucleic Acids Research, 39(20):8833–8843, November 2011.
  * [5] Weinan E, Weiqing Ren, and Eric Vanden-Eijnden.  String method for the study of rare events.  Physical Review B, 66(5):052301, August 2002.
  * [6] John F Marko Elnaz Alipour.  Self-organization of domain structures by DNA-loop-extruding enzymes.  October 2012.
  * [7] E. Feinstein, C. Danilowicz, A. Conover, R. Gunaratne, N. Kleckner, and M. Prentiss.  Single-molecule studies of the stringency factors and rates governing the polymerization of RecA on double-stranded DNA.  Nucleic Acids Research, 39(9):3781–3791, May 2011.
  * [8] Geoffrey Fudenberg and Leonid A. Mirny.  Higher-order chromatin structure: bridging physics and biology.  Current Opinion in Genetics & Development, 22(2):115–124, April 2012.
  * [9] Mahipal Ganji, Indra A. Shaltiel, Shveta Bisht, Eugene Kim, Ana Kalichava, Christian H. Haering, and Cees Dekker.  Real-time imaging of DNA loop extrusion by condensin.  Science, 360(6384):102–105, April 2018.
  * [10] Anton Goloborodko, Maxim V Imakaev, John F Marko, and Leonid Mirny.  Compaction and segregation of sister chromatids via active loop extrusion.  eLife, 5:e14864, May 2016.
  * [11] Yunyan He, Josh Lawrimore, Diana Cook, Elizabeth Erin Van Gorder, Solenn Claire De Larimat, David Adalsteinsson, M. Gregory Forest, and Kerry Bloom.  Statistical mechanics of chromosomes: in vivo and in silico approaches reveal high-level organization and structure arise exclusively through mechanical feedback between loop extruders and chromatin substrate properties.  Nucleic Acids Research, 48(20):11284–11303, November 2020.
  * [12] Caitlin Hult, David Adalsteinsson, Paula A. Vasquez, Josh Lawrimore, Maggie Bennett, Alyssa York, Diana Cook, Elaine Yeh, Mark Gregory Forest, and Kerry Bloom.  Enrichment of dynamic chromosomal crosslinks drive phase separation of the nucleolus.  Nucleic Acids Research, 45(19):11159–11173, November 2017.
  * [13] Alexander A. Kukreja, Sisira Kavuri, and Ajit P. Joglekar.  Microtubule Attachment and Centromeric Tension Shape the Protein Architecture of the Human Kinetochore.  Current Biology, 30(24):4869–4881.e5, December 2020.
  * [14] Yang Li, Shengyuan Xu, Jinqiao Duan, Xianbin Liu, and Yuming Chu.  A machine learning method for computing quasi-potential of stochastic dynamical systems.  Nonlinear Dynamics, 109(3):1877–1886, August 2022.
  * [15] Bo Lin, Qianxiao Li, and Weiqing Ren.  A Data Driven Method for Computing Quasipotentials.  In Proceedings of Machine Learning Research, volume 145, pages 652–670, 2022.
  * [16] Andrew D. McAinsh and Geert J. P. L. Kops.  Principles and dynamics of spindle assembly checkpoint signalling.  Nature Reviews Molecular Cell Biology, 24(8):543–559, August 2023\.
  * [17] Jay Newby, Jennifer L. Schiller, Timothy Wessler, Jasmine Edelstein, M. Gregory Forest, and Samuel K. Lai.  A blueprint for robust crosslinking of mobile species in biogels with weakly adhesive molecular anchors.  Nature Communications, 8(1):833, October 2017.
  * [18] Matthew Robert Paul, Andreas Hochwagen, and Sevinç Ercan.  Condensin action and compaction.  Current genetics, 65(2):407–415, April 2019.
  * [19] Grigorios A. Pavliotis and Andrew M. Stuart.  Multiscale Methods: Averaging and Homogenization.  Springer, New York, NY, 2008.
  * [20] Stephanie Andrea Schalbetter, Anton Goloborodko, Geoffrey Fudenberg, Jon-Matthew Belton, Catrina Miles, Miao Yu, Job Dekker, Leonid Mirny, and Jonathan Baxter.  SMC complexes differentially compact mitotic chromosomes according to genomic context.  Nature Cell Biology, 19(9):1071–1080, September 2017.
  * [21] Tsuyoshi Terakawa, Shveta Bisht, Jorine M. Eeftens, Cees Dekker, Christian H. Haering, and Eric C. Greene.  The condensin complex is a mechanochemical motor that translocates along DNA.  Science (New York, N.Y.), 358(6363):672–676, November 2017.
  * [22] Nynke L. van Berkum, Erez Lieberman-Aiden, Louise Williams, Maxim Imakaev, Andreas Gnirke, Leonid A. Mirny, Job Dekker, and Eric S. Lander.  Hi-C: A Method to Study the Three-dimensional Architecture of Genomes.  Journal of Visualized Experiments : JoVE, (39):1869, May 2010.
  * [23] Paula A Vasquez and Kerry Bloom.  Polymer models of interphase chromosomes.  Nucleus, 5(5):376–390, September 2014.
  * [24] Paula A. Vasquez, Caitlin Hult, David Adalsteinsson, Josh Lawrimore, Mark G. Forest, and Kerry Bloom.  Entropy gives rise to topologically associating domains.  Nucleic Acids Research, 44(12):5540–5549, July 2016.
  * [25] Paula A. Vasquez, Ben Walker, Kerry Bloom, Daniel Kolbin, Neall Caughman, Ronit Freeman, Martin Lysy, Caitlin Hult, Katherine A. Newhall, Micah Papanikolas, Christopher Edelmaier, and M. Gregory Forest.  The power of weak, transient interactions across biology: A paradigm of emergent behavior.  Physica D: Nonlinear Phenomena, 454:133866, November 2023.
  * [26] Benjamin Walker, Dane Taylor, Josh Lawrimore, Caitlin Hult, David Adalsteinsson, Kerry Bloom, and M. Gregory Forest.  Transient crosslinking kinetics optimize gene cluster interactions.  PLOS Computational Biology, 15(8):e1007124, August 2019.
  * [27] Benjamin L. Walker and Katherine A. Newhall.  Numerical computation of effective thermal equilibria in stochastically switching Langevin systems.  Physical Review E, 105(6):064113, June 2022.
  * [28] Timothy Wessler, Alex Chen, Scott A. McKinley, Richard Cone, M. Gregory Forest, and Samuel K. Lai.  Using Computational Modeling To Optimize the Design of Antibodies That Trap Viruses in Mucus.  ACS Infectious Diseases, 2(1):82–92, January 2016.
  * [29] Joseph Xu Zhou, M. D. S. Aliyu, Erik Aurell, and Sui Huang.  Quasi-potential landscape in complex multi-stable systems.  Journal of the Royal Society Interface, 9(77):3539–3553, December 2012.