# Is Anisotropy Inherent to Transformers?

Nathan Godey$\mkern 6.0mu{}^{1,2}$ Éric de la Clergerie1 Benoît Sagot1
1Inria, Paris, France
2Sorbonne Université, Paris, France
<EMAIL_ADDRESS>
###### Abstract

The representation degeneration problem is a phenomenon that is widely
observed among self-supervised learning methods based on Transformers. In NLP,
it takes the form of anisotropy, a singular property of hidden representations
which makes them unexpectedly close to each other in terms of angular distance
(cosine-similarity). Some recent works tend to show that anisotropy is a
consequence of optimizing the cross-entropy loss on long-tailed distributions
of tokens. We show in this paper that anisotropy can also be observed
empirically in language models with specific objectives that should not suffer
directly from the same consequences. We also show that the anisotropy problem
extends to Transformers trained on other modalities. Our observations tend to
demonstrate that anisotropy might actually be inherent to Transformers-based
models.

##  1 Introduction

In recent years, deep learning models based on Transformers have led to
significant breakthroughs in the field of natural language processing (NLP).
These models have demonstrated state-of-the-art performance across a range of
tasks, such as language modeling, machine translation, and sentiment analysis.
However, despite their successes, these models suffer from a phenomenon known
as the representation degeneration problem. Specifically, this problem
manifests as anisotropy, a property of hidden representations that makes them
all close to each other in terms of angular distance (cosine-similarity).

Anisotropy has been widely observed among self-supervised learning methods
based on Transformers, and recent research suggests that it may be a
consequence of optimizing the cross-entropy loss on long-tailed distributions
of tokens (Biś et al., 2021). However, it remains uncertain whether anisotropy
is a fundamental property of Transformers-based models or a consequence of the
pre-training process.

In this paper, we investigate the anisotropy problem in depth, and we make
several contributions:

  * •

We demonstrate empirically that anisotropy can be observed in language models
with character-aware architectures that should not suffer directly from the
same consequences as token-based NLP;

  * •

We extend our observations to Transformers trained on other modalities, such
as image and audio data;

  * •

We provide empirical observations on the anisotropy problem in the Transformer
block and argue that it may be an inherent property of the self-attention
mechanism.

##  2 Related Work

Figure 1: Average cosine-similarity between hidden representations across
layers for token-level NLP models. For T5-base, we concatenate encoder and
decoder results.

The general phenomenon of anisotropy in token-based Transformers for language
models has been shown in Ethayarajh (2019). Figure 1 extends one of their
experiment to more architectures. Gao et al. (2019) shows that the
degeneration of representations come from the distributions of subwords in
natural language, namely the existence of unused and rare tokens, that tend to
push all representations away from the origin towards a specific direction.

Other works have established a connection between word frequency and
distortions of the latent spaces (Yu et al., 2022; Puccetti et al., 2022;
Rajaee and Pilehvar, 2022). Biś et al. (2021) have shown that anisotropy in
LMs could be explained by a global drift of the representations in the same
direction, thus unifying conclusions from Ethayarajh (2019) and Gao et al.
(2019). The authors propose that this drift is caused by the persistent
updating of the representation of rare and unused tokens in a consistent
direction, due to the nature of the softmax operation in the cross-entropy
loss. They show that removing the average component to all representations
leads to a nearly perfect isotropy.

Several methods have been proposed to reduce anisotropy in Transformers-based
LMs at token-level (Rajaee and Pilehvar, 2021; Wang et al., 2020), or at
sentence-level (Gao et al., 2021; Yan et al., 2021; Su et al., 2021). They
usually consist in post-processing the representations, and lead to downstream
performance boosts. We argue that these positive results are paving the way
for the search of pre-training objectives that do not introduce anisotropy in
the first place, in the hope that the resulting models will also be more
performant without any post-processing, and potentially be trained more
efficiently.

##  3 Anisotropy in pre-trained Transformers

###  3.1 Character-based NLP

Figure 2: Average cosine-similarity between hidden representations across
layers for character-level models.

To validate that the cross-entropy objective applied on vocabularies
containing rare tokens is the main cause for the common drift issue, we
explore anisotropy in character-based models. We study different
architectures:

  * •

CharacterBERT (El Boukkouri et al., 2020) is constructing whole word
representations from character embeddings put through convolutions and highway
layers, before feeding them to a Transformers architecture.

  * •

CANINE (Clark et al., 2022) is downsampling contextualized character
representations via a strided convolution before feeding them to a
Transformers. It can be trained either with a subword-based objective
(CANINE-s) or with a character-level one (CANINE-c).

  * •

MANTa-LM (Godey et al., 2022) is based on a differentiable segmentation and
embedding module added before an encoder-decoder model in the style of T5
(Raffel et al., 2020). It takes bytes as inputs and outputs, but builds
internal representations that are usually based on several bytes.

  * •

ByT5 (Xue et al., 2022) is a version of T5 that is trained at byte-level. To
afford for more complex encoding, the authors resize the encoder-decoder
architecture.

Neither of these architectures should suffer from out-of-vocabulary in the
process of creating representations. The models that predict at word or sub-
word level (CharacterBERT and CANINE-s) could have the cross-entropy loss
systematically pushing away rare item representations. However, it is rather
unclear why it would imply an embedding drift at deeper layers. Hence, if
anisotropy was only caused by the presence of unused or rare subwords, those
character-level models should be much less prone to this issue.

To verify this hypothesis, we compute hidden representations for the
validation set of the WikiText-103 corpus (Merity et al., 2016). We then
compute the average cosine-similarity between two representations, uniformly
taken in the whole validation corpus.

In fact, as we can see in Figure 2, those models all display significant
levels of anisotropy at least in one of their layers. Interestingly, the
models that are based solely on characters or bytes for input and prediction
(ByT5 and CANINE-c) seem to display even higher levels of anisotropy. We note,
as it is the case for the T5 model, that the ByT5 decoder displays extremely
high levels of anisotropy.

###  3.2 Other modalities

(a) Speech

(b) Vision

Figure 3: Average cosine-similarity between hidden representations across
layers for Speech and Vision modalities. We observe that across both
modalities, several models display significant levels of anisotropy.

We’ve shown in the previous section that character-level language models
suffer from anisotropy similarly to token-level ones. We proceed to explore
the anisotropy problem for Transformers-based models in other modalities,
specifically speech and vision.

For speech models, we consider wav2Vec 2.0 (Baevski et al., 2020), HuBERT (Hsu
et al., 2021), and Whisper (Radford et al., 2022) with the Common Voice 11.0
dataset (Ardila et al., 2020). For vision models, we use ViT (Wu et al.,
2020), BEiT (Bao et al., 2021), MiT (Xie et al., 2021), and DEiT (Touvron et
al., 2021) on the ImageNet dataset (Russakovsky et al., 2015).

As in subsection 3.1, we infer hidden representations on the validation sets
for each modality. We then uniformly sample pairs of vectors to get cosine-
similarity values for every layer of every model. The results are displayed in
Figure 3.

Once again, almost every model shows a significant level of anisotropy on some
of its layers. Notably, speech models show very anisotropic representations,
as every layer of every model outputs an average cosine-similarity of at least
$0.2$. We find some exceptions among vision models, since the MiT model seems
to use only isotropic representations and the ViT model has a low average
cosine-similarity for all its layers.

We also conduct the same experiment for convolution-based networks in the
vision modality. The models at glance are ResNet (He et al., 2016),
EfficientNet (Tan and Le, 2019), CvT (Wu et al., 2021), ConvNeXt (Liu et al.,
2022), and VAN (Guo et al., 2022). For these networks, we flatten convolution
maps before computing the cosine-similarity.

Figure 4: Average cosine-similarity between hidden representations across
layers for convolution-based vision models.

We observe in Figure 4 that most of the convolution-based models are
isotropic. Interestingly, the only exception is ResNet-50, whose
representations become more and more isotropic as one explores deeper layers.
This could partially be explained by the fact that the batch normalization
(Ioffe and Szegedy, 2015) used in some of these models mitigates a posteriori
the drift effect by removing the mean component of the representations.
However, the ConvNeXt model also seems to use isotropic representations while
not using batch normalization, which shows that this is not the only factor in
the isotropic behavior of these models.

###  3.3 To drift or not to drift?

Related works (Biś et al., 2021; Gao et al., 2019) show that anisotropy in
subword-level language models is caused by a drift of the hidden
representations in a shared direction. In this section, we try to extend this
observation to other modalities.

We study the correlation between the uniformly measured cosine-similarity, and
the norm of the average hidden representation $||\bar{x}||_{2}$ for each
layer. If anisotropy could be explained by the drift effect, we would expect a
monotonic relation between $||\bar{x}||_{2}$ and the average cosine-
similarity. To verify this, we apply a Spearman correlation test on these two
metrics for every model from subsection 3.1 and subsection 3.2, along with
some token-level language models, namely T5 (Raffel et al., 2020), BERT
(Devlin et al., 2019), RoBERTa (Liu et al., 2019), and GPT-2 (Radford et al.,
2019).

Figure 5: p-value of the Spearman correlation test between the norm of the
average representation and the cosine-similarity averaged over all layers,
across modalities. Models above the red dotted line are not significantly
affected by the drift effect.

In Figure 5, we observe that we cannot correlate the anisotropy level and the
drift component for several models. We notice that the anisotropy affecting
some CNN-based vision models can be explained by the drift effect, but we
underline that this is not the case for Tranformers-based models in the same
modality. Some speech models (HuBERT and Whisper-base) also display signs of
anisotropy that cannot be correlated with the drift effect. On the other hand,
the anisotropy of subword-based models can generally be explained by the drift
effect, except for GPT-2 for which the Spearman correlation metric may not be
appropriate. We provide a similar analysis based on the Pearson correlation
test and discuss the relevance of each statistic in Appendix A.

Interestingly, Figure 5 displays a correlation for all character-based models
but Canine-C and MANTa-base. All in all, if Figure 5 confirms the conclusions
of related works when it comes to token-based language models, it contradicts
the hypothesis that the anisotropy in Transformers-based models can be
generally be explained by the drift effect.

##  4 Explaining the representation drift

In this section, we expose some properties of the Transformer block as used in
BERT (Devlin et al., 2019). We analyze experimentally the behavior of the
Transformer block $T$ when a common bias term $b$ is added to the input
representations $\mathbf{x}$. This allows us to mimic the common drift as
mentioned in Biś et al. (2021) and to identify some properties induced by this
artificial drift on the output representations.

###  4.1 Experimental setup

We consider an embedding lookup table $E$ and a Transformer block $T$ built
and initialized as in BERT (Devlin et al., 2019). We then draw 16 input
embedding sequences $\mathbf{x}$ of length 512 uniformly from $E$. To account
for a drift component of norm $N\in\mathbb{R}$, we generate a vector
$b_{u}\sim\mathcal{N}(0,I_{d})$, which we normalize into
$b=\frac{b_{u}}{||b_{u}||_{2}}\times N$. We finally compute
$T(\mathbf{x}_{i}+b)$ for every sequence $x_{i}$, and study the resulting
distributions.

Specifically, we study the average norm of the input representations
$E(||\mathbf{x}_{i}+b||_{2})$ against the average norm of the output
representations $E(||T(\mathbf{x}_{i}+b)||_{2})$ in 6(b). We also retrieve the
self-attention scores before the softmax operation, namely
$\frac{QK^{T}}{\sqrt{d_{k}}}$, along with the corresponding $Q$ and $K$
matrices. We study some of their properties in Figure 7 and Figure 8.

###  4.2 Input vs. output analysis

(a) Cosine similarity

(b) Norm

Figure 6: Input/Output comparison of a Transformer block from BERT-base as the
bias norms increases.

In 6(a), we observe that the output representations have an average cosine-
similarity value that is slightly higher than the one of the input
representations, no matter the level of input bias. We also notice that while
the norm of the average output representation increases with the bias norm, it
seems to meet the corresponding input measure for a given bias norm.

Interestingly, this shows that there is a fixed point in terms of norm in the
Transformers function with biased input. More formally, there seems to exist a
bias norm $N^{*}\in\mathbb{R}_{+}$ such that

$E_{x,b_{N^{*}}}(||x_{i}+b_{N^{*}}||)=E_{x,b_{N^{*}}}(||T(x_{i}+b_{N^{*}})||)$

Moreover, this fixed point level $N^{*}$ is in the order of magnitude of the
average hidden state norms of the layers of the trained BERT model. This hints
that the model’s representations stabilize when their norm is close to this
fixed point.

###  4.3 Exploring the Transformer block

To understand the effect of the drift effect on the inner workings of the
Transformer layer, we take a closer look at the self-attention operation as
the average input representation drifts away.

Figure 7: Histograms of the pre-softmax attention scores as the input bias
norm increases. Other initializations of the layer and of the bias direction
$b_{u}$ led to a general increase of the attention scores instead.

Figure 7 shows that the attention scores tend to move away from zero as the
input bias norm increases. Indeed, as the norm of the average $\bar{x}$ of the
input embeddings increases, we can expect the query and key vectors $Q$ and
$K$ to also display signs of anisotropy. Actually, for each self-attention
head, and for all position $i\in[1,L]$, one should observe:

$\begin{cases}E_{x}(Q_{i})=W_{Q}\bar{x}+b_{Q}\\\
E_{x}(K_{i})=W_{K}\bar{x}+b_{K}\end{cases}$ (1)

We can observe in Figure 8 that query and key representations indeed increase
in norm with the input bias norm. We also notice that the corresponding
distributions are anisotropic even when no bias is added, which may be a
consequence of BERT’s initialization parameters.

(a) Cosine sim.

(b) Norm

Figure 8: Analysis of the self-attention query and key distributions

###  4.4 Impact of the drift

After exploring the consequences of the drift of input representations on the
query-key product in self-attention, we identify in this section the
implications of this drift at a more explainable level.

Figure 9: Evolution of the self-attention softmax values as the input bias
norm increases.

In Figure 9, we retrieve softmax values in the self-attention block and for
each position, we extract the maximum, the median and the minimum. We then
average these values over the whole batch, and repeat for various input bias
norm levels. We notice that as the input bias norm increases, the self-
attention softmax distributions become more and more categorical.

(a) Maximum

(b) Minimum

Figure 10: Comparison of the extreme values of each sequence averaged over the
batch as the bias norm increases.

This becomes even clearer if we take the maximum and minimum values over the
whole sequences, as in Figure 10.

However, at low anisotropy levels, i.e. when the bias norm is low, we see that
the effect is not very important. Figure 9 and Figure 10 only hint at the fact
that the drift of embeddings may help the self-attention to be more
categorical. Another explanation could be that training favors categorical
self-attention patterns, as has been pointed out in previous works (Clark et
al., 2019), which in turn induces a drift in the models’ representations. We
leave these questions for future work.

##  5 Discussion

In this work, we hypothesize that the drift effect in Transformers-based
models is not solely responsible for the anisotropy observed in most hidden
representations across modalities. As section 4 shows, untrained Transformers
layers display a tendency towards anisotropy. Biased inputs tend to increase
the variance of the attention scores and thus facilitate the emergence of
categorical patterns in the self-attention mechanisms.

To strengthen our hypothesis, an analysis of anisotropy via drift during the
pre-training process would help understand how the query and key distributions
evolve and verify that the progressive drift of the input representations
leads to a progressive increase of the sharpness of the self-attention
behavior. We found that such a study would require a more complex experimental
setup, as the input data distribution would then need to be considered
carefully, and as the Gaussian sample of a bias term becomes less relevant
since a specific drift direction naturally emerges in the representations.

We hope that this work paves the way for a clearer understanding of the
inherent anisotropy in Transformer-based models. Even though anisotropy has
not been shown to be an issue in language modeling, previous works have
advocated that removing anisotropy in output representations leads to better
sense disambiguation abilities (Bihani and Rayz, 2021; Biś et al., 2021). It
appears to us that re-thinking the self-attention operation to facilitate the
emergence of categorical patterns with minimal implications for the input
representations could be an interesting continuation of our work.

## Conclusion

In this paper, we investigated the anisotropy problem through the lens of the
drift effect, and made several contributions to the understanding of this
phenomenon. We demonstrated that anisotropy can be observed in language models
with character-aware architectures, extended our observations to Transformers
trained on other modalities, and studied anisotropy in untrained Transformers
layers.

We conclude that anisotropy almost systematically affects Transformers on all
modalities, in a way that is not always correlated with the drift of the
representations. We also provide empirical evidence that anisotropy may be an
inherent property of the self-attention mechanism when modeling categorical
attention patterns. We hypothesize that a revision of the self-attention
operation could help reduce anisotropy by facilitating the emergence of sharp
self-attention softmax distributions without modifying the geometry of the
hidden representations.

## Limitations

As mentioned in the Discussion section, we acknowledge that section 4 does not
take into account the training dynamics, and only exposes some properties of
the Transformer layer at initialization. We also notice that the Spearman
correlation test used in Figure 5 may not be well-suited for such noisy
observations, as the high p-value of the GPT-2 model shows. We provide a
similar graph based on the Pearson correlation in Appendix A.

## Ethics Statement

To the best of our knowledge, our work does not raise any ethical concern.
However, as noted in Zhou et al. (2021), we believe that distortions in the
embedding space may be related to bias in the training data, whether it is
inherent to the structure of the modality (e.g. the Zipfian distribution of
words), or due to human factors (e.g. geographical considerations).

## Acknowledgements

This work was funded by the last authors’ chair in the PRAIRIE institute
funded by the French national agency ANR as part of the “Investissements
d’avenir” programme under the reference ANR-19-P3IA-0001.

We would like to thank Roman Castagné for useful discussions that led to
focusing on observing the effect of anisotropy in the self-attention process.

## References

  * Ardila et al. (2020) R. Ardila, M. Branson, K. Davis, M. Henretty, M. Kohler, J. Meyer, R. Morais, L. Saunders, F. M. Tyers, and G. Weber. 2020.  Common voice: A massively-multilingual speech corpus.  In _Proceedings of the 12th Conference on Language Resources and Evaluation (LREC 2020)_ , pages 4211–4215.
  * Baevski et al. (2020) Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed, and Michael Auli. 2020.  wav2vec 2.0: A framework for self-supervised learning of speech representations.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 12449–12460. Curran Associates, Inc.
  * Bao et al. (2021) Hangbo Bao, Li Dong, and Furu Wei. 2021.  Beit: BERT pre-training of image transformers.  _CoRR_ , abs/2106.08254.
  * Bihani and Rayz (2021) Geetanjali Bihani and Julia Rayz. 2021.  Low anisotropy sense retrofitting (LASeR) : Towards isotropic and sense enriched representations.  In _Proceedings of Deep Learning Inside Out (DeeLIO): The 2nd Workshop on Knowledge Extraction and Integration for Deep Learning Architectures_ , pages 81–95, Online. Association for Computational Linguistics.
  * Biś et al. (2021) Daniel Biś, Maksim Podkorytov, and Xiuwen Liu. 2021.  Too much in common: Shifting of embeddings in transformer language models and its implications.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 5117–5130, Online. Association for Computational Linguistics.
  * Clark et al. (2022) Jonathan H. Clark, Dan Garrette, Iulia Turc, and John Wieting. 2022.  Canine: Pre-training an efficient tokenization-free encoder for language representation.  _Transactions of the Association for Computational Linguistics_ , 10:73–91.
  * Clark et al. (2019) Kevin Clark, Urvashi Khandelwal, Omer Levy, and Christopher D. Manning. 2019.  What does BERT look at? an analysis of BERT’s attention.  In _Proceedings of the 2019 ACL Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP_ , pages 276–286, Florence, Italy. Association for Computational Linguistics.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * El Boukkouri et al. (2020) Hicham El Boukkouri, Olivier Ferret, Thomas Lavergne, Hiroshi Noji, Pierre Zweigenbaum, and Jun’ichi Tsujii. 2020.  CharacterBERT: Reconciling ELMo and BERT for word-level open-vocabulary representations from characters.  In _Proceedings of the 28th International Conference on Computational Linguistics_ , pages 6903–6915, Barcelona, Spain (Online). International Committee on Computational Linguistics.
  * Ethayarajh (2019) Kawin Ethayarajh. 2019.  How contextual are contextualized word representations? Comparing the geometry of BERT, ELMo, and GPT-2 embeddings.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 55–65, Hong Kong, China. Association for Computational Linguistics.
  * Gao et al. (2019) Jun Gao, Di He, Xu Tan, Tao Qin, Liwei Wang, and Tie-Yan Liu. 2019.  Representation degeneration problem in training natural language generation models.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net.
  * Gao et al. (2021) Tianyu Gao, Xingcheng Yao, and Danqi Chen. 2021.  SimCSE: Simple contrastive learning of sentence embeddings.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 6894–6910, Online and Punta Cana, Dominican Republic. Association for Computational Linguistics.
  * Godey et al. (2022) Nathan Godey, Roman Castagné, Éric de la Clergerie, and Benoît Sagot. 2022.  MANTa: Efficient gradient-based tokenization for end-to-end robust language modeling.  In _Findings of the Association for Computational Linguistics: EMNLP 2022_ , pages 2859–2870, Abu Dhabi, United Arab Emirates. Association for Computational Linguistics.
  * Guo et al. (2022) Meng-Hao Guo, Cheng-Ze Lu, Zheng-Ning Liu, Ming-Ming Cheng, and Shi-Min Hu. 2022\.  Visual attention network.  _arXiv preprint arXiv:2202.09741_.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. 2016.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778.
  * Hsu et al. (2021) Wei-Ning Hsu, Benjamin Bolte, Yao-Hung Hubert Tsai, Kushal Lakhotia, Ruslan Salakhutdinov, and Abdelrahman Mohamed. 2021.  Hubert: Self-supervised speech representation learning by masked prediction of hidden units.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 29:3451–3460.
  * Ioffe and Szegedy (2015) Sergey Ioffe and Christian Szegedy. 2015.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pages 448–456, Lille, France. PMLR.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized BERT pretraining approach.  _CoRR_ , abs/1907.11692.
  * Liu et al. (2022) Zhuang Liu, Hanzi Mao, Chao-Yuan Wu, Christoph Feichtenhofer, Trevor Darrell, and Saining Xie. 2022.  A convnet for the 2020s.  _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_.
  * Merity et al. (2016) Stephen Merity, Caiming Xiong, James Bradbury, and Richard Socher. 2016.  Pointer sentinel mixture models.  _CoRR_ , abs/1609.07843.
  * Puccetti et al. (2022) Giovanni Puccetti, Anna Rogers, Aleksandr Drozd, and Felice Dell’Orletta. 2022\.  Outlier dimensions that disrupt transformers are driven by frequency.  In _Findings of the Association for Computational Linguistics: EMNLP 2022_ , pages 1286–1304, Abu Dhabi, United Arab Emirates. Association for Computational Linguistics.
  * Radford et al. (2022) Alec Radford, Jong Wook Kim, Tao Xu, Greg Brockman, Christine McLeavey, and Ilya Sutskever. 2022.  Robust speech recognition via large-scale weak supervision.
  * Radford et al. (2019) Alec Radford, Jeff Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019.  Language models are unsupervised multitask learners.
  * Raffel et al. (2020) Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu. 2020.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _Journal of Machine Learning Research_ , 21(140):1–67.
  * Rajaee and Pilehvar (2021) Sara Rajaee and Mohammad Taher Pilehvar. 2021.  A cluster-based approach for improving isotropy in contextual embedding space.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 2: Short Papers)_ , pages 575–584, Online. Association for Computational Linguistics.
  * Rajaee and Pilehvar (2022) Sara Rajaee and Mohammad Taher Pilehvar. 2022.  An isotropy analysis in the multilingual BERT embedding space.  In _Findings of the Association for Computational Linguistics: ACL 2022_ , pages 1309–1316, Dublin, Ireland. Association for Computational Linguistics.
  * Russakovsky et al. (2015) Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei. 2015.  ImageNet Large Scale Visual Recognition Challenge.  _International Journal of Computer Vision (IJCV)_ , 115(3):211–252.
  * Su et al. (2021) Jianlin Su, Jiarun Cao, Weijie Liu, and Yangyiwen Ou. 2021.  Whitening sentence representations for better semantics and faster retrieval.  _arXiv preprint arXiv:2103.15316_.
  * Tan and Le (2019) Mingxing Tan and Quoc V. Le. 2019.  Efficientnet: Rethinking model scaling for convolutional neural networks.  _ArXiv_ , abs/1905.11946.
  * Touvron et al. (2021) Hugo Touvron, Matthieu Cord, Matthijs Douze, Francisco Massa, Alexandre Sablayrolles, and Herve Jegou. 2021.  Training data-efficient image transformers and distillation through attention.  In _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pages 10347–10357. PMLR.
  * Wang et al. (2020) Lingxiao Wang, Jing Huang, Kevin Huang, Ziniu Hu, Guangtao Wang, and Quanquan Gu. 2020.  Improving neural language generation with spectrum control.  In _International Conference on Learning Representations_.
  * Wu et al. (2020) Bichen Wu, Chenfeng Xu, Xiaoliang Dai, Alvin Wan, Peizhao Zhang, Masayoshi Tomizuka, Kurt Keutzer, and Péter Vajda. 2020.  Visual transformers: Token-based image representation and processing for computer vision.  _ArXiv_ , abs/2006.03677.
  * Wu et al. (2021) Haiping Wu, Bin Xiao, Noel Codella, Mengchen Liu, Xiyang Dai, Lu Yuan, and Lei Zhang. 2021.  Cvt: Introducing convolutions to vision transformers.  _arXiv preprint arXiv:2103.15808_.
  * Xie et al. (2021) Enze Xie, Wenhai Wang, Zhiding Yu, Anima Anandkumar, Jose M. Alvarez, and Ping Luo. 2021.  Segformer: Simple and efficient design for semantic segmentation with transformers.  _CoRR_ , abs/2105.15203.
  * Xue et al. (2022) Linting Xue, Aditya Barua, Noah Constant, Rami Al-Rfou, Sharan Narang, Mihir Kale, Adam Roberts, and Colin Raffel. 2022.  ByT5: Towards a token-free future with pre-trained byte-to-byte models.  _Transactions of the Association for Computational Linguistics_ , 10:291–306.
  * Yan et al. (2021) Yuanmeng Yan, Rumei Li, Sirui Wang, Fuzheng Zhang, Wei Wu, and Weiran Xu. 2021.  ConSERT: A contrastive framework for self-supervised sentence representation transfer.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 5065–5075, Online. Association for Computational Linguistics.
  * Yu et al. (2022) Sangwon Yu, Jongyoon Song, Heeseung Kim, Seongmin Lee, Woo-Jong Ryu, and Sungroh Yoon. 2022.  Rare tokens degenerate all tokens: Improving neural text generation via adaptive gradient gating for rare token embeddings.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 29–45, Dublin, Ireland. Association for Computational Linguistics.
  * Zhou et al. (2021) Kaitlyn Zhou, Kawin Ethayarajh, and Dan Jurafsky. 2021.  Frequency-based distortions in contextualized word embeddings.  _CoRR_ , abs/2104.08465.

##  Appendix A Pearson correlation of the drift norm and anisotropy

Figure 11: p-value of the Pearson correlation test between the norm of the
average representation and the cosine-similarity averaged over all layers,
across modalities. Models above the red dotted line are not significantly
affected by the drift effect.

The Pearson test measures a linear correlation between random variables, while
the Spearman test measures a monotonic correlation. As there is no specific
argument in favor of a linear relationship between the measured distributions
(average cosine-similarity and norm of the average representation), we decided
to use the Spearman correlation test in order to take into account more
complex relation patterns.

Nevertheless, this metric is based on the rank of each observation, and is
thus not robust to fluctuations due to sample variance, specifically when
working with such small samples. This is reflected by the discrepancy between
Pearson and Spearman p-values for some models (e.g. GPT-2).


# Identifying stable communities in Hi-C data using a multifractal null model

Lucas Hedström<EMAIL_ADDRESS>Integrated Science Lab, Department of
Physics, Umeå University, Umeå, Sweden  Antón Carcedo Martínez  Integrated
Science Lab, Department of Physics, Umeå University, Umeå, Sweden  Ludvig
Lizana  Integrated Science Lab, Department of Physics, Umeå University, Umeå,
Sweden

###### Abstract

Chromosome capture techniques like Hi-C have expanded our understanding of
mammalian genome 3D architecture and how it influences gene activity. To
analyze Hi-C data sets, researchers increasingly treat them as DNA-contact
networks and use standard community detection techniques to identify mesoscale
3D communities. However, there are considerable challenges in finding
significant communities because the Hi-C networks have cross-scale
interactions and are almost fully connected. This paper presents a pipeline to
distill 3D communities that remain intact under experimental noise. To this
end, we bootstrap an ensemble of Hi-C datasets representing noisy data and
extract 3D communities that we compare with the unperturbed dataset. Notably,
we extract the communities by maximizing local modularity (using the
Generalized Louvain method), which considers the multifractal spectrum
recently discovered in Hi-C maps. Our pipeline finds that stable communities
(under noise) typically have above-average internal contact frequencies and
tend to be enriched in active chromatin marks. We also find they fold into
more nested cross-scale hierarchies than less stable ones. Apart from
presenting how to systematically extract robust communities in Hi-C data, our
paper offers new ways to generate null models that take advantage of the
network’s multifractal properties. We anticipate this has a broad
applicability to several network applications.

##  I Introduction

Mammalian genomes fold into complex 3D structures that help regulate genetic
processes such as transcription, DNA replication and repair, and epigenetics
[1, 2, 3, 4, 5]. Several of these insights derive from chromosome capture
techniques, such as Hi-C [6, 7, 8], that quantifies the number of physical
contacts between all DNA segment pairs in the genome across a cell population.
Data from these experiments allowed constructing comprehensive chromosome-wide
3D interaction maps, highlighting local 3D structures with high internal
frequencies, such as Topologically Associated Domains (TADs) and the large-
scale binary division into A/B compartments generally associated with active
(A) and inactive (B) chromatin [9, 10].

However, these 3D structures represent but two instances on a wide spectrum of
structures that fold into each other, forming complex, often nested,
hierarchies [11, 12, 13, 14, 15]. These hierarchies make it difficult to
design robust TAD-finding algorithms (”TAD callers”) and for the field to
agree on common TAD definitions [16, 17]. This ambiguity is likely one reason
there exists a large collection of TAD callers [18], each associated with
their own parameters and arbitrary choices. In practice, this means that two
callers may disagree on the TAD boundaries using the same dataset, let alone
two datasets deriving from theoretically identical replicate experiments,
including experimental and biological noise [19].

To help resolve some of these issues, we recently proposed a method to extract
robust 3D communities in Hi-C data [20]. In that paper, we used a network
community detection algorithm (Generalized Louvain) to generate an ensemble of
feasible network partitions from a public Hi-C data set [21] and quantified
the most conserved node-community relationships. This approach distinguishes
robust from variable communities and tracks how the robustness changes with
the network scale (and chromatin state). We found that robustness is highly
scale-dependent and typically worse for small communities, including TADs.
When robustness is low, there are many ways to split the network into
mathematically feasible community partitions, which is likely one reason TAD
callers disagree, as each focuses on slightly different data features [22].

This paper extends our study [20] by addressing how community robustness
changes under experimental noise rather than variability in the community
detection method. To this end, we first estimate noise levels from actual Hi-C
maps (from humans [21]) and construct an ensemble of blurred maps by adding
noise to each contact (or pixel). We call these bootstrapped maps [23]. Next,
we extract 3D communities by considering the Hi-C data as a network and using
the Generalized Louvain algorithm. Finally, we quantify the overlap between
the bootstrapped maps relative to the unperturbed map. We find that stable
communities typically have above-average internal contact frequencies, tend to
be enriched in active chromatin marks, and fold into more nested cross-scale
hierarchies than less stable ones.

Figure 1: A schematic showcasing the methodology that computes the community
stability for a Hi-C network. Using the original Hi-C data, we bootstrap by
sampling data from the contact distributions associated with each distance
(Fig. 2). We then run a community detection algorithm on the original Hi-C
data to get the original partition (GenLouvain). To reduce the number of
samples required for good statistics, we use the original partition as an
initial partition for GenLouvain. Each sample gives a list of partitions, each
with unique communities, shown as boxes along the diagonal colored by the
community number ID. Based on the overlap (calculated with the Jaccard index),
we count the number of times each community in the original Hi-C data appears
in the bootstrapped data, which we use as a proxy for stability. (rightmost
panel) Original partition, where a darker color indicates communities with
higher stability.

In addition, we made yet another improvement that reaches beyond Hi-C data.
This improvement takes advantage of the network’s multifractal spectrum [24],
in a way we anticipate useful for other community detection approaches. When
using Generalized Louvain, or other methods resting on maximum modularity
[25], users must specify a background network connectivity or null model.
Simple null models like the random Newman-Girvan model [26] are common choices
in network science due to limited alternatives. However, this model does not
apply to Hi-C data since it disregards the spatial constraints associated with
the DNA segments (represented as nodes) belonging to a long linear molecule.
One way to deal with this specific problem is to introduce a distant-dependent
null model, where the contacts decay as a power law (as discovered in Hi-C)
[14, 27]. But as the power-law exponent changes with distance, too, (from
$\approx$0.75 to $\approx$1.08 [28]) this introduces a new arbitrary parameter
as to when to switch exponent. While we explored both exponents in [20], here
we remove this arbitrary choice by implementing a one-parameter null model
that we fit to the Hi-C map’s bifractal spectrum using techniques from [24].
This approach is not limited to Hi-C data and should apply to any complex
network.

##  II Methods

Our computational method has several steps (Fig. 1). First, we generate
bootstrapped samples representing noisy Hi-C data. Second, we find communities
using a method that considers the bifractal structure in Hi-C data. Third, we
step through the pipeline to determine the stable communities that are
invariant under noise.

###  II.1 Collecting and preprocessing experimental data

Chromosome contacts. We downloaded Hi-C data for the B-lymphoblastoid human
cell line (GM12878) [21] from the GEO database (MAPQG0 dataset, 100 kb
resolution) [29]. We only consider intra-chromosome contacts in our analysis.
We interpret the data file as a weighted network (”Hi-C network”) with
elements $A_{ij}$ in sparse form, where each node represents 100 kb of DNA,
and the link weight is the measured contact count.

Before constructing the Hi-C network, we normalize the data using the standard
Knight-Ruiz matrix balancing algorithm [30] so that the sum of all weights to
a node equals unity. Note that the centromere (with no measured contact
counts) must be removed before KR-normalization to avoid convergence issues.
However, since the indices on the Hi-C matrix represent physical distances, we
reintroduce the centromere before performing community detection (see Sec.
II.3).

Chromatin states. To relate our results to different functional genomic
regions, we downloaded a dataset from ENCODE [31] of 13 different chromatin
types created from a multivariate hidden Markov model (denoted HMM states)
[32, 33]. To reduce the number of states, we grouped the data into five
effective categories: Promoters, Enhancers, Transcribed regions,
Heterochromatin, and Insulators (see App. C for complete definitions in terms
of the standard HMM states). Moreover, this data set allows us to calculate
the chromatin state’s folds of enrichment (FE) for each Hi-C bin. We calculate
enrichment relative to the chromosome-wide average (App. C).

###  II.2 Bootstrapping networks

To estimate community stability across replicates, we bootstrapped data that
simulated noise from a typical Hi-C experiment. This noise comes from multiple
sources [19], such as randomness of ligations between loci [34] or random
contacts between loci due to fluctuations [6]. In Fig 2, we show histograms of
the logged Hi-C contacts $\log(A_{ij})$ for four fixed distances in chromosome
10. We note that the distributions at all these distances are log-normal,
albeit with their unique mean and standard deviation We use these empirical
distributions as a proxy for experimental noise.

To generate bootstrapped Hi-maps with noise, we calculate the standard
deviation $\sigma_{d}$ of the logged contacts for a fixed distance $d=|i-j|$
as

$\sigma_{d}=\sqrt{\langle\log^{2}(A_{ij})\rangle_{d}-\langle\log(A_{ij})\rangle_{d}^{2}},$
(1)

where $\langle\ldots\rangle_{d}$ denotes the average at the distance $d$.
Next, we bootstrap new noisy contacts $A_{ij}^{\prime}$ for each bin using

$A_{ij}^{\prime}=A_{ij}\times e^{N(0,(\epsilon\sigma_{d})^{2})}$ (2)

where $\epsilon$ is the noise amplitude, and $N(0,(\epsilon\sigma_{d})^{2})$
denotes the normal distribution with mean $0$ and variance
$(\epsilon\sigma_{d})^{2}$. As shown below, we tune the amplitude $\epsilon$
to not completely disrupt the chromosome’s community structure (e.g., the
TADs), yet yield maps with realistic noise.

###  II.3 Detecting communities using a multifractal null model

Figure 2: Histogram of the logged Hi-C contacts $A_{ij}$ at four distances for
chromosome 10 (bars). Overlaying the histogram as solid lines, we depict log-
normal distributions fitted with the mean and standard deviation from the
data.

To find Hi-C network communities, or ”3D communities” [14], we use the
Generalized Louvain method (GenLouvain) implemented in MATLAB [35]. GenLouvain
searches for network partitions that maximize the modularity function $Q$,
capturing local deviations from the expected background connectivity or null
model. In its general form, the modularity is expressed as

$Q=\frac{1}{2m}\sum_{i\neq j}\left(A_{ij}-\gamma
P_{ij}^{(0)}\right)\delta_{C_{i},C_{j}}$ (3)

where $A_{ij}$ are entries in the weighted adjacency (Hi-C) matrix, $m$ is the
total weight of this matrix, $\gamma$ is a scale parameter setting the
effective community size, $C_{i}$ is node $i$’s community assignment, and
$P_{ij}^{(0)}$ is the null model. The last factor, $\delta_{ij}$, denotes the
Kronecker-delta function and indicates that the sum calculates the partial sum
of all local modularities within a community.

To calculate $Q$ and identify communities, one must specify a null model that
captures the characteristic contact patterns in the data. The most common
choice is random connections, known as the Newman-Girvan null model [26], but
it does not account for structural aspects associated with long DNA polymer
chain folded in 3D inside the cell nucleus [36, 37]. Previous work generalized
the Newman-Girvin model to include a power-law with linear node separation $d$
as observed in Hi-C maps and predicted by theoretical polymer models [14].
However, the power-law exponent tends to change with distance [6, 28]. For
example, Within TADs, the contacts decay on average as $d^{-0.75}$, whereas it
is closer to $d^{-1.08}$ for longer distances. Therefore, if using a single
decay function as the null hypothesis, one must introduce an arbitrary cutoff
at which the decay exponent should change. To circumvent this problem, we
propose another null model embracing the multifractal properties of Hi-C maps
[24]. Admittedly, this model also has a parameter, but the fitting to Hi-C
data is more systematic.

The model reads

$P_{ij}^{(0)}=2m\frac{H_{a}(i,j)}{\sum_{i,j}H_{a}(i,j)}.$ (4)

where $H_{a}(i,j)$ represents the so-called Hierarchical domain model (HDM)
from [24], calculated as

$\displaystyle H_{a}(i,j)=$ $\displaystyle\left\\{\begin{aligned}
&\frac{1}{2^{n}-d}\sum_{k=0}^{n-1}2^{n}g(d2^{k-n})\frac{a^{k}b^{1-\delta_{kn}}}{4^{\max(n-k-1,0)}}\,\,{\rm
if}\,d<2^{n}\\\
&\sum_{k=0}^{n-1}2^{n}g(2^{k}-2^{k-n})\frac{a^{k}b^{1-\delta_{kn}}}{4^{\max(n-k-1,0)}}\,\,{\rm
otherwise}\end{aligned}\right.$ (5)

where

$g(x)=\left\\{\begin{aligned} &x\,&{\rm if}\,0\leq x\leq 0.5\\\ &1-x\,&{\rm
if}\,0.5\leq x\leq 1\\\ &0\,&{\rm otherwise.}\end{aligned}\right.$ (6)

Here, $n=\lfloor\log_{2}(N)\rfloor$, $N$ is the size of the Hi-C matrix, $a$
and $b$ are fitting paramters ($b=0.5-a$). We use a matrix formulation of
$H_{a}(i,j)$ and calculate the mean contact probability at each distance
$d=|i-j|$, serving as $P_{ij}^{(0)}$ (see App. A). Finally, we fitted the
parameter $a$ against our Hi-C data using least-squares and found $a=0.464(3)$
with minor variation between chromosomes. This value model accurately yields
the decay exponent in Hi-C contact maps at all scales (Fig. 8).

###  II.4 Reducing community detection variance

Since the GenLouvin algorithm tries to optimize the modularity function in a
Monte Carlo-like fashion, it often finds different local minima, producing
slightly different partitions with similar quality [20] (also discussed in
App. A). To reduce the number of samples required to get good statistics, we
use the partition detected in the original data as the starting partition and
then generate a collection of bootstrapped samples, see Fig. 1. This procedure
reduces the number of required samples since we do not have to account for
varying initial conditions.

##  III Results

###  III.1 Fine-tuning noise amplitudes $\epsilon$

To find realistic noise levels $\epsilon$ when bootstrapping our noisy Hi-C
maps, we benchmarked against measured TAD stabilities between replicate
experiments [16]. These experiments found the overlap of TAD boundaries
between replicates to be 62%. We generated communities using GenLouvain, where
we set the scale parameter $\gamma$ to achieve an effective community size
$\hat{s}\sim 0.88$ Mb that best agrees with the TADs considered in [16] (see
App. B). Then we generated sets of communities for $n_{\rm samp}$ bootstrapped
samples for different values of $\epsilon$. To calculate the TAD boundary
similarity between two sets of community boundaries $b_{i}$ and $b_{j}$ from
two bootstrapped samples $i$ and $j$, we use the Jaccard index

$J_{\rm TAD}(b_{i},b_{j})=\frac{|b_{i}\cap b_{j}|}{|b_{i}\cup b_{j}|},$ (7)

which varies between 1 (identical) and 0 (dissimilar).

However, Eq. (7) has a few problems when using it to quantify boundary
overlaps. For example, under extensive noise, the communities shrink until
each node is its community. This situation yields $J_{\rm
TAD}(b_{i},b_{j})\sim 1$, indicating a large overlap between community
borders, but where the partitions no longer represent TAD-sized communities.
To solve this problem, we multiply the Jaccard index in Eq. (7) by a ”penalty
factor” containing the fraction of the number of communities in the
bootstrapped sample compared to the original. This gives

$J_{\rm TAD}^{*}(b_{i},b_{j})=\frac{|b_{i}\cap b_{j}|}{|b_{i}\cup
b_{j}|}\times\left(1-\frac{||b_{i}\cup b_{j}|-|b_{o}||}{|b_{i}\cup
b_{j}|}\right)$ (8)

where $b_{o}$ are the TAD boundaries calculated from the original dataset
without noise.

Figure 3: TAD boundary stability for different noise amplitudes for all
chromosomes (color bar). The dash-dotted lines show the modified Jaccard index
in Eq. 9 (averaged over 100 samples per chromosome). As the noise $\epsilon$
increases, the overlap decreases. We note that every chromosome follows a
similar decaying pattern, where the stability drops and reaches a plateau. The
horizontal dashed line at $J_{\rm TAD}^{*}=0.62$ indicates the experimental
boundary overlap from [16]. This intersects all $J_{\rm TAD}^{*}$ lines at
slightly different critical noise amplitudes. The vertical dashed line shows
the mean of all these amplitudes, $\epsilon^{*}=0.11$.  Figure 4: Boxplots
showing the distribution of community stability at an effective community size
of $\sim 0.88$Mb across all chromosomes using $\epsilon^{*}=0.11$ and a
Jaccard index cutoff of 1.0. The kernel density estimation (KDE) for the
distribution of all stabilities over all chromosomes is shown to the right.
With these parameters, the stability median over all chromosomes is 0.40. We
note that the total distribution is bimodal, albeit with a slight maximum of
around 0.5.

To better understand how the penalty factor works, consider two illustrative
examples. One is where the bootstrapped partitions contain the same number of
communities as the original dataset but shifted along the chromosome. Here,
the overlap is low, making the term $|b_{i}\cap b_{j}|/|b_{i}\cup b_{j}|$
small. Another case is when the bootstrapped samples have many small
communities (in the extreme case, where every node represents a separate
community). Then, the overlap between bootstrapped samples will be large,
indicating a high $J_{\rm TAD}(b_{i},b_{j})$, but the number of communities
compared to the original partition differs significantly. Therefore, the last
term $||b_{i}\cup b_{j}|-|b_{o}||/|b_{i}\cup b_{j}|$ approaches $0$, leading
to small $J_{\rm TAD}^{*}(b_{i},b_{j})$, suggesting poor overlap.

Finally, to provide a global measure of TAD stability in an ensemble of
bootstrapped samples, we calculate the mean $J_{\rm TAD}^{*}(b_{i},b_{j})$
across all boundaries

$J_{\rm TAD}^{*}=\frac{2}{n_{\rm samp}(n_{\rm samp}-1)}\sum_{b_{i}\neq
b_{j}}J_{\rm TAD}^{*}(b_{i},b_{j}),$ (9)

where $n_{\rm samp}$ is the number of bootstrapped samples.

We show $J_{\rm TAD}^{*}$ over a range of noise amplitudes $\epsilon$ across
all chromosomes in Fig. 3. As expected, we find perfect overlap ($J_{\rm
TAD}^{*}\rightarrow 1$) when the noise level becomes low ($\epsilon\rightarrow
0$). Also, as $\epsilon$ increases, the overlap reduces below the experimental
benchmark (0.62, dashed horizontal line) until reaching a plateau where each
node represents separate communities. Using the 0.62-line, we may find
numerically the critical noise amplitude $\epsilon^{*}$ for each chromosome.
We note that $\epsilon^{*}$ varies between $\sim 0.05-0.15$. This parameter
presents an analog to stability since a higher $\epsilon^{*}$ also indicates
stronger resilience to noise, as the original TADs are still extractable from
the dataset. From Fig. 3, we choose $\epsilon^{*}=0.11$ as the critical noise
amplitude for all chromosomes (vertical dashed line).

###  III.2 Identifying stable communities

Once we tuned the noise amplitude, we generated several bootstrapped Hi-C maps
(see Sec. II.2) and extracted the parts of the community partition that stayed
consistent under noise relative to the original dataset. The specific
procedure reads as follows (inspired by [23]):

  * •

Get the original community partition $P_{o}$.

  * •

Generate a $n_{\rm samp}$ bootstrapped samples from the original dataset.

  * •

Get communities from each bootstrapped sample $P_{n}$, based on the starting
partition $P_{0}$.

  * •

Identify the communities from the original partition $P_{o}$ that are also in
the bootstrapped samples. To do this, we calculate the Jaccard indices between
all community pairs in each sample and pick the ones with the largest Jaccard
index.

  * •

Filter all pairs with a Jaccard index smaller than some cutoff $J_{\rm
cutoff}$ and calculate the fraction of bootstrapped samples that each
community in $P_{o}$ appears in.

Based on these steps,we calculate the stability fraction $S_{C_{i}}$, for
community $C_{i}$ in the original partition $P_{o}$ as

$S_{C_{i}}=\frac{1}{n_{\rm samp}}\sum_{n=1}^{n_{\rm
samp}}\theta\left(\max_{C_{j}\in P_{n}}J(C_{i},C_{j})-J_{\rm cutoff}\right),$
(10)

Here, $\theta(\cdot)$ denotes the Heaviside step function, and
$J(C_{i},C_{j})$ is the Jaccard index between two communities as formulated in
Eq. (7) but not specific to TADs.

The cutoff $J_{\rm cutoff}$ is a free parameter. We set it to 1.0, meaning
that we define a community $C_{i}$ as stable if all of its nodes remain in the
original dataset and the bootstrapped samples. After setting this cutoff, we
generated several bootstrapped Hi-C maps and calculated the stability fraction
$S_{C_{i}}$ (Eq. (10)) for each community $C_{i}$ associated with the original
dataset ($P_{o}$). As before, we set $\gamma$ to get the effective community
size $\hat{s}\sim 0.88Mb$ for each chromosome.

We show the distributions of the stabilities $S_{C_{i}}$ as boxplots in Fig 4,
one box per chromosome, sorted by the $S_{C_{i}}$ medians. By reading the
chromosome numbers and the x-axis, we note no strong correlation between
chromosome size and stability, as both large and small chromosomes appear on
either end. To compare this result to experiments, we note that the same work
that calculated the median TAD boundary similarity to be 62% also used an
advanced similarity metric TADsim to measure the structural similarity of TAD
sets between replicates [16]. They found this overlap to be lower, about 50%
between replicates. This is close to our median stability between all
chromosomes, which turns out to be $\sim 0.4$.

Until this point, we studied each chromosome separately. Below, we aggregate
all the data and consider genome-wide averages. We also used $n_{\rm
samp}=1000$ bootstrapped samples per chromosome to calculate the stability,
which we bin to a resolution of 1%.

###  III.3 Understanding stability and community connectedness

Figure 5: The mean internal node strength $\hat{s}_{C_{i}}^{\rm int}$ versus
community stability $S_{C_{i}}$ at an effective community size of $\sim
0.88$Mb across all chromosomes. The mean internal node strength is binned by
the stability. We observe a nearly linear relationship between the two
quantities, indicating that the internal strength within a community is a good
predictor of community stability between bootstrapped samples. Figure 6:
Community stability versus folds of enrichment (FE) for several chromatin
groups (colored markers, defined in App. C). The number inside each panel
shows the Spearman correlation coefficient. We note a convincing positive
correlation with enrichment in the three leftmost panels (promoters,
enhancers, and transcribed regions). We note no or weak negative correlation
for the two rightmost panels (heterochromatin & repressive states,
insulators).

The boxplots and the distribution in Fig. 4 foreshadow that some of the
communities have a more resilient community structure than others under noise.
Here, we ask how much this stability depends on the local structure of the
underlying network. In particular, we study the communities’ internal node
strength. This metric differs from general node strength—the sum of all edge
weights of a node—by restricting the sum to nodes within a community. Node
strength is a good predictor of, for example, node centrality in a fully
connected network [38] or search times to a random node [39, 40, 41, 42].

However, in our case, the stability of a community is not related to how well
the nodes are connected to other parts of the network. Instead, we quantify
the internal community connectedness. We define this as the sum of all edge
weights within a community (including self-loops) normalized by the total node
strength

$\hat{s}_{C_{i}}^{\rm int}=\frac{\sum_{j\in C_{i}}\sum_{i\in
C_{i}}A_{ij}}{\sum_{j}\sum_{i\in C_{i}}A_{ij}},$ (11)

where ”int” stands for internal node strength. Plotting $\hat{s}_{C_{i}}^{\rm
int}$ for all communities across all chromosomes against the community
stability $S_{C_{i}}$ (Eq. (10)), we find that high internal node strength
suggests high community stability under noise (Fig. 5). This means that strong
inter-connected communities are less affected by random noise during
bootstrapping.

###  III.4 Relating community stability to chromatin states

The previous section studied how stability is connected to network structure.
Next, we perform a similar analysis using biologically relevant markers:
chromatin states. In particular, we ask if our stability metric calculated on
TAD-sized communities is associated with the enrichment or depletion of active
or inactive genomic regions, which could suggest underlying mechanisms that
form strongly connected communities. For example, there are several
indications of insulation enrichment at TAD borders and facilitated promoter-
enhancer interactions within TADs [43, 44].

To study this, we used data that classify genomic regions into 15 different
functional states (see Tab. 2), which we grouped into five categories (defined
in App. C). We plotted the folds of enrichment (FE) of these five groups
against community stability for TAD-sized communities ($\hat{s}\sim 0.38$ Mb)
in Fig. 6. For three of the groups associated with active genomic regions
(promoters, enhancers, and transcribed regions), there is a clear positive
relation between stability and FE. However, repressed regions (heterochromatin
& repressive states and insulators) show no such correlation. This indicates
that the active genomic regions are associated with strong stability towards
noise or communities having strong internal connectivity.

###  III.5 Quantifying cross-scale nestedness of stable communities

Thus far, we have focused on the stability of communities of a particular
size. However, by tuning the GenLouvain parameter $\gamma$ we can get
communities of varying sizes. This allows us to study how communities at one
scale nest with communities at larger scales. Such studies complement previous
work finding that active chromatin congregates in communities with a
significant hierarchical nestedness across community sizes $\hat{s}$ [15], in
particular, hierarchical TADs [45, 46]. Inspired by this, we ask: do nodes in
stable communities of one size stay within stable communities of a larger
size? Or is community nesting related to stability?

To this end, we first calculated the community stability $S_{C}$ for different
effective community sizes. We choose four commonly studied chromosome scales:
TADs ($\hat{s}\sim 0.38$MB), A/B segments ($\hat{s}\sim 1.04$MB), A1,2, B1,2,3
structures ($\hat{s}\sim 49$MB) and A/B compartments ($\hat{s}\sim 60$MB) (see
App. B). Next, we assign a node stability $s_{i}$ to each node $i$, which we
choose as the community stability $S_{C}$. Thus, as $\gamma$ varies, we obtain
$s_{i}$ values for each size $\hat{s}$. The change in $s_{i}$ reflects how
stability changes over different network scales.

We show the change in $s_{i}$ over three different scale transitions as a
boxplot in Fig. 7 (a) (green boxes). For relatively small communities (around
the size of TADs), we note that the stability difference is smaller in smaller
communities. To make a fair comparison, we calculated the same difference
between two randomly chosen nodes (orange boxes). These boxes show that the
stability difference is smaller in the actual data (green) than in the random
case, at least for smaller communities. This indicates that stable communities
at one size are generally supported by stable communities at a larger size,
indicating the cross-scale nestedness of stable communities.

However, this relationship holds best as communities become smaller. At large
sizes, the communities occupy a considerable part of the chromosome, implying
that randomization does not change much since the probability of selecting two
random nodes from the same community is high. This does not necessarily imply
that the stability is not nested at larger sizes but rather that a much larger
partition sample is required to get good statistics.

To better illustrate how nested stability changes between two select
chromosome scales ($\hat{s}=0.38\mathrm{Mb}\to\hat{s}=1.04\rm{Mb}$), we
plotted $s_{i}$ in an alluvial diagram (Fig. 7 (b)). To increase readability,
we rounded the $s_{i}$ values to a precision of 0.2. We note that most flows
from $\hat{s}=0.38\mathrm{Mb}$ to $\hat{s}=1.04\mathrm{Mb}$ go towards between
nodes with similar stabilities or just above or below. This gives additional
support to our conclusion that stability is nested across network scales. We
also depict the chromatin state enrichment FE as bars to the left and right of
the alluvial diagram (same data as in Fig. 6). We note that the FE levels for
active chromatin (three leftmost bars) remain conserved across both scales for
all six stability groups. The inactive chromatin shows similar conservation,
but this like stem from them being uncorrelated to node stability (Fig. 6).

Figure 7: Connection between stability, nestedness, and active genomic
regions. (a) Difference in node stability between different network scales
(community sizes). Each group in the boxplots shows the distribution of
absolute differences in node stability $s_{i}$ between two different scales.
The left boxes (green) show the difference from the actual data, and the right
(orange) boxes show $s_{i}$ after randomization. The x-labels indicate the
effective community sizes, and the vertical lines show the median. For smaller
communities, we note smaller $s_{i}$ in the actual data compared to the random
case. This indicates that the stability is hierarchical, where a stable node
in a smaller community belongs to a similarly stable community with a larger
size (or scale). However, this trend does not seem to hold for larger
communities. (b) Cross-scale stability and chromatin state. The alluvial
diagram shows the individual node stabilities between community sizes of 0.38
Mb and 1.04 Mb, rounded with a precision of $0.2$. We see that most of the
flows from left to right go from between similar stabilities (e.g., most of
the dark red at 0.38 Mb remains dark red at 1.04 Mb or slightly lighter). The
bars on the sides show the five coarse-grained chromatin states (same colors
as in Fig. 6, indicated by colors). We see that correlated HMM states stay
correlated between different community sizes, where groups 4 and 5 are notable
exceptions.

##  IV Discussion

Hi-C represents a promising method to uncover essential relationships between
gene activity and chromosome 3D organization. However, extracting biologically
meaningful 3D structures from Hi-C data sets poses several challenges. One
challenge concerns the 3D contacts forming a nearly fully connected complex
network with cross-scale interactions. This leads to inconsistent network
partitions when comparing data clustering or TAD-finding methods [18, 13, 20].
Another challenge comes from experimental noise, where the same clustering
method yields slightly different network partitions between theoretically
identical replicate Hi-C experiments. Here, we propose a new method that
systematically helps deal with the noise problem. We also propose a new null
model for the distant-dependent average contact frequency that builds on Hi-C
maps’ multifractal properties. Using our method, we find that stable 3D
communities typically have a high internal connectedness, tend to be enriched
in active chromatin marks, and have a more nested cross-scale hierarchy than
less stable ones.

As mentioned previously and discussed before [15], genome organization does
not likely follow a perfect hierarchical tree structure. This is partly due to
scale-dependent folding mechanisms, such as loop extrusion and chromatin-
chromatin interactions, which are critical in forming TADs and A/B
compartments [47]. Our work focuses on extracting stable communities under
noise and suggests new folding mechanisms that are worthy of further
experimental inquiry. We base this conclusion on Fig. 7(b), where highly
stable communities maintain chromatin state enrichment across chromosome
scales. This foreshadows the same mechanisms being responsible for
hierarchical folding in select parts of the 3D hierarchy.

We acknowledge that our results stem from a specific Hi-C data resolution
(100kb). However, our method is not limited to this choice and can handle any
resolution and Hi-C data set. In fact, we do not anticipate drastic changes to
our findings because previous research finds persistent correlations between
TAD boundaries and structure in the 40-100kb range [16]. When applying the
approach in this paper to a different Hi-C resolution, it is essential to
double-check the noise amplitude $\epsilon^{*}$, as it is not a universal
parameter.

Furthermore, $\epsilon^{*}$ also depends on the specific choice of community
detection method as we calibrate against the communities associated with the
unperturbed Hi-C data set. For example, if choosing one of the TAD-finding
methods in [18]. Our pipeline also allows experimenters to assess how stable
3D communities are under noise from varying sampling depths [19]. This case
likely requires a revised model of the noise (Sec. II.2) and subsequent
recalibration of $\epsilon^{*}$. Otherwise, the pipeline is the same.

We choose the GenLouvin algorithm to find communities. While it has known
shortcomings [48, 20], it is one of the most widely used methods in network
science because it is relatively simple and allows specifying a theoretical
null. Here, we use a null model that takes advantage of the multifractal
spectrum of the Hi-C map. We believe this approach is useful for community
detection in any network setting representing a complex system because it
provides a straightforward way of parameterizing an explicit null model that
otherwise, in lack of better options, often becomes the random Newman-Girvan
model. While there exist methods to construct multifractal networks (e.g.,
[49]), our approach addresses the inverse problem.

To summarize, our work expands the techniques employed for community detection
in Hi-C data and introduces additional methods and metrics to extract
communities that survive realistic noise levels. We anticipate having better
access to robust 3D communities will help future research uncover causal
connections between chromosome contact networks and genetic processes, such as
gene expression and repression.

##  V Author contributions

L.L. and L.H. devised the study. L.H. and A.C.M. performed the analysis. L.H.
made the visualizations and wrote the manuscript. All authors edited and
accepted the manuscript in its final form.

##  VI Competing interests

The authors declare no competing interests.

###### Acknowledgements.

We thank Moa Lundkvist and Juhee Lee for providing valuable discussion and
feedback. We also thank Moa Lundkvist for her help creating the alluvial plot
in Fig. 7. We acknowledge financial support from the Swedish Research Council
(grant no. 2017-03848). The computations were enabled by resources provided by
the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and
the Swedish National Infrastructure for Computing (SNIC) at High-Performance
Computing Center North (HPC2N), partially funded by the Swedish Research
Council through grant agreements no. 2022-06725 and no. 2018-05973.

##  Appendix A The hierarchical domain model as a null model

As discussed in Methods, the null model used for community detection has to be
chosen to represent the underlying structure of the data. Mimicking the power-
law dropoff in contact probability observed in Hi-C, the null model for Hi-C
is commonly chosen to follow $\sim d^{-\alpha}$, where $d=|i-j|$ is the
distance between two bins on Hi-C and $\alpha$ is a scale parameter [14].
However, later work showed that to capture communities at different scales on
Hi-C $\alpha$ had to be fine-tuned to avoid degenerate partitions, i.e.
partitions that are very different in structure and quality [20]. To remove
the need for fine-tuning a scale parameter, we look at recent work where the
contact dropoff in Hi-C was fitted to a bifractal contact map [24].

In that paper, they devise a hierarchical domain model (HDM) as a more
representative metric for the distance-dependent structure of Hi-C. In short,
the algorithm starts from a $2\times 2$ matrix, which self-repeats $n$ times
to a size $2^{n}$. We use the analytical expression of this contact map to
calculate the average contact probability dependent on the distance between
two Hi-C bins show in Eq. (II.3). Since the Hi-C matrix of size $N$ is not
necessarily a multiple of $2$, we let $H_{a}(i,j)$ repeat for $|i-j|\geq
2^{n}$, which is a constant value. We show $H_{a}(i,j)$ versus the genomic
distance, along with commonly used contact probabilities [14] in Fig. 8 (a).
By capturing the different fractal scales of the chromosomes using this null
model, we remove the need to fine-tune the scale parameter $\alpha$ to capture
communities at different sizes. Instead, we only fine-tuned the parameter $a$
to match our Hi-C data, which we calculated to have the mean $a=0.464(3)$.

Figure 8: Comparing the HDM null-model versus two different power-law
exponents. (a) The HDM model shows a similar power-law scaling compared to the
slope using $\alpha=0.75$. However, the HDM has a slight change in curvature,
flattening out for large distances. (b) Degeneracy evaluation plot recreated
with permission from [20] using chromosome 10. Two metrics versus the
effective community size are shown, the mean distance ${\rm MD}(d)$ between
pairwise partitions for different community detection runs, and the
coefficient of variation ${\rm CV}(Q)$ of the modularity function between
partitions. Higher values in both metrics indicate separated partitions with
differing qualities, which are degenerate states. We see that the HDM avoids
local minima by having a low ${\rm CV}(Q)$, albeit with slightly larger
distances between partitions for smaller communities.

To evaluate the HDM null model against the power-law scaling based on the
fractal globule assumption, we recreate some of the results in [20], shown in
Fig. 8 (b). In short, the metric ${\rm MD}(d)$ quantifies the distance between
partitions, and a large value indicates that the partitions between community
detection runs are structurally significantly different. ${\rm CV}(Q)$
quantifies the coefficient-of-variation between partitions, where a large
value indicates that the partition varies much in quality, indicating strong
local minima. We see that ${\rm MD}(d)$ for the HDM null model falls somewhere
in between the two power-law relationships. However, the HDM has on average a
smaller variation in quality $Q$, indicating that even if the partitions have
a slightly larger distance between each other, they are comparatively similar
in quality, indicating that there is little risk of getting stuck in local
minima.

##  Appendix B Picking $\gamma$ for different sized communities on Hi-C

We picked the parameter $\gamma$ in the community detection pipeline to
approximately represent structural elements at different scales on the Hi-C.
We calculated the median effective community size of each partition for all
chromosomes as in [20] and compared them against their values reported in the
same paper. The parameter $\gamma$ and the corresponding structural scales are
shown in Table. 1.

Table 1:  The value of $\gamma$ that reproduce certain scaled communities on the Hi-C. The sizes of the structures and how to calculate the characteristic community size can be found in [20]. Here we used 100kb data using $a=0.464$ in the hierarchical domain model. Structure | Size (Mb) | $\gamma$ | Comm. size (Mb)
---|---|---|---
TADs | $0.33$ | $1.75$ | $0.38$
TADs | $0.88$ | $1.4$ | $0.90$
A/B segments | $1.1$ | $1.35$ | $1.04$
A1,…,B3 | $27$ | $0.55$ | $49$
A/B compartments | $66$ | $0.40$ | $60$

##  Appendix C Preprocessing chromatin states to Hi-C bins

To compare our results of community stability on Hi-C with different
functional regions on the genome, we used data dividing the genome into 15
separate functional divisions [32, 33], taken from the ENCODE database [31].
Due to the functional similarity between different states, we group the states
into five groups. The different states, along with their names and groupings
can be found in Table. 2 below.

Table 2:  The 15 functional states identified on the genome grouped by their function [33, 32, 50]. Here we show the state numbers, their functional names and the groups which we combine the states into. Group & name | HMM state | Name
---|---|---
(A) Promoters | S1 | active promoter
S2 | weak promoter
(B) Enhancers | S4 | strong enhancer
S5 | strong enhancer
S6 | weak/poised enhancer
S7 | weak/poised enhancer
(C) Transcribed regions | S9 | transcriptional transition
S10 | transcriptional elongation
S11 | weakly transcribed
(D) Heterochromatin & repressive states | S3 | inactive/poised promoter
S13 | heterochromatin
S14 | repetitive/copy number variation
S15 | repetitive/copy number variation
(E) Insulators | S8 | insulator

The dataset contains start and stop regions for each of the state. To map how
much of each state is expressed in each Hi-C bin we calculate the folds of
enrichment (FE) relative to the genome-wide average following the same
procedure as in [50]

  1. 1.

We count the number of peaks $k_{X}$ per bin, where $X$ corresponds to one of
the states. We only count the start of the peaks.

  2. 2.

We calculate the peak’s expected frequency value using the hypergeometric test
as $k^{\prime}_{X}=K_{X}\times(n/N)$. Here $n$ is the number of peaks of any
state in a bin, $N$ is the total number of peaks per chromosome, and $K_{X}$
is the total number of peaks for the state $X$.

  3. 3.

Using this, we calculate the folds of enrichment for each state $Z$ by
dividing the observed peak by the expected $FE_{X}=k_{X}/k^{\prime}_{X}$.

## References

  * Dixon  _et al._ [2016] J. R. Dixon, D. U. Gorkin, and B. Ren, Molecular cell 62, 668 (2016).
  * Schwartz and Cavalli [2017] Y. B. Schwartz and G. Cavalli, Genetics 205, 5 (2017).
  * Bonev and Cavalli [2016] B. Bonev and G. Cavalli, Nature Reviews Genetics 17, 661 (2016).
  * Denker and De Laat [2016] A. Denker and W. De Laat, Genes & development 30, 1357 (2016).
  * Marchal  _et al._ [2019] C. Marchal, J. Sima, and D. M. Gilbert, Nature Reviews Molecular Cell Biology 20, 721 (2019).
  * Lieberman-Aiden  _et al._ [2009] 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, _et al._ , Science 326, 289 (2009).
  * Sexton  _et al._ [2012] T. Sexton, E. Yaffe, E. Kenigsberg, F. Bantignies, B. Leblanc, M. Hoichman, H. Parrinello, A. Tanay, and G. Cavalli, Cell 148, 458 (2012).
  * Dekker  _et al._ [2013] J. Dekker, M. A. Marti-Renom, and L. A. Mirny, Nature Reviews Genetics 14, 390 (2013).
  * MacKay and Kusalik [2020] K. MacKay and A. Kusalik, Briefings in functional genomics 19, 292 (2020).
  * Liu  _et al._ [2021] Y. Liu, L. Nanni, S. Sungalee, M. Zufferey, D. Tavernari, M. Mina, S. Ceri, E. Oricchio, and G. Ciriello, Nature communications 12, 1 (2021).
  * Weinreb and Raphael [2016] C. Weinreb and B. J. Raphael, Bioinformatics 32, 1601 (2016).
  * Fraser  _et al._ [2015] J. Fraser, C. Ferrai, A. M. Chiariello, M. Schueler, T. Rito, G. Laudanno, M. Barbieri, B. L. Moore, D. C. Kraemer, S. Aitken, _et al._ , Molecular systems biology 11, 852 (2015).
  * Cabreros  _et al._ [2016] I. Cabreros, E. Abbe, and A. Tsirigos, in  _2016 Annual Conference on Information Science and Systems (CISS)_ (IEEE, 2016) pp. 584–589.
  * Lee  _et al._ [2019] S. H. Lee, Y. Kim, S. Lee, X. Durang, P. Stenberg, J.-H. Jeon, and L. Lizana, Scientific reports 9, 1 (2019).
  * Bernenko  _et al._ [2022] D. Bernenko, S. H. Lee, P. Stenberg, and L. Lizana, bioRxiv (2022).
  * Sauerwald  _et al._ [2020] N. Sauerwald, A. Singhal, and C. Kingsford, NAR genomics and bioinformatics 2, lqz008 (2020).
  * de Wit [2020] E. de Wit, Journal of molecular biology 432, 638 (2020).
  * Sefer [2022] E. Sefer, BMC bioinformatics 23, 127 (2022).
  * Yardımcı  _et al._ [2019] G. G. Yardımcı, H. Ozadam, M. E. Sauria, O. Ursu, K.-K. Yan, T. Yang, A. Chakraborty, A. Kaul, B. R. Lajoie, F. Song, _et al._ , Genome biology 20, 1 (2019).
  * Holmgren  _et al._ [2023] A. Holmgren, D. Bernenko, and L. Lizana, Scientific Reports 13, 12979 (2023).
  * Rao  _et al._ [2014] S. S. 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, _et al._ , Cell 159, 1665 (2014).
  * Eres and Gilad [2021] I. E. Eres and Y. Gilad, Trends in Genetics 37, 216 (2021).
  * Rosvall and Bergstrom [2010] M. Rosvall and C. T. Bergstrom, PloS one 5, e8694 (2010).
  * Pigolotti  _et al._ [2020] S. Pigolotti, M. H. Jensen, Y. Zhan, and G. Tiana, Physical Review Research 2, 043078 (2020).
  * Sarnataro  _et al._ [2017] S. Sarnataro, A. M. Chiariello, A. Esposito, A. Prisco, and M. Nicodemi, PLoS One 12, e0188201 (2017).
  * Newman and Girvan [2004] M. E. Newman and M. Girvan, Physical review E 69, 026113 (2004).
  * Yan  _et al._ [2017] K.-K. Yan, S. Lou, and M. Gerstein, PLoS computational biology 13, e1005647 (2017).
  * Sanborn  _et al._ [2015] A. L. Sanborn, S. S. Rao, S.-C. Huang, N. C. Durand, M. H. Huntley, A. I. Jewett, I. D. Bochkov, D. Chinnappan, A. Cutkosky, J. Li, _et al._ , Proceedings of the National Academy of Sciences 112, E6456 (2015).
  * Edgar  _et al._ [2002] R. Edgar, M. Domrachev, and A. E. Lash, Nucleic acids research 30, 207 (2002).
  * Knight and Ruiz [2013] P. A. Knight and D. Ruiz, IMA Journal of Numerical Analysis 33, 1029 (2013).
  * Raney  _et al._ [2024] B. J. Raney, G. P. Barber, A. Benet-Pagès, J. Casper, H. Clawson, M. S. Cline, M. Diekhans, C. Fischer, J. Navarro Gonzalez, G. Hickey, _et al._ , Nucleic Acids Research 52, D1082 (2024).
  * Ernst  _et al._ [2011] J. Ernst, P. Kheradpour, T. S. Mikkelsen, N. Shoresh, L. D. Ward, C. B. Epstein, X. Zhang, L. Wang, R. Issner, M. Coyne, _et al._ , Nature 473, 43 (2011).
  * Ernst and Kellis [2010] J. Ernst and M. Kellis, Nature biotechnology 28, 817 (2010).
  * Lajoie  _et al._ [2015] B. R. Lajoie, J. Dekker, and N. Kaplan, Methods 72, 65 (2015).
  * Jeub  _et al._ [2019] L. Jeub, M. Bazzi, I. Jutla, and P. Mucha, https://github.com/GenLouvain/GenLouvain (2011-2019).
  * Grosberg  _et al._ [1988] A. Y. Grosberg, S. K. Nechaev, and E. I. Shakhnovich, Journal de physique 49, 2095 (1988).
  * Ghosh and Jost [2018] S. K. Ghosh and D. Jost, PLoS computational biology 14, e1006159 (2018).
  * Opsahl  _et al._ [2010] T. Opsahl, F. Agneessens, and J. Skvoretz, Social networks 32, 245 (2010).
  * Nyberg  _et al._ [2021] M. Nyberg, T. Ambjörnsson, P. Stenberg, and L. Lizana, Physical Review Research 3, 013055 (2021).
  * Hedström and Lizana [2023] L. Hedström and L. Lizana, New Journal of Physics 25, 033024 (2023).
  * Noh and Rieger [2004] J. D. Noh and H. Rieger, Physical review letters 92, 118701 (2004).
  * Tejedor  _et al._ [2009] V. Tejedor, O. Bénichou, and R. Voituriez, Physical Review E 80, 065104 (2009).
  * Gong  _et al._ [2018] Y. Gong, C. Lazaris, T. Sakellaropoulos, A. Lozano, P. Kambadur, P. Ntziachristos, I. Aifantis, and A. Tsirigos, Nature communications 9, 542 (2018).
  * Qu  _et al._ [2019] J. Qu, G. Yi, and H. Zhou, Epigenetics & chromatin 12, 31 (2019).
  * Liu  _et al._ [2022] E. Liu, H. Lyu, Q. Peng, Y. Liu, T. Wang, and J. Han, Communications Biology 5, 608 (2022).
  * Cresswell  _et al._ [2020] K. G. Cresswell, J. C. Stansfield, and M. G. Dozmorov, BMC bioinformatics 21, 1 (2020).
  * Nuebler  _et al._ [2018] J. Nuebler, G. Fudenberg, M. Imakaev, N. Abdennur, and L. A. Mirny, Proceedings of the National Academy of Sciences 115, E6697 (2018).
  * Fortunato and Hric [2016] S. Fortunato and D. Hric, Physics reports 659, 1 (2016).
  * Palla  _et al._ [2010] G. Palla, L. Lovász, and T. Vicsek, Proceedings of the National Academy of Sciences 107, 7640 (2010).
  * Bernenko  _et al._ [2023] D. Bernenko, S. H. Lee, and L. Lizana, Journal of Physics: Complexity 4, 035004 (2023).