# Multiplicity and Transverse Spherocity dependence of $\langle p_{\rm
T}\rangle$ fluctuations of charged particles in p$-$p collisions at $\sqrt{s}$
= 7 and 13 TeV

Subhadeep Roy<EMAIL_ADDRESS>Indian Institute of Technology Bombay,
Mumbai 400076, India  Tulika Tripathy<EMAIL_ADDRESS>Indian
Institute of Technology Bombay, Mumbai 400076, India  Sadhana Dash
<EMAIL_ADDRESS>Indian Institute of Technology Bombay, Mumbai 400076,
India

###### Abstract

The multiplicity dependence of event-by-event fluctuations in mean transverse
momentum, $\langle p_{\rm T}\rangle$, of charged particles has been studied in
p$-$p collisions at $\sqrt{s}$ = 7 TeV and 13 TeV using the PYTHIA 8 event
generator. The charged particles were selected in kinematic range of
$0.15<p_{\rm T}<2$ GeV/c and $|\eta|<0.8$. The dynamical fluctuations would
indicate towards the correlated emission of particles. The measurements in
A$-$A and p$-$p collisions has shown a decrease in the strength of $\langle
p_{\rm T}\rangle$ fluctuations with the average charged particle multiplicity.
The effects of various microscopic processes like color reconnection and
multi-partonic interactions has been studied. A minimal dependency on the
collision energy is also observed. Furthermore, the fluctuation observables
are investigated in the intervals of transverse spherocity in order to
comprehend the relative contributions resulting from hard scattering and
underlying events. The present study would act as a baseline for future
measurements in A$-$A as well as p$-$p collisions at the LHC.

##  I Introduction

The search for event-by-event fluctuations in dynamical quantities, such as
mean transverse momentum ($\langle p_{\rm T}\rangle$) in heavy-ion collisions
is primarily motivated by the search for the presence and characterization of
the phase transition between a quark-gluon plasma (QGP) and hadron gas state
2, 3, 4. One can quantify this search by looking for excess fluctuations of
thermodynamic quantities like temperature. The event-wise $\langle p_{\rm
T}\rangle$ can be taken as a proxy for local temperature and hence the study
of its event-by-event fluctuation has a thermodynamic context. The temperature
fluctuations are related to the heat capacity of the system. The non-monotonic
behaviour of heat capacity characterizes a phase transition and hence the
study of event-by-event $\langle p_{\rm T}\rangle$ fluctuations could reveal
important information about the formation of QGP 5.

The recent measurements by the STAR 2, 6 and ALICE 3 experiment at RHIC and
LHC reported that the strength of dynamical $\langle p_{\rm T}\rangle$
fluctuations decreased with event multiplicity (or centrality) in heavy-ion
collisions. The observation was attributed to the onset of thermalization,
collectivity, jet suppression, presence of di-jets and minijets and other
processes. The interesting studies showing the $\langle p_{\rm T}\rangle$
fluctuations and their correlations with radial flow and elliptic flow in
heavy-ion collisions can be found in the references 7, 8. It remains
questionable whether the system created in small systems like p$-$p or p$-$A
collisions exhibit collectivity as observed in A$-$A collisions. This can be
addressed by studying particle production mechanisms, their correlations etc.
as a function of the event shape and particle multiplicity. Recent
experimental results of $\langle p_{\rm T}\rangle$ in p$-$p collisions 3
suggest that, high-multiplicity events are mostly produced by multiple parton
interactions (MPIs) 9. The data exhibited a power-law trend with charged
particle multiplicity density, as expected from the independent superposition
scenario of independent multi-partonic interactions. However, the slight
deviation of the power-law index from $0.5$ indicated towards non-trivial
fluctuations. In PYTHIA 8 10, the observed strong correlation between the
particle multiplicity and $\langle p_{\rm T}\rangle$ has been attributed to
the mechanism of color reconnections (CR) between hadronizing strings 12, 13,
14. The color strings can also overlap with each other to form color ropes
that act coherently in high-multiplicity events 15, 16. This overlap region
causes a pressure gradient, thereby shoving the strings in the outward
direction creating similar effects as observed in hydrodynamic scenario. The
color reconnections together with the rope hadronization could explain some of
the high-multiplicity observables in p$-$p collisions17, 18, 19. A similar
mechanism of collective relativistic string hadronization is also implemented
in the EPOS model which describes a wealth of LHC data in p$-$p, p$-$Pb, and
Pb$-$Pb collisions 20.

In this work, an attempt has been made to study the $\langle p_{\rm T}\rangle$
fluctuations in p$-$p collisions at collision energies, $\sqrt{s}$ = 7 TeV and
13 TeV as a function of multiplicity and different tansverse spherocity
classes using PYTHIA 8 generator. The fluctuation is generally quantified by
the the second moment of the distribution and the two-particle correlator has
been used to study the interplay of various processes. The higher moments,
namely the third and fourth moment were also studied as a function of
multiplicity to shed additional lights on higher order correlations in multi-
particle production mechanism. The additional effects of color reconnection as
well as the formation of ropes were also studied. This study would act as a
baseline for future measurement in A$-$A and p$-$p systems.

##  II Transverse Spherocity

The study of various observables in terms of event shape variables such as
transverse spherocity, thurst etc. helps to disentangle the contributions from
hard and soft (scattering) processes of particle production. The transverse
spherocity is a tool which characterises the event into jet-like or isotropic
events 21, 22, 23 . It is defined for a unit transverse vector ,
${\hat{n}(n_{\rm T},0)}$, which minimizes the ratio

$S_{0}=\frac{\pi^{2}}{4}\underset{\hat{n}}{\rm
min}\left({\frac{\sum\limits_{i}|\overrightarrow{p_{\rm
T_{i}}}\times\hat{n}|}{\sum\limits_{i}p_{\rm T_{i}}}}\right)^{2}$ (1)

This unit vector coincides with the transverse direction of the momentum,
$\overrightarrow{p_{\rm T}}$ which makes the $S_{0}$ infrared and collinear
safe. The sum runs over all the particles considered in the kinematic
acceptance. Events with $S_{\rm 0}\rightarrow 0$ are likely to be dominated by
a single hard scattering, producing a jet-like structure (Jetty). In contrary,
$S_{\rm 0}\rightarrow 1$ implies that particles are produced through several
softer interactions, giving an isotropic distribution. An illustration of the
two limits in the azimuthal plane is presented in Figure 1.

Figure 1: Illustration of the typical topologies for the two $S_{0}$ limits in
the azimuthal plane.

To replicate the same circumstances as the ALICE experiment at the LHC, the
transverse spherocity distributions are chosen in the pseudo-rapidity range of
$|\eta|<$ 0.8 with a minimum constraint of $5$ charged particles with $p_{\rm
T}>$ 0.15 GeV/c. The various percentiles of the whole $S_{0}$ distribution are
used for the selection of Jetty (0-20%) and isotropic (80-100%) event classes.
Transverse spherocity will henceforth be referred to as spherocity for the
sake of simplicity.

## Observables

The mean value of the transverse momentum, obtained for a particular
multiplicity class can be defined as

$\left\langle\\!\left\langle p_{\rm
T}\right\rangle\\!\right\rangle=\left\langle\frac{\sum_{i=1}^{N_{\rm
ch}}p_{{\rm i}}}{N_{\rm ch}}\right\rangle$ (2)

where $N_{\rm ch}$ denotes the number of charged particles in a single event
and $p_{{\rm i}}$ represents the transverse momentum of the $i^{\rm th}$
particle in the event. The first average is determined over all of the
$p_{{\rm i}}$ in an event, while the second average is carried out over all of
the events included in a certain multiplicity class. From this point on, the
mean transverse momentum of the particles will be denoted by
$\left\langle\\!\left\langle p_{\rm T}\right\rangle\\!\right\rangle$.

The two-particle transverse momentum correlator, $\left\langle\Delta p_{\rm
i}\Delta p_{\rm j}\right\rangle$, is constructed as $-$

$\left\langle\Delta p_{\rm i}\Delta p_{\rm
j}\right\rangle=\left\langle\frac{\sum_{i,j\neq i}(p_{\rm
i}-\left\langle\left\langle p_{\rm T}\right\rangle\right\rangle)(p_{\rm
j}-\left\langle\left\langle p_{\rm T}\right\rangle\right\rangle)}{N_{\rm
ch}(N_{\rm ch}-1)}\right\rangle$ (3)

Similarly, the three-particle and four-particle correlators are given by

$\left\langle\Delta p_{\rm i}\Delta p_{\rm j}\Delta p_{\rm
k}\right\rangle=\left\langle\frac{\sum_{i,j\neq i,k\neq i,j}(p_{\rm
i}-\left\langle\left\langle p_{\rm T}\right\rangle\right\rangle)(p_{\rm
j}-\left\langle\left\langle p_{\rm T}\right\rangle\right\rangle)(p_{\rm
k}-\left\langle\left\langle p_{\rm T}\right\rangle\right\rangle)}{N_{\rm
ch}(N_{\rm ch}-1)(N_{\rm ch}-2)}\right\rangle$ (4) $\left\langle\Delta p_{\rm
i}\Delta p_{\rm j}\Delta p_{\rm k}\Delta p_{\rm
l}\right\rangle=\left\langle\frac{\sum_{i,j\neq i,k\neq i,j,l\neq
i,j,k}(p_{\rm i}-\left\langle\left\langle p_{\rm
T}\right\rangle\right\rangle)(p_{\rm j}-\left\langle\left\langle p_{\rm
T}\right\rangle\right\rangle)(p_{\rm k}-\left\langle\left\langle p_{\rm
T}\right\rangle\right\rangle)(p_{\rm l}-\left\langle\left\langle p_{\rm
T}\right\rangle\right\rangle)}{N_{\rm ch}(N_{\rm ch}-1)(N_{\rm ch}-2)(N_{\rm
ch}-3)}\right\rangle$ (5)

These correlators can be written in a more simpler form when expressed in
terms of raw moments as follows 24 :

$\displaystyle\langle\\!\langle p_{\rm T}\rangle\\!\rangle$
$\displaystyle=\left\langle\frac{Q_{1}}{N_{\rm ch}}\right\rangle,$
$\displaystyle\left\langle\Delta p_{\rm i}\Delta p_{\rm j}\right\rangle$
$\displaystyle=\left\langle\frac{Q_{1}^{2}-Q_{2}}{N_{\rm ch}\left(N_{\rm
ch}-1\right)}\right\rangle-\left\langle\frac{Q_{1}}{N_{\rm
ch}}\right\rangle^{2},$ $\displaystyle\left\langle\Delta p_{\rm i}\Delta
p_{\rm j}\Delta p_{\rm k}\right\rangle$
$\displaystyle=\left\langle\frac{Q_{1}^{3}-3Q_{2}Q_{1}+2Q_{3}}{N_{\rm
ch}\left(N_{\rm ch}-1\right)\left(N_{\rm
ch}-2\right)}\right\rangle-3\left\langle\frac{Q_{1}^{2}-Q_{2}}{N_{\rm
ch}\left(N_{\rm ch}-1\right)}\right\rangle\left\langle\frac{Q_{1}}{N_{\rm
ch}}\right\rangle+2\left\langle\frac{Q_{1}}{N_{\rm ch}}\right\rangle^{3}.$
$\left\langle\Delta p_{\rm i}\Delta p_{\rm j}\Delta p_{\rm k}\Delta p_{\rm
l}\right\rangle=\left\langle\frac{Q_{1}^{4}-6Q_{4}+8Q_{1}Q_{3}-6Q_{1}^{2}Q_{2}+3Q_{2}^{2}}{N_{\rm
ch}(N_{\rm
ch}-1)(N_{ch}-2)(N_{ch}-3)}\right\rangle-4\left\langle\frac{Q_{1}^{3}-3Q_{2}Q_{1}+2Q_{3}}{N_{\rm
ch}\left(N_{\rm ch}-1\right)\left(N_{\rm
ch}-2\right)}\right\rangle\left\langle\frac{Q_{1}}{N_{\rm ch}}\right\rangle\\\
+6\left\langle\frac{Q_{1}^{2}-Q_{2}}{N_{\rm ch}\left(N_{\rm
ch}-1\right)}\right\rangle\left\langle\frac{Q_{1}}{N_{\rm
ch}}\right\rangle^{2}-3\left\langle\frac{Q_{1}}{N_{\rm ch}}\right\rangle^{3}.$
(6)

where, $Q_{\rm n}$s are the moments of the transverse momentum distributions
in an event and is written as

$Q_{\rm n}=\sum_{i=1}^{N_{\rm ch}}(p_{\rm i})^{n},$ (7)

$p_{\rm i}$ denotes the transverse momentum of the particle $i$, and the sum
runs over all the charged particles in an event. The different value of $n$
corresponds to different order of the moment.

In order to make the analysis less sensitive to the range in $p_{\rm i}$, we
normalise the parameters with $\left\langle\\!\left\langle p_{\rm
T}\right\rangle\\!\right\rangle$. Thus, the final expressions for two-particle
and four-particle correlator under the study takes the following forms $-$

$\displaystyle\left\langle\Delta p_{\rm i}\Delta p_{\rm
j}\right\rangle_{final}=\frac{\sqrt{\left\langle\Delta p_{\rm i}\Delta p_{\rm
j}\right\rangle}}{\left\langle\\!\left\langle p_{\rm
T}\right\rangle\\!\right\rangle}$ (8) $\displaystyle\left\langle\Delta p_{\rm
i}\Delta p_{\rm j}\Delta p_{\rm k}\Delta p_{\rm
l}\right\rangle_{final}=\frac{(\left\langle\Delta p_{\rm i}\Delta p_{\rm
j}\Delta p_{\rm k}\Delta p_{\rm
l}\right\rangle)^{1/4}}{\left\langle\\!\left\langle p_{\rm
T}\right\rangle\\!\right\rangle}$ (9)

Likewise, in order to make the three-particle correlator dimensionless, it is
normalised as follows

$\displaystyle\gamma_{p_{\rm T}}=\frac{\left\langle\Delta p_{\rm i}\Delta
p_{\rm j}\Delta p_{\rm k}\right\rangle}{{\left\langle\Delta p_{\rm i}\Delta
p_{\rm j}\right\rangle}^{3/2}}$ (10)

Now the obtained version of the three-particle correlator is called
standardized skewness , and further has a system size or centrality
dependence, which is measured by the number of participant nucleons in the
collision process. So, in an effort to eliminate the trivial size dependence,
a second measure of the skewness can be obtained.

$\displaystyle\Gamma_{p_{\rm T}}=\frac{\left\langle\Delta p_{\rm i}\Delta
p_{\rm j}\Delta p_{\rm k}\right\rangle\left\langle\\!\left\langle p_{\rm
T}\right\rangle\\!\right\rangle}{{\left\langle\Delta p_{\rm i}\Delta p_{\rm
j}\right\rangle}^{2}}$ (11)

$\Gamma_{p_{\rm T}}$ is called intensive skewness and is independent of the
number of participant nucleons or centrality.

##  III Results and Discussion

The present analysis has been carried out with 100 million events, generated
for p$-$p collisions at $\sqrt{s}=$ 7 TeV and 13 TeV using PYTHIA 8 event
generator with the standard Monash 2013 tune. The charged particles were
accepted in a pseudo-rapidity window of $|\eta|<0.8$. Further, charged
particles having transverse momentum range $0.15<p_{\rm T}<2.0$ GeV/c were
considered for this investigation. The event multiplicity estimation was done
by classifying the events based on the charged particle multiplicity obtained
in the forward (backward) pseudo-rapidity ranges : $2.8<\eta<5.1$ and
$-3.7<\eta<-1.7$. The obtained multiplicity values scale linearly with the
charged particle multiplicity in the accepted range ( $|\eta|<0.8$), used for
the analysis. For each event class, the corresponding mean charged
multiplicity, $\langle N_{\rm ch}\rangle$ was obtained in $|\eta|<0.8$. This
was done to avoid auto-correlation biases.

Figure 2: The top panel shows the variation of $\langle\\!\langle p_{\rm
T}\rangle\\!\rangle$ with $\langle N_{\rm ch}\rangle$ for p$-$p collisions at
$\sqrt{s}=$ 7 TeV and 13 TeV. The effect of different modes of color
reconnection (CR) on the variation of $\langle\\!\langle p_{\rm
T}\rangle\\!\rangle$ with multiplicity for collision energy $\sqrt{s}=$ 13 TeV
is shown in the middle panel. The bottom panel shows the variation of
$\langle\\!\langle p_{\rm T}\rangle\\!\rangle$ with $\langle N_{\rm
ch}\rangle$ for p$-$p collisions at $\sqrt{s}=$ 13 TeV for different
spherocity classes.

Figure 2 shows the variation of $\langle\\!\langle p_{\rm T}\rangle\\!\rangle$
as a function of mean charged particle multiplicity $\langle N_{ch}\rangle$
within $|\eta|<0.8$ for p$-$p collisions at $\sqrt{s}=$ 7 TeV and 13 TeV. The
effect of different modes of color reconnections has been shown in the middle
panel of Figure 2 for 13 TeV beam energy. One can observe that the
$\langle\\!\langle p_{\rm T}\rangle\\!\rangle$ values are systematically
higher for higher collison energy. It can also be observed that switching on
the color reconnection (CR) modes increases the $\langle\\!\langle p_{\rm
T}\rangle\\!\rangle$ value of the charged particles. The implementation of CR
models in PYTHIA 8 was found to mimic several collective-like behaviour seen
in heavy-ion collisons and thus, one is prompted to study the effects of
different CR models on the fluctuation observables. The mean transverse
momentum, $\langle\\!\langle p_{\rm T}\rangle\\!\rangle$ of the charged
particles shows a clear rise with the average charged particle multiplicity
($\langle N_{\rm ch}\rangle$), for all CR modes, compared to the case when
colour reconnection is turned off (middle panel). The QCD-based CR model
(CR=1) with the rope formation exhibits the strongest upward trend. It should
be mentioned that the mechanism of color reconnection was originally
introduced to explain the increase of $\langle\\!\langle p_{\rm
T}\rangle\\!\rangle$ with charged particle multiplicity. It allows the
partonic interactions through color strings from different semi-hard
scatterings. There is a prior fusion of strings from different multi-partonic
systems before hadronization, which leads to a decrease in overall
multiplicity but the particles are more energetic and consequently have a
higher value of $\langle\\!\langle p_{\rm T}\rangle\\!\rangle$.

The bottom plot in the Figure 2 shows the variation of $\langle\\!\langle
p_{\rm T}\rangle\\!\rangle$ as a function of average charged particle
multiplicity for different spherocity classes. One can observe that the values
are significantly higher for the events representing the lower 0-20%
spherocity class (Jetty). These are the events dominated by hard scatterings
and give rise to di-jet topology. The values for isotropic events (80-100%
spherocity class) are consistently lower than those for spherocity integrated
(0-100% spherocity) event class. The behaviour is qualitatively similar to the
one measured by ALICE experiment 11. This can be attributed to the dominance
of hard scattering events in the lower 0-20% spherocity class, while the
80-100% are dominated by underlying events.

As mentioned earlier, the event-by-event mean transverse momentum fluctuations
can be quantified using two-, three- and four- particle correlators. These
fluctuations mainly arises from the initial state fluctuations and contain
relevant information about the underlying physical processes of particle
production 24.

Figure 3:  The variation of two-particle correlator at $\sqrt{s}=$ 7 and 13
TeV as a function of $\langle N_{\rm ch}\rangle$. The middle panel shows the
effect of different color reconnection models at $\sqrt{s}=$13 TeV. The bottom
panel shows the variation of two-particle correlator with $\langle N_{\rm
ch}\rangle$ in different spherocity classes.

Figure 3 depicts the fluctuation of the two-particle correlator for p$-$p
collisions at $\sqrt{s}=$ 7 and 13 TeV as a function of charged particle
multiplicity. It can be observed that the strength of the two-particle
correlator decreases with an increase of the average multiplicity. This can be
attributed to the independent superposition of multi-partonic interactions in
p$-$p collisions. However, the values at $\sqrt{s}=$13 TeV are systematically
higher than that observed at 7 TeV, indicating an energy dependence of the
observable. The middle panel of the Figure 3 shows the effect of various modes
of color reconnections on the strength of the correlator. The mechanism of
color reconnection increases the correlation strength due to enhanced partonic
interactions via color strings. The interplay of rope formation in different
CR modes also causes a difference in the strength of the correlator. The
variation of the two-particle correlator as a function of multiplicity was
also compared for three different spherocity classes as shown in Figure 3. One
can observe that the values are higher for the lower 0-20% spherocity class
(Jetty), indicating that significant contribution comes from the di-jets in
low spherocity classes.

Recently, it was shown in the reference 24 that event-by-event fluctuations of
the $\langle\\!\langle p_{\rm T}\rangle\\!\rangle$ in heavy-ion collisions
exhibit a positive skew, as predicted by hydrodynamic evolution of the medium.
The idea stemmed from the strong correlation of transverse momentum
fluctuations with the fluctuations of the initial energy of the fluid. An
estimation of the three-particle correlation using standardized and intensive
skewness for p$-$p collisions at $\sqrt{s}=7$ and 13 TeV is shown in the
Figure 4 and 5, respectively. Assuming that the dynamical and statistical
fluctuations scale inversely with multiplicity 24, the standardized skewness
is expected to be larger for lower multiplicity classes. It is indeed observed
that $\gamma_{p_{\rm T}}$ exhibits a sharp diminishing trend upto a charged
particle multiplicity range of 15 and saturates thereafter. However, the
magnitude is considerably higher when compared to hydrodynamic predictions for
heavy-ion collisions. This trend is also seen for intensive skewness,
$\Gamma_{p_{\rm T}}$ which was introduced to eliminate the volume dependence
in heavy-ion collisions. Although, p$-$p system is not affected by volume
fluctuations, studying the variation of $\Gamma_{p_{\rm T}}$ is interesting to
form a baseline. Like $\gamma_{p_{\rm T}}$, the values of $\Gamma_{p_{\rm T}}$
are higher for low multiplicity events and gradually saturates after 15. In
addition, the comparisons of various CR models (middle panels of Fig. 4, 5)
revealed a decrease in the correlation strength when CR was introduced. This
could indicate that the color string fusion dominantly affects (or enhances)
the two-particle correlation and does not have any significant effect on
higher order correlations. The spherocity comparisons of the three-particle
correlator indicates higher correlation for jetty events (0-20% spherocity
class), whilst the isotropic events (80-100% spherocity class) have lower
values. (Fig. 4 and 5).

Figure 4:  The energy dependence (top) and the effect of different color
reconnection (CR) model comparison (middle) of $\gamma_{p_{\rm T}}$ in p$-$p
collisions at $\sqrt{s}=$ 7 TeV and 13 TeV. The bottom panel shows its
variation in different spherocity classes at $\sqrt{s}=$ 13 TeV.

Figure 5:  The energy dependence (top) and the effect of different color
reconnection (CR) model comparison (middle) of $\Gamma_{p_{\rm T}}$ in p$-$p
collisions at $\sqrt{s}=$ 7 TeV and 13 TeV. The bottom panel shows its
variation in different spherocity classes at $\sqrt{s}=$13 TeV.

Figure 6:  The four-particle correlator as a function of $\langle N_{\rm
ch}\rangle$ for different energies (top) and color reconnection (CR) models
(middle). The bottom panel shows the correlator strength in different
spheocity classes at $\sqrt{s}=13$ TeV.

The study was further extended to observe the variation of the four-particle
correlator with mean charged particle multiplicity for p$-$p collisions at
$\sqrt{s}=7$ and 13 TeV as shown in Figure 6. These types of studies are
subjected to large statistics and can reveal crucial information on the higher
order correlation strength. It is observed that the correlator shows a smooth
decreasing trend towards higher multiplicity classes. The increase in parton
activity that causes the parton’s terminal strings to recombine$-$ known as
color reconnection, produces a weaker trend in the correlation function. This
trend is weakest in case of a color reconnection model CR=2 which is based on
the movement of gluons in order to produce a string of minimum tension.
Further, spherocity investigation shows that the correlations are strongest
for events dominated by hard scattering (lower 0-20% spherocity class).

##  IV Summary

Event-by-event fluctuation of the mean transverse momentum ($\langle p_{\rm
T}\rangle$) has been studied as a function of average charged particle
multiplicity ($\langle N_{\rm ch}\rangle$) in p$-$p collisions at $\sqrt{s}$ =
7 TeV and 13 TeV using the PYTHIA 8 event generator. The charged particles
were selected from a kinematic range of $0.15<p_{\rm T}<2$ GeV/c and
$|\eta|<0.8$. The dynamical fluctuations of the mean transverse momentum is
caused by fluctuation in the early stages of the collision systems and have
been measured using two-, three-, and four-particle correlators. A
characteristic decrease of the two-particle particle correlator with the
charged particle multiplicity is observed. This phenomenon is attributed to
independent superposition of multi-partonic interaction in p$-$p collisions.
The effects of various microscopic processes like color reconnection and
multi-partonic interactions has also been investigated. The mechanism of color
reconnections with rope formation scenario showed higher correlation. Higher
order (third and fourth) correlation functions, which have been reported for
the first time, did not, however, exhibit the same behaviour. The study also
suggests an energy dependence of the mean transverse momentum fluctuation in
p$-$p collisions at $\sqrt{s}$ = 7 TeV and 13 TeV. Furthermore, the
fluctuation observable has also been studied in different intervals of
transverse spherocity to understand the relative contributions originating
from hard scattering and underlying events. It was observed that the
correlations were always stronger for events with di-jet-like topology. The
present study in p$-$p collisions would act as a baseline for the upcoming
measurements in A$-$A as well as p$-$p collisions at the LHC.

##  V Acknowledgements

The authors would like to thank the Department of Science and Technology
(DST), India for supporting the present work.

## References

  * [1]
  * [2] J. Adams et al, STAR Collaboration, Phys. Rev. C 72, 044902 (2005).
  * [3] B. B. Abelev et al, ALICE Collaboration , Eur. Phys. J. C 74, 3077 (2014).
  * [4] T. A. Trainor, Phys. Rev. C 92, 024915 (2015).
  * [5] L. Stodolsky, Phys. Rev. Lett. 75, 1044 (1995).
  * [6] J. Adams er al, STAR Collaborationl, Phys. Rev. C 99, 044918 (2019).
  * [7] S. Voloshin, Phys. Lett. B 632, 490-494 (2006).
  * [8] B. Schenke, Chun Shen and Derek Teaney, Phys. Rev. C 102, 034905 (2020).
  * [9] P. Bartalini, E. Berger, B. Blok, G. Calucci, R. Corke, et al., arXiv:1111.0469 [hep-ph]
  * [10] Torbjörn Sjöstrand, Stefan Ask, Jesper R Christiansen, Richard Corke, Nishita Desai, Philip Ilten, Stephen Mrenna, Stefan Prestel, Christine O Rasmussen, and Peter Z Skands, Comput. phys. commun. 191, 159–177 (2015).
  * [11] B. B. Abelev et al, ALICE Collaboration, Phys. Lett. B 727, 371 (2013).
  * [12] Torbjörn Sjöstrand, arXiv:1310.8073,(2013).
  * [13] Christian Bierlich and Jesper Roy Christiansen, Phys. Rev D 92,094010 (2015).
  * [14] Jesper Roy Christiansen and P. Z.Skands , JHEP, 08, 003 (2015).
  * [15] T. S. Biro, H. B. Nielson and J. Knoll, Nucl. Phys. B 245, 449, (1984).
  * [16] C.Bierlich, G. Gustafson, L. Lonnblad and A. Tarasov, J. High Energ. Phys. 2015: 148, (2015).
  * [17] Ranjit Nayak, Subhadip Pal, and Sadhana Dash, Phys. Rev. D 100, 074023 (2019).
  * [18] Pritam Chakraborty and Sadhana Dash, Phys. Rev. C 102, 055202 (2020).
  * [19] Ankita Goswami, Ranjit Nayak, Basanta Kumar Nandi, and Sadhana Dash, Eur. Phys. J. C 81, 988 (2021).
  * [20] T. Pierog, Iu. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, Phys. Rev C 92, 034906 (2015).
  * [21] A. Banfi, G. P. Salam, and G. Zanderighi, JHEP 06, 038 (2010).
  * [22] A. Ortiz, G. Paic, and E. Cuautle, Nucl. Phys. A 941 , 78-86 (2015).
  * [23] S. Acharya et al, ALICE Collaboration , Eur. Phys. J. C 79, 857 (2019).
  * [24] Giuliano Giacalone, Fernando G. Gardim, Jacquelyn Noronha-Hostler, and Jean-Yves Ollitrault, Phys. Rev. C 103, 024910 (2021).


# Topological data analysis of truncated contagion maps

Florian Klimm<EMAIL_ADDRESS>Department of Computational Molecular
Biology, Max Planck Institute for Molecular Genetics, Ihnestraße 63-73,
D-14195, Berlin, Germany  Department of Computer Science, Freie Universität
Berlin, Arnimallee 3, D-14195 Berlin, Germany

###### Abstract

The investigation of dynamical processes on networks has been one focus for
the study of contagion processes. It has been demonstrated that contagions can
be used to obtain information about the embedding of nodes in a Euclidean
space. Specifically, one can use the activation times of threshold contagions
to construct contagion maps as a manifold-learning approach. One drawback of
contagion maps is their high computational cost. Here, we demonstrate that a
truncation of the threshold contagions may considerably speed up the
construction of contagion maps. Finally, we show that contagion maps may be
used to find an insightful low-dimensional embedding for single-cell RNA-
sequencing data in the form of cell-similarity networks and so reveal
biological manifolds. Overall, our work makes the use of contagion maps as
manifold-learning approaches on empirical network data more viable.

††preprint: AIP/123-QED

> It is known that the analysis of spreading processes on networks may reveal
> their hidden geometric structures. These techniques, called contagion maps,
> are computationally expensive, which raises the question of whether they can
> be methodologically improved. Here, we demonstrate that a truncation (i.e.,
> early stoppage) of the spreading processes leads to a substantial speed-up
> in the computation of contagion maps. For synthetic networks, we find that a
> carefully chosen truncation may also improve the recovery of hidden
> geometric structures. We quantify this improvement by comparing the
> topological properties of the original network with the constructed
> contagion maps by computing their persistent homology. Lastly, we explore
> the embedding of single-cell transcriptomics data and show that contagion
> maps can help to distinguish different cell types.

##  I Introduction

The study of spreading processes has a rich tradition in epidemiology Kermack
and McKendrick (1991); Keeling and Rohani (2011); Schlosser  _et al._ (2020),
the social sciences Valente (1996); Christakis and Fowler (2013); Guilbeault,
Becker, and Centola (2018); Granovetter (1978); Wiedermann  _et al._ (2020);
Oh and Porter (2018), and applied mathematics Borge-Holthoefer  _et al._
(2013); Porter and Gleeson (2014); Gleeson (2008). Due to improved
availability of interaction data, the description of contagion processes on
networks has been a focus in recent years Watts and Strogatz (1998). It is
known that spreading phenomena are influenced by the geometry of a networks’s
underlying embedding Brockmann and Helbing (2013), for example for the spread
of pandemics Marvel  _et al._ (2013) and the adaption of technologies Ammerman
and Cavalli-Sforza (1971). By constructing _diffusion maps_ , the interplay
between dynamical processes and underlying geometry can be used to detect
manifold structure in high-dimensional data Rohrdanz  _et al._ (2011);
Ferguson  _et al._ (2011); Coifman  _et al._ (2008); Kim, Eykholt, and Salas
(1999). More recently, _contagion maps_ have been proposed as a technique to
detect manifold structure in networks Taylor  _et al._ (2015). Contagion maps
may reliable detect manifold structure in noisy data, even if standard
approaches, such as isomap fail Mahler (2020). Generalisation of these
approaches, for example to contagions on simplicial complexes Kilic and Taylor
(2022) and Kleinberg-like networks Mahler (2021), is an active field of
research. The construction of contagion maps, however, is computationally
expensive, which makes the exploration of strategies for their speed-up a
pressing issue. In the context of random walks on networks, it has been
established that stopping a dynamical process and restarting it may improve
the performance of machine-learning algorithms Tong, Faloutsos, and Pan
(2006); Kim, Lee, and Lee (2008); Shun  _et al._ (2016), raising the questions
whether similar approaches might be able to enhance contagion maps.

Figure 1: Truncated contagion maps can improve recovery of topological
features of a networks’ underlying manifold while having reduced computational
cost. (a) We show a schematic representation of the analysis pipeline in this
manuscript: With contagion maps we can compute a node embedding for a network.
With persistent homology, we compute barcode diagrams to quantify the recovery
of topological features in the learned embeddings. This approach is used for
the full contagion map and the newly introduced truncated contagion map. (b)
We compare the topological recovery of full contagion maps with that of
truncated contagion maps. For appropriate parameter choices, truncated
contagion maps show higher performance than full contagion maps while taking
considerably less time to compute (numbers for a noisy ring lattice with
$N=1000$ nodes).

In this work, we introduce a generalisation of contagion maps that we call
_truncated contagion maps_ (see Fig. 1). The core idea is that an early
stoppage of the contagion processes that are used to detect the manifold may
drastically reduce the computational cost of the embedding algorithm. In
particular, our numerical observations indicate a quadratic scaling in problem
size for full contagion maps and a sub-quadratic scaling for truncated
contagion maps. Furthermore, as we will show, the activations in the early
stages of a contagion often follow the underlying manifold of the networks
more clearly than the later stages. Omitting these late-stage dynamics results
in an improved manifold detection, which we verify with techniques from
_topological data analysis_.

Topological data analysis is a field that aims to develop tools that allow the
extraction of qualitative features in data that are hidden to standard
approaches. A wide range of tools and methods for the topological analysis of
data have been developed Otter  _et al._ (2017); Bianconi (2015); Bick  _et
al._ (2021) and applied to problems in biologyStolz, Harrington, and Porter
(2017); Giusti, Ghrist, and Bassett (2016), social sciences Stolz, Harrington,
and Porter (2016); Feng and Porter (2021, 2020), and many others Stolz  _et
al._ (2020); Salnikov  _et al._ (2018); Sizemore  _et al._ (2018). _Persistent
homology_ Edelsbrunner, Harer  _et al._ (2008); Edelsbrunner (2013) in
particular is a topological data analysis method that received a lot of
interest from practitioners because it enables the extraction of topological
features in high-dimensional point-cloud data. In particular, topological data
analysis is considering higher-order interactions between data points by
constructing simplicial complexes.

Advances in _single-cell transcriptomics_ enable the measurement and analysis
of gene expression at a single-cell resolution Aldridge and Teichmann (2020).
One challenge in the analysis of single-cell transcriptomics data is that
experimental advances lead to an explosion in data set sizes Lähnemann  _et
al._ (2020). As complex organisms have thousands to ten thousands of genes,
their gene-expression space is high-dimensional, making single-cell
transcriptomics data notoriously hard to visualise. Yet, the vast majority of
this space is empty and cells follow few low-dimensional manifolds that
represent, for example, cell types or differentiation trajectories. The
presence of low-dimensional structures embedded in high-dimensional gene-
expression spaces makes single-cell transcriptomics a promising field for the
application of manifold-learning approaches Moon  _et al._ (2018). In this
work, we apply truncated contagion maps to single-cell transciptomics data and
show that the constructed embeddings are fruitful in representing the complex
biological structure, which represents cell types and developmental
trajectories.

The remainder of this article is organised as follows. In Section II, we
summarise the existing methods and introduce truncated contagion maps. In
Section III, we investigate the temporal development of the contagion
processes. In Section IV, we demonstrate that truncated contagion maps may
improve the recovery of topological features in contagion maps. In Section V,
we discuss the improved computational complexity of truncated contagion maps
versus full contagion maps. In Section VI, we apply truncated contagion maps
to single-cell transcriptomics data of developing mouse oocytes. In Section
VII, we discuss our findings.

##  II Methods

###  II.1 Noisy geometric networks

Noisy geometric networksTaylor  _et al._ (2015) are geometric networks
Barthélemy (2011) to which ‘noisy’ non-geometric edges have been added. They
are defined as follows. Consider a manifold $\mathcal{M}$ that is embedded in
an ambient space $\mathcal{A}$ such that $\mathcal{M}\subset\mathcal{A}$.
First, we construct nodes $V$ that are embedded in the manifold $\mathcal{M}$
such that $\mathbf{w}^{(i)}\in\mathcal{M}$ for nodes $i\in V$. Second, we
construct edges which are of one of two types, either _geometric_ or _non-
geometric_. The geometric edges $E^{(G)}$ follow the manifold such that
$i,j\in E^{(G)}$ if the distance along the manifold between nodes $i$ and $j$
is below some distance threshold. The non-geometric edges $E^{(NG)}$ are
constructed in a random process that does not consider the manifold. We obtain
the noisy geometric network as $G=(V,E^{(G)}\cup E^{(NG)})$. Noisy geometric
networks are a type of small-world network Watts and Strogatz (1998) which
consist of long-range and short-range connections. The interplay of these two
types of interactions has been discussed in the context of many types of
networks, for example, social networks Granovetter (1973) and traffic networks
Brockmann and Helbing (2013).

In this manuscript, we focus on noisy ring lattices $NRL(N,d^{(G)},d^{(NG)})$
as prototypical examples of noisy geometric networks. Their parameters are the
number $N$ of nodes, the geometric degree $d^{(G)}$, and the non-geometric
degree $d^{(NG)}$. The NRL networks are defined as follows. First, we place
$N$ nodes uniformly spaced along the unit circle
$\mathcal{M}=\\{(a,b)|a^{2}+b^{2}=1\\}\subset\mathbb{R}^{2}$ such that
$\mathbf{w}^{(i)}=(\cos(2\pi i/N),\sin(2\pi i/N))$ is the location of node
$i$. Second, we construct geometric edges such that each node is connected to
its $d^{(G)}$ nearest-neighbour nodes. Third, we construct $N\times d^{(NG)}$
non-geometric edges uniformly at random such that each node is has exactly
$d^{(NG)}$ non-geometric edges. In practice, the non-geometric edges are
constructed with _stub-matching_ : We create a list of stubs that contains
each node exactly $d^{(NG)}$ times. Then we construct edges iteratively by
drawing without replacement uniformly at random from this list while rejecting
self-edges or parallel edges. In rare cases this procedure might fail to
terminate (e.g., if only stubs of the same node are remaining) but this
problem can be overcome by restarting the procedure. The ratio
$\alpha=d^{(NG)}/d^{(G)}$ of non-geometric to geometric edges quantifies the
‘noisiness’ of a noisy ring lattice with $\alpha=0$ indicating that there are
exclusively geometric edges, and $\alpha=1$ indicating that there are as many
non-geometric edges as geometric edges. In Fig. 2, we show two example noisy
geometric lattices with $\alpha=0$ and $\alpha=1$, respectively.

Figure 2: Noisy geometric networks contain geometric edges that follow the
underlying manifold and non-geometric edges that ignore the underlying
manifold. Here, we show noisy geometric lattices of size $N=50$, whose
underlying manifold $\mathcal{M}$ is the unit circle
$\\{(a,b)|a^{2}+b^{2}=1\\}$ with ambient space
$\mathcal{A}\subset\mathbb{R}^{2}$. The geometric lattice contains exclusively
geometric edges and has $d^{(G)}=2$. After adding $d^{(NG)}=2$ non-geometric
edges per node, we obtain a noisy geometric lattice with $\alpha=1$.

###  II.2 Watts threshold model

The Watts threshold model (WTM) is a well-established deterministic binary-
state discrete-time model of social contagions on networks Watts (2011) and is
a modification of Mark Granovetter’s threshold model Granovetter (1978). The
WTM is defined as follows. Each node is in one of two states, either active or
inactive. At time $t=0$ all nodes but a set of _seed nodes_ $S$ is inactive.
All nodes have the same activation threshold $T\in[0,1]$ (homogenous
threshold), which determines how easily a node adapts their neighbours’
activation. An active node stays active but an inactive node is activated if
more than a threshold fraction $T$ of its $d$ neighbours are active. We update
all node states synchronously at each time point and stop the contagion once a
steady state is reached (for a discussion of synchronous vs asynchronous
update procedures, see Porter and Gleeson (2016)).

###  II.3 Contagion maps

As introduced by Taylor  _et al._ , we use the WTM to construct deterministic
embeddings for networks. These embeddings are based on the nodes’ activation
times under different initial conditions. Specifically, we define a contagion
map as $V\to\\{\mathbf{x}^{(i)}\\}_{i\in V}$, where
$\mathbf{x}^{(i)}=[x_{1}^{(i)},x_{2}^{(i)},\dots,x_{j}^{(i)}]$ with
$x_{j}^{(i)}$ indicating the activation time of node $i$ for contagion $j$.
For the initial conditions, we use _cluster seeding_ such that the
neighbourhood $\mathcal{N}(j)$ of node $j$ is active at $T=0$. After a steady
state of the WTM is reached some nodes might still be inactive. The infinite
activation time of these inactive nodes is set to $2N\ll\infty$ because $2N$
is much larger than the largest possible activation time of $N-1$ (An
activation time of $N-1$ would be reached by a single node if in each time
step exactly one additional node is activated, for example, along a line
graph.)

The obtained matrix $\mathbf{X}$ of activation times is not necessarily
symmetric. To study an embedding that defines a semimetric, we investigate the
symmetric contagion map
$\mathbf{X}_{\text{symmetric}}=\mathbf{X}+\mathbf{X}^{\top}$. See Algorithms 1
and 2 for pseudocode for the contagion maps and contagion dynamics,
respectively.

For clarity, we will refer to the contagion maps as defined by Taylor  _et
al._ as ’full contagion maps’ in contrast to the ‘truncated contagion maps’
that we will introduce now.

input : Network $G=(V,E)$ with $N=|V|$ nodes

Threshold $T$

Number of steps $s$ // Ignored for full contagion maps

output : Node embedding $\mathbf{X}_{\text{symmetric}}\in\mathbb{N}^{N\times
N}$

begin

for  _each node $j\in V=\\{1,\dots,N\\}$_ do

Initialise cluster seeding by $S$ infecting $j$ and its neighbours;

$\mathbf{X[:,j]}\leftarrow$ Compute (truncated) contagion dynamics;

end for

$\mathbf{X}_{\text{symmetric}}=\mathbf{X}+\mathbf{X}^{\top}$;

return $\mathbf{X}_{\text{symmetric}}$

end

Algorithm 1 Computing (truncated) contagion maps

input : Network $G=(V,E)$ with $N=|V|$ nodes

seed nodes $S\subset V$

threshold $T$

max number of steps $s$

output : Activation time $\vec{x}\in\mathbb{N}^{N}$

begin

Initialise activation times $\vec{x}=\vec{\text{NaN}}$

Initialise step number $t=0$

while  _$\vec{x}\neq\vec{x}$ _ do

for  _each node $i\in V=\\{1,\dots,N\\}$_ do

Check if node $i$ activated

Update $\vec{x}$

end for

// Break condition only checked for truncated contagion maps

if  _$t==s$_ then

Set activation time to $2\times s$ for inactive nodes

break while loop

end if

end while

return $\vec{x}$

end

Algorithm 2 Compute (truncated) contagion dynamics

###  II.4 Truncated contagion maps

In this work, we introduce _truncated contagion maps_ as a generalisation of
full contagion contagion maps. Intuitively, they represent full contagion maps
in which contagions are not run until a steady state is reached. Rather, we
stop the contagion after $s\in\mathbb{N}_{>0}$ steps. Specifically, we define
a truncated contagion map as $V\to\\{\mathbf{x}^{(i)}\\}_{i\in V}$, where
$\mathbf{x}^{(i)}=[x_{1}^{(i)},x_{2}^{(i)},\dots,x_{j}^{(i)}]$ with

$\displaystyle x_{j}^{(i)}=\begin{cases}a_{j}^{(i)}&\text{for}\
a_{j}^{(i)}\leq s\\\ 2s&\text{else}\,.\end{cases}$ (1)

where $a_{j}^{(i)}$ indicates the activation time of node $i$ for contagion
$j$. For the initial conditions, we use _cluster seeding_ such that the
neighbourhood $\mathcal{N}(j)$ of node $j$ is active at $T=0$. We set the
activation time of unactivated nodes to $2s$, such that full contagion maps
are a special case of truncated contagion maps for $s=N$, but other choices
for the handling of infinite activation times are possible. In practice,
however, we anticipate an approximately similar behaviour of truncated and
full contagion maps for $s\ll N$ because most contagions reach a steady state
after a small number $s$ of steps (in comparison to the network size $N$). As
for full contagion maps, we study exclusively the symmetric contagion maps
$\mathbf{X}_{\text{symmetric}}=\mathbf{X}+\mathbf{X}^{\top}$ such that we
obtain a semimetric.

###  II.5 Topological data analysis

Persistent homology is a widely adapted method from topological data analysis.
It aims to identify stable (i.e., persistent) topological features of point
cloud data across a range of resolutions. To obtain these topological
features, we must construct simplicial complexes from the point cloud data.
Specifically, we employ an Vietoris–Rips simplicial filtration, which is a
sequence of simplicial complexes that constructs topological structure based
on the data. Given a data set $X$, a metric $d:Y\times Z\to\mathbb{R}$, and a
scale parameter $\alpha\in\mathbb{R}$, we can construct a Vietoris–Rips
complex as

$\displaystyle\\{{VR}(X,d,\alpha):=\\{\sigma\subset X:d(p,q)\leq\alpha\
\forall\ p,q\in\sigma\\}\,.$ (2)

Computing the Vietoris–Rips complex for a selection of scale parameters
$\alpha_{1}\leq\alpha_{2}\leq\dots\leq\alpha_{a}$, then yields a Vietoris–Rips
simplicial filtration

$\displaystyle\emptyset\subseteq\text{VR}(X,d,\alpha_{1})\subseteq\text{VR}(X,d,\alpha_{2})\subseteq\dots\subseteq\text{VR}(X,d,\alpha_{a})\,,$
(3)

which is a sequence of embedded simplicial complexes. This filtration can be
though of as a dynamical process on a simplicial complex in which faces are
added stepwise. The most common way to visualise the persistence of
topological features across a filtration is a _barcode_ , which indicates in
which filtration step a topological feature is born and in which step it
ceases to exist. We use the Ripser implementation Tralie, Saul, and Bar-On
(2018) for the computation of Vietoris–Rips persistence barcodes and use the
Euclidean norm. To summarise the 1D persistence barcodes (i.e., the presence
of loops in the data), we compute _ring stability_ as $\Delta=l_{1}-l_{2}$,
where $l_{1}$ and $l_{2}$ are the persistence of the longest and second
longest 1ifetime of loops in the data. This procedure has previously been used
to identify whether contagions are predominantly spreading by wavefront
propagation versus new cluster appearance Taylor  _et al._ (2015). In Fig. 1,
we show two example Vietoris–Rips persistence barcodes constructed from
contagion maps.

##  III Temporal development of contagions reveals suitable truncation values

Figure 3: Bifurcation diagram of WTM on noisy ring lattices matches well the
observed ring stability of truncated contagion maps. We plot the critical
thresholds of WFP (solid line, Eq. 4) and ANC (dashed line, Eq. 5) versus the
noisiness $\alpha$. These curves divide the $(T,\alpha)$ parameter space into
four qualitatively different contagion regimes: exclusively WFP, exclusively
ANC, WFP and ANC, and neither. We compare the theoretical bifurcation diagram
with the empirically observed ring stability $\Delta$ in truncated contagion
maps with $s=20$ of noisy ring lattices of size $N=400$ with geometric degree
$d^{(G}=40$ and varying non-geometric degrees $d^{(NG}\in[0,40]$ such that
their noisiness $\alpha\in[0,1]$. We find that only in the parameter regime of
exclusive WFP the ring stability $\Delta$ is high, indicating a good recovery
of the underlying manifold’s topology. We highlight the four parameter choices
that we show in Figs. 4 and 5 with crosses.

Figure 4: The temporal development of the contagion is strongly dependent on
the threshold $T$. (a) The size $q(t)$ of the contagion is growing over time
for thresholds $T\in\\{0.05,0.2,0.3\\}$ and reaches $q(t)=1$, which indicates
a global contagion (i.e., all nodes are active). The size $q(t)$ of the
contagion remains constant for $T=0.45$, indicating that no cascade is
triggered. (b) We show the number $E(t)$ of edges involved in the activation
process and distinguish between geometric edges (solid lines) and non-
geometric edges (dashed lines). The inlay shows a zoom onto the contagion with
threshold $T=0.3$ in which we observe that for the early (i.e., $t<20$)
contagions steps almost no non-geometric edges are involved.

In this section, we will describe the spread of WTM contagions on noisy ring
lattices. In particular, we will demonstrate that in certain parameter
regimes, the initial spread is mainly occurring along the manifold (called
‘wave front propagation’, WFP) and at later stages of the contagion the spread
is less restricted to the manifold and leads to new contagion clusters at
distant parts of the manifold (‘appearance of new clusters’, ANC).

Taylor  _et al._ Taylor  _et al._ (2015) proved two critical thresholds for
noisy ring lattices of noisiness $\alpha$:

$\displaystyle T^{\text{(WFP)}}$ $\displaystyle=1/(2+2\alpha)\,,\ \text{and}$
(4) $\displaystyle T^{\text{(ANC)}}$ $\displaystyle=\alpha/(1+\alpha)\,,$ (5)

which indicate the thresholds above which no WFP and no ANC, respectively,
occur. We show these critical thresholds, which intersect at
$(\alpha,T)=(1/2,1/3)$, in Fig. 3 and observe four qualitatively different
contagion regimes: exclusively WFP, exclusively ANC, WFP and ANC, and neither.
As we will show, these regimes generalise from full contagion maps to
truncated contagion maps. In the following, we investigate four parameter
choices in more detail. Specifically, we summarise the Watts threshold model
contagions for thresholds $T\in\\{0.05,0.2,0.3,0.45\\}$ on a noisy ring
lattice with $N=400$ nodes and noisiness of $\alpha=d^{(NG)}/d^{(G)}=2/6=1/3$
(see Fig. 4).

In Fig. 4a, we show the size $q(t)$ of the contagion over time (i.e., the
number of nodes that are active). First, we note that for
$T\in\\{0.05,0.2,0.3\\}$ a global cascade is triggered as $q(t)\to 1$, whereas
for $T=0.45$ the contagion size remains constant. Second, a smaller threshold
$T$ results in a quicker contagion spread which in turn yields a smaller time
until a global cascade is reached. By comparing the chosen thresholds with the
critical thresholds of WFP (Eq. 4) and ANC (Eq. 5) we observe that $T=0.05$
and $T=0.2$ fall in the regime in which WFP and ANC occur, $T=0.3$ falls in
the regime in which exclusively WFP occurs, and $T=0.45$ results in no
contagion (see coloured crossed in Fig. 3).

To study whether the contagion spreads along the manifold or independent of
it, we investigate along which type of edges the contagion spreads
predominantly. In Fig. 4b, we show the number $E(t)$ of activating edges
(i.e., edges incident to a node that is active at $t$ and a node that becomes
active at $t+1$). We distinguish between geometric edges (solid lines) and
non-geometric edges (dashed lines). We observe that for thresholds
$T\in\\{0.05,0.2\\}$ the contagion spreads quickly and to similar extent along
geometric and non-geometric edges. For a threshold of $T=0.3$ in contrast, we
find a slow contagion process. We also find that initially the contagion
process is strongly dominated by contagion along geometric edges (see inlay).
After time $t\approx 20$, however, the influence of non-geometric edges
increases. Our observation indicates that for slow contagion processes, the
contagion initially spreads along the manifold as WFP. Yet, a spreading along
non-geometric edges becomes more important over time. This highlights that the
derived critical thresholds (Eqs. 4 and 5) are strictly valid only for the
early contagion steps and become less appropriate the longer the contagion
progresses. These observations indicate that a truncation choice that removes
these late-stage contagion steps might improve the contagion map. In this
case, for example, a truncation parameter $s\in[20,60]$ might be appropriate.
As a simple, heuristic choice, we suggest $s^{\ast}=N/10=40$, which we will
explore in more detail in Sec. IV. As we will show in Sec. IV, such
truncations allow us to improve the recovery of topological features with
truncated contagion maps in comparison with full contagion maps.

##  IV Topological data analysis of truncated contagion maps

Figure 5: Truncated contagion maps applied to noisy ring lattices can recover
the ring manifold for appropriate thresholds $T$. We applied contagion maps to
noisy ring lattices with $N=400$ nodes degrees $=(d^{(G)},d^{(NG)})=(6,2)$. In
each column, we show the results for one of the four thresholds
$T\in\\{0.05,0.2,0.3,0.45\\}$. (top row) The network with nodes embedded on
the ring manifold. Each node’s colour indicates the activation time from small
(blue) to large (yellow) in one realisation of the WTM. Green nodes are seed
nodes and grey nodes are inactive at the steady state. (other rows) We show
two-dimensional projections of the $N$-dimensional contagion maps. A node’s
colour indicate its activation time as in the top row. The rows represent
truncated contagion maps with $s=10$, $s=40$, and $s=60$, and a full contagion
map, respectively. We observe that the threshold $T=0.3$ the contagion maps
most closely resemble the circular embedding manifold. In this regime,
intermediate step numbers $s=40$ leads to a best recovery of the circular
manifold (up to rotation).

Figure 6: Truncated contagion maps enable an improved recovery of topological
features, in comparison with a full contagion map. We quantify the recovery of
the topological features by computing the ring stability $\Delta$ of contagion
maps on noisy ring lattices with $N=400$. (a) For all contagion maps, the ring
stability is largest in the threshold regime
$\frac{\alpha}{1+\alpha}<T<\frac{1}{2+2\alpha}$ in which WFP dominates the
spreading process. The truncated contagion maps with $s=10$ (red line) and
$s=20$ (blue line) have a larger ring stability $\Delta$ than the full
truncation map (dashed black line). (b) The ring stability $\Delta$ varies
with the number $s$ of included contagion steps. For small number $s$ of
contagion steps, the truncated contagion map results in an improved ring
stability $\Delta$ in comparison for the full contagion map (dashed black
line). For $s>60$, when ANC occurs, the ring stability $\Delta$ diminishes. As
expected the performance of the full contagion map matches the truncated
contagion map for large $s$. Results for the threshold $T=0.3$ (c) We
investigate the truncation parameter $s^{\ast}$ that yields the best recovery
of the ring manifold for noisy ring lattice networks of sizes $N\in[50,550]$.
We find that the optimal truncation increases approximately linearly with
network size $N$, which suggests a simple heuristic choice of $s^{\ast}=N/10$
(dashed line).

In this section, we analyse the topological properties of contagion maps on
noisy ring lattices. This allows us to identify parameter regimes in which the
contagion processes follow the underlying manifold. In particular, we compare
the behaviour of truncated contagion maps with that of full contagion maps.

In Fig. 6, we study contagion maps for noisy ring lattices with $N=400$ nodes
and noisiness of $\alpha=d^{(NG)}/d^{(G)}=2/6=1/3$. Each column represents a
distinct threshold $T\in\\{0.05,0.2,0.3,0.45\\}$, which correspond to the
example thresholds in Fig. 4 and distinct regimes in the bifucation analysis
of the critical thresholds critical thresholds (Eqs. 4 and 5). In the top row,
we show the noisy ring lattice with nodes placed at their locations
$\mathbf{w}^{(i)}=(\cos(2\pi i/N),\sin(2\pi i/N)$. The colour of each node
indicates the node’s activation time during one realisation of the WTM in
which the green nodes are the seed nodes. For all other nodes, the colour map
ranges from blue (small activation time) to yellow (large activation time).
Grey nodes never adopt the contagion. As expected, we observe that for two
smallest thresholds $T\in\\{0.05,0.2\\}$ the contagion quickly spreads to the
whole network and that for the largest threshold $T=0.45$ no contagion beyond
the seed nodes occurs. For intermediary threshold $T=0.3$, however, we observe
that the contagion spreads approximately along the circular manifold.

In the remaining rows of Fig. 6, we show the point clouds $\mathbf{x}^{(i)}$
obtained from contagion maps. Specifically, rows two, three, and four show
truncated contagion maps for $s=10$ steps, $s=40$ steps, and $s=60$ steps and
the last row shows the full contagion map, which is equivalent to $s=N$. To
visualise the $N$-dimensional point clouds, we use a two-dimensional
projection through principal component analysis (PCA). Similar to the networks
in the top row, we colour each point with the activation time under one
realisation of the WTM. Yet, as we truncate the contagion after $s$ steps a
node might be grey (i.e., never activated) either because a steady state is
reached or because the contagion stopped earlier. Therefore, the number of
grey (i.e., never activated) nodes varies with the number $s$ of considered
contagion steps although the seed nodes were identical.

In general, we find that the contagion maps $\mathbf{x}^{(i)}$ for threshold
$T=0.3$, which is in a regime in which there is exclusively WFP predicted,
most closely resemble the two-dimensional ring topology of the original
network (up to rotation). Furthermore, in this regime, the truncation
parameter $s$ has a strong influence on the obtained embedding. For a small
number $s=10$ of contagion steps, the ring structure is perturbed and appears
folded. Increasing the number of contagion steps to $s=20$ makes the ring
structure more clearly visible. Increasing the number $s$ of contagion steps
further, however, results in a ‘broadening’ of the circular point clouds. Our
observation suggest that for the threshold $T=0.3$ a truncated contagion map
might better recover the underlying circular manifold than the full contagion
map.

After the bifurcation analysis in Sec. III and the qualitative investigation
in Fig. 5, we now quantify the recovery of the manifolds topological features
with persistent homology. In particular, we know that the underlying manifold
$\mathcal{M}$ of a noisy ring lattice is the unit circle. By comparing how
closely the persistent homology of the contagion maps resemble the topological
properties of the unit circle, we can quantify the recovery of its topological
properties. A unit circle has one circular hole and thus its first Betti
number is $b_{1}=1$. Therefore, we compute the ring stability $\Delta$ as a
quantification of how ‘ring-like’ a contagion map is. In Fig. 6a, we show the
ring stability $\Delta$ of contagion maps on noisy ring lattices of size
$N=400$ and noisiness $\alpha=1/3$. We vary the threshold $T\in[0,0.6]$ and
highlight the critical thresholds $T^{\text{(ANC)}}=0.25$ and
$T^{\text{(WFP)}}=0.375$ (Eqs. 4 and 5) as vertical dotted lines. All
contagion maps have their largest ring stability $\Delta$ in the regime
between these thresholds (i.e., when only WFP is possible). Furthermore, the
truncated contagion maps with $s=10$ and $s=20$ have a higher ring stability
$\Delta$ than the full contagion map, which matches our observations in Fig. 5
that an appropriate truncation of the contagion map may improve the recovery
of the ring-like topology. In other threshold regimes, the ring stability
$\Delta$ is almost indistinguishable between truncated and full contagion
maps.

In Fig. 6b, we investigate the influence of the number $s$ of contagion steps
on the ring stability $\Delta$ of contagion maps with a threshold of $T=0.3$.
We find that for a large number $s\to N$ of steps, the ring stability is the
same as for the full contagion map (dashed line), which matches our
theoretical expectation (see ‘Truncated contagion maps’ in the Methods
section). For a small number $s<60$ of steps, the truncated contagion map has
a larger ring stability $\Delta$ than the full contagion map, whereas for
$s>60$ steps, the full contagion map outperforms the truncated contagion map.
We explain this behaviour as a trade-off between data quality and data
quantity: The early-stage contagions more closely follow the underlying
manifold, yielding better embeddings for the truncated contagion maps. Yet,
full contagions have more data as less entries in the contagion map have
infinite activation times.

The dependence of the ring stability $\Delta$ on the truncation parameter $s$
raises the question whether the optimal choice varies with the network size
$N$. In Sec. III, we suggested $s^{\ast}=N/10$ as a simple, heuristic choice
because in this regime, the contagion is predominantly spreading along then
manifold. In Fig. 6c, we compare this heuristic with the actual optimal
$s^{\ast}$ for noisy ring-lattice networks of sizes $N\in[50,550]$. We find
that the optimal truncation parameter $s^{\ast}$ is increasing with the
network size $N$ and is in a decent agreement with the suggested heuristic.

Thus far, we fixed the noisiness $\alpha$ of the networks. In Fig. 3, we
compare the analytically predicted bifurcation diagram with the ring stability
$\Delta$ of truncated contagion maps with $s=20$ steps for various thresholds
$T\in[0,0.6]$ for noisy ring lattices with $N=400$ and various noisiness
values $\alpha\in[0,1]$. Parameter combinations $(\alpha,T)$ with large ring
stability $\Delta$ closely match the regime in which the bifurcation analysis
predicts exclusively WFP. This indicates that in regimes which are dominated
by WFP, the truncated contagion map consistently recovers the underlying
topology of the unit circle, whereas in other parameter regimes this is not
possible.

##  V Computational complexity of (truncated) contagion maps

Figure 7: Truncated contagion maps have a smaller computational complexity
than full contagion maps. (a) We show the computation time $\delta t$ for
contagion maps of noisy ring lattices of various sizes $N\in\\{32,10^{4}\\}$,
noisiness $\alpha=d^{(NG)}/d^{(G)}=6/2$, threshold $T=0.3$. The truncated
contagion map (red disks; $s=20$ steps considered) is for small network sizes
$N<100$ of similar computational cost as the full contagion map (blue squares)
but is considerably faster for large networks. We estimate the computational
complexity with least-squares fits of a linear equation to the logarithm of
the computation time. The computation time $\delta t$ for full and truncated
contagion maps scale super-quadratic ($\gamma\approx 2.24$) and sub-quadratic
($\gamma\approx 1.79$) with the network size $N$, respectively. (b) We show
the computation time $\delta t$ for contagion maps of noisy ring lattices of
size $N=10,000$, noisiness $\alpha=d^{(NG)}/d^{(G)}=6/2$, and threshold
$T=0.3$. We vary the number $s$ of steps and identify an approximately
$\gamma\approx 0.58$, with truncated contagion maps approaching the
computational time of full contagion maps for $s\to N$.

Taylor  _et al._ Taylor  _et al._ (2015) established that the typical
computational cost of full contagion maps is $\mathcal{O}(N^{2}d)$, where $N$
is the number of nodes in the input network and $d$ the average degree of its
nodes . As truncated contagion maps omit some activation steps, we anticipate
that less computational time is needed than for the full contagion map.

We empirically verify the reduction of computational cost for contagion maps
with thresholds $T=0.3$ on noisy ring lattices with noisiness $\alpha=1/3$ and
consider $s=20$ steps in the truncated contagion maps, which is a parameter
combination which yields a good recovery of the ring manifold (see Sec. IV).
We vary the network sizes $N\in\\{32,\dots,10^{4}\\}$ to identify the
computational complexity (see Fig. 7a). We find that for small network sizes
$N\leq 100$, the computation time $\delta t$ of truncated and full contagion
maps are indistinguishable because a steady state is reached before a
truncation after $s=20$ influences the contagions. Yet, for networks of sizes
$N>100$, the truncated contagion map is quicker to compute than the full
contagion map. Furthermore, we identify that the gap between computational
costs widens with increasing network size $N$. To obtain estimates of the
computational complexity of the form $\delta t=\zeta N^{\gamma}$, we fit its
linearised form $\log(\delta t(\log(N)))=\log{\zeta}+{\gamma}\log(N)$ to the
logarithm of the data. We find that the computational time of the full
contagion map and the truncated contagion map scale superquadratic ($\delta
t\propto N^{2.24}$) and subquadratic ($\delta t\propto N^{1.79}$) with the
problem size $N$, respectively, yielding an improved scaling of almost
$\sqrt{N}$ for the truncated contagion maps. The reduced computational
complexity of the truncated contagion maps can be understood by investigating
the time until a steady state is reached. If comparing networks of the same
class (i.e., similar characteristics, such as characteristic path length and
mean degree), in larger networks the time until a global cascade has been
reached tends to be longer than in smaller networks. Thus, a restriction to a
fixed number $s$ of steps by the truncated contagion maps, omits a larger
number of steps for larger networks. This increasing omission reduces the
computational time for larger networks, leading to a substantial reduction of
the computation time for truncated contagion maps.

We further investigate the influence of the number $s$ of steps on the
computational time $\delta t$ in Fig. 7b. Empirically, we find that the
computational time grows sublinear with $s$ such that $\delta t\propto
s^{0.58}$. As expected, for large number $s$ of steps, the computational time
of the truncated contagion maps approaches the computational time of the full
contagion maps, highlighting the need for truncation to reduce the
computational complexity.

##  VI Contagion maps of single-cell transcriptomics data

Figure 8: Truncated contagion map uncovers differentiation trajectory in
single-cell transcriptomics data of mouse oocytes. (a) We show a truncated
contagion map $\mathbf{x}^{(i)}$ with threshold $T=0.3$ and truncation $s=10$
for a cell-similarity network obtained from single-cell transcriptomics data
Gu  _et al._ (2019). The UPGMA method with Euclidean distance was used for the
hierarchical clustering. (b) A two-dimensional embedding of the truncated
contagion map reveals the oocyte developmental trajectory from primordial
oocytes over primary, secondary, and tertiary oocytes to antral oocytes
(colour indicates cell identity). (c) The cell-group separation $C$ is highest
for a truncation parameter $s=15$ and outperforms a full contagion map
(dashed, horizontal line). As expected, for large number $s$ of steps, the
performance of the truncated contagion map approaches the performance of the
full contagion map.

In this section, we investigate a single-cell transcriptomics dataset of
growing mouse oocytes Gu  _et al._ (2019); Zhang  _et al._ (2021), in which
the oocyte undergo several key stages to generate mature oocytes.
Specifically, Gu  _et al._ measure gene expression in five successive stages:
primordial follicle, primary oocytes, secondary oocytes, tertiary oocytes, and
antral oocytes. The raw gene expression matrix
$\mathbf{Count}\in\mathbb{N}^{(c\times g)}$ in which entry
$\mathbf{Count}_{i,j}$ indicates the strength of expression of gene $i$ in
cell $j$ is available from GEO under GSE114822. For further analysis, we use
scanpyWolf, Angerer, and Theis (2018) and follow current best practices
Luecken and Theis (2019). For example, we filter out 5937 genes that are
detected in less than 3 cells (for details, see online material). We construct
a cell-similarity graph as a $k$-nearest-neighbour graph $G=(V,E)$ in which
the $N=223$ nodes represent the cells $V\in\\{1,\dots,c\\}$ and each node is
connected to its $k=20$ nearest neighbours in gene-expression space.

Then, we construct a truncated contagion map with threshold $T=0.3$ and
truncation $s=15$ from the cell-similarity graph (see Fig. 8). The contagion
map $\mathbf{x}^{(i)}$ in Fig. 8a has a block-diagonal structure in which the
the non-growing oocytes (i.e., primordial) are clearly separated from the
growing oocyte (primary, secondary, and tertiary) and the fully-grown oocyte
(antral). A two-dimensional PCA projection of the contagion map
$\mathbf{x}^{(i)}$ in Fig. 8b revels the oocyte development trajectory in
gene-expression space. The primary oocytes are developing into secondary and
tertiary oocytes and finally in antral oocytes. The separation of the
developmental stages is not sharp but rather a continuum, yet a pseudo-
temporal ordering that matches the cell identities is observable.

To investigate the performance of the truncated contagion maps in separating
the different cell types, we compute the _cell-group separation_ $C$, as the
ration of median distance between cells of different cell type to the median
distance between cells of the same cell type. We investigate this ration $C$
in dependence of the number $s$ of steps in Fig. 8c. We find that for most
truncation parameters $s$, the cell-group separation $C$ is improved for the
truncated contagion maps versus the full contagion map (dashed line). The
maximum separation (i.e., best performance) is reached for $s=15$. We note
that this optimal choice is close to the heuristic choice of
$s^{\ast}=N/10\approx 22$.

Our study of the single-cell transcriptomics data shows that contagion maps
can reveal a manifold structure in a high-dimensional gene expression space by
constructing a cell-similarity graph. In particular, we observe that
primordial oocytes are clearly separated from developing oocytes and within
the developing oocytes, we observe an ordering in gene-expression space.

##  VII Conclusion

In this work, we introduced truncated contagion maps, an extension of
contagion maps that has improved computational complexity and is thus
computationally less expensive. This is achieved by truncating (i.e.,
stopping) the contagion processes before they reach a steady state, a
technique that is well-established in the context of random walks Tong,
Faloutsos, and Pan (2006). Furthermore, we demonstrate by quantifying the
persistent homology that, for appropriate parameter choices, truncated
contagion maps might improve the recovery of the embedding manifold of noisy
geometric networks. We explain this behaviour with the dominance of spreading
along the manifold in early steps of a contagion, whereas at later steps the
spread along non-geometric edges increases. Lastly, we demonstrated that
contagion maps are a manifold-learning technique that can reveal low-
dimensional structure in a cell similarity network constructed from single-
cell transcriptomics data.

Our findings suggest various avenues for further studies. For example, it is
an open question whether these truncation approaches can also be extended to
the simplical complex variant of the contagion maps Kilic and Taylor (2022),
which is particularly relevant as higher-order interactions are important in
cellular interactions Klamt, Haus, and Theis (2009); Klimm, Deane, and Reinert
(2021). Furthermore, we hypothesise that a subsampling approach with selected
seeding, as used in the k-means++ clustering algorithm Vassilvitskii and
Arthur (2006), or landmark selection, as used in L-Isomap De Silva and
Tenenbaum (2002), might further reduce the computational complexity of
contagion maps. Furthermore, the establishment of more sophisticated
heuristics for the choice of the truncation parameter $s$ will be fruitful, as
its optimal choice is likely to depend on the structure of the analysed
network in a complex, non-linear way, which is not reflected in the simple
heuristic we proposed in this paper.

###### Acknowledgements.

FK is supported as an Add-on Fellow for Interdisciplinary Life Science by the
Joachim Herz Stiftung and by the Max Planck Society. We thank Martin Vingron
for fruitful discussions.

## Data Availability Statement

Data openly available in a public repository that does not issue DOIs. Python
code to reproduce the results in this paper are available under
https://github.com/floklimm/contagionMap.

## References

  * Kermack and McKendrick (1991) W. O. Kermack and A. G. McKendrick, “Contributions to the mathematical theory of epidemics–i. 1927.” Bulletin of Mathematical Biology 53, 33–55 (1991).
  * Keeling and Rohani (2011) M. J. Keeling and P. Rohani,  _Modeling infectious diseases in humans and animals_ (Princeton University Press, 2011).
  * Schlosser  _et al._ (2020) F. Schlosser, B. F. Maier, O. Jack, D. Hinrichs, A. Zachariae, and D. Brockmann, “Covid-19 lockdown induces disease-mitigating structural changes in mobility networks,” Proceedings of the National Academy of Sciences of the United States of America 117, 32883–32890 (2020).
  * Valente (1996) T. W. Valente, “Social network thresholds in the diffusion of innovations,” Social Networks 18, 69–89 (1996).
  * Christakis and Fowler (2013) N. A. Christakis and J. H. Fowler, “Social contagion theory: examining dynamic social networks and human behavior,” Statistics in Medicine 32, 556–577 (2013).
  * Guilbeault, Becker, and Centola (2018) D. Guilbeault, J. Becker, and D. Centola, “Complex contagions: A decade in review,” Complex spreading phenomena in social systems , 3–25 (2018).
  * Granovetter (1978) M. Granovetter, “Threshold models of collective behavior,” American Journal of Sociology 83, 1420–1443 (1978).
  * Wiedermann  _et al._ (2020) M. Wiedermann, E. K. Smith, J. Heitzig, and J. F. Donges, “A network-based microfoundation of Granovetter’s threshold model for social tipping,” Scientific Reports 10, 1–10 (2020).
  * Oh and Porter (2018) S.-W. Oh and M. A. Porter, “Complex contagions with timers,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 033101 (2018).
  * Borge-Holthoefer  _et al._ (2013) J. Borge-Holthoefer, R. A. Banos, S. González-Bailón, and Y. Moreno, “Cascading behaviour in complex socio-technical networks,” Journal of Complex Networks 1, 3–24 (2013).
  * Porter and Gleeson (2014) M. A. Porter and J. P. Gleeson, “Dynamical systems on networks: A tutorial,” arXiv preprint arXiv:1403.7663 (2014).
  * Gleeson (2008) J. P. Gleeson, “Cascades on correlated and modular random networks,” Physical Review E 77, 046117 (2008).
  * Watts and Strogatz (1998) D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’ networks,” Nature 393, 440–442 (1998).
  * Brockmann and Helbing (2013) D. Brockmann and D. Helbing, “The hidden geometry of complex, network-driven contagion phenomena,” Science 342, 1337–1342 (2013).
  * Marvel  _et al._ (2013) S. A. Marvel, T. Martin, C. R. Doering, D. Lusseau, and M. E. Newman, “The small-world effect is a modern phenomenon,” arXiv preprint arXiv:1310.2636 (2013).
  * Ammerman and Cavalli-Sforza (1971) A. J. Ammerman and L. L. Cavalli-Sforza, “Measuring the rate of spread of early farming in europe,” Man , 674–688 (1971).
  * Rohrdanz  _et al._ (2011) M. A. Rohrdanz, W. Zheng, M. Maggioni, and C. Clementi, “Determination of reaction coordinates via locally scaled diffusion map,” The Journal of chemical physics 134, 03B624 (2011).
  * Ferguson  _et al._ (2011) A. L. Ferguson, A. Z. Panagiotopoulos, I. G. Kevrekidis, and P. G. Debenedetti, “Nonlinear dimensionality reduction in molecular simulation: The diffusion map approach,” Chemical Physics Letters 509, 1–11 (2011).
  * Coifman  _et al._ (2008) R. R. Coifman, I. G. Kevrekidis, S. Lafon, M. Maggioni, and B. Nadler, “Diffusion maps, reduction coordinates, and low dimensional representation of stochastic systems,” Multiscale Modeling & Simulation 7, 842–864 (2008).
  * Kim, Eykholt, and Salas (1999) H. Kim, R. Eykholt, and J. Salas, “Nonlinear dynamics, delay times, and embedding windows,” Physica D: Nonlinear Phenomena 127, 48–60 (1999).
  * Taylor  _et al._ (2015) D. Taylor, F. Klimm, H. A. Harrington, M. Kramár, K. Mischaikow, M. A. Porter, and P. J. Mucha, “Topological data analysis of contagion maps for examining spreading processes on networks,” Nature Communications 6, 1–11 (2015).
  * Mahler (2020) B. I. Mahler, “Contagion dynamics for manifold learning,” arXiv preprint arXiv:2012.00091 (2020).
  * Kilic and Taylor (2022) B. U. Kilic and D. Taylor, “Simplicial cascades are orchestrated by the multidimensional geometry of neuronal complexes,” arXiv preprint arXiv:2201.02071 (2022).
  * Mahler (2021) B. I. Mahler, “Analysis of contagion maps on a class of networks that are spatially embedded in a torus,” SIAM Journal on Applied Mathematics 81, 1416–1440 (2021).
  * Tong, Faloutsos, and Pan (2006) H. Tong, C. Faloutsos, and J.-Y. Pan, “Fast random walk with restart and its applications,” in  _Sixth international conference on data mining (ICDM’06)_ (IEEE, 2006) pp. 613–622.
  * Kim, Lee, and Lee (2008) T. H. Kim, K. M. Lee, and S. U. Lee, “Generative image segmentation using random walks with restart,” in  _European conference on computer vision_ (Springer, 2008) pp. 264–275.
  * Shun  _et al._ (2016) J. Shun, F. Roosta-Khorasani, K. Fountoulakis, and M. W. Mahoney, “Parallel local graph clustering,” arXiv preprint arXiv:1604.07515 (2016).
  * Otter  _et al._ (2017) N. Otter, M. A. Porter, U. Tillmann, P. Grindrod, and H. A. Harrington, “A roadmap for the computation of persistent homology,” EPJ Data Science 6, 1–38 (2017).
  * Bianconi (2015) G. Bianconi, “Interdisciplinary and physics challenges of network theory,” EPL (Europhysics Letters) 111, 56001 (2015).
  * Bick  _et al._ (2021) C. Bick, E. Gross, H. A. Harrington, and M. T. Schaub, “What are higher-order networks?” arXiv preprint arXiv:2104.11329 (2021).
  * Stolz, Harrington, and Porter (2017) B. J. Stolz, H. A. Harrington, and M. A. Porter, “Persistent homology of time-dependent functional networks constructed from coupled time series,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 047410 (2017).
  * Giusti, Ghrist, and Bassett (2016) C. Giusti, R. Ghrist, and D. S. Bassett, “Two’s company, three (or more) is a simplex,” Journal of Computational Neuroscience 41, 1–14 (2016).
  * Stolz, Harrington, and Porter (2016) B. J. Stolz, H. Harrington, and M. A. Porter, “The topological ’shape’of Brexit,” Available at SSRN 2843662 (2016).
  * Feng and Porter (2021) M. Feng and M. A. Porter, “Persistent homology of geospatial data: A case study with voting,” SIAM Review 63, 67–99 (2021).
  * Feng and Porter (2020) M. Feng and M. A. Porter, “Spatial applications of topological data analysis: Cities, snowflakes, random structures, and spiders spinning under the influence,” Physical Review Research 2, 033426 (2020).
  * Stolz  _et al._ (2020) B. J. Stolz, J. Tanner, H. A. Harrington, and V. Nanda, “Geometric anomaly detection in data,” Proceedings of the National Academy of Sciences of the United States of America 117, 19664–19669 (2020).
  * Salnikov  _et al._ (2018) V. Salnikov, D. Cassese, R. Lambiotte, and N. S. Jones, “Co-occurrence simplicial complexes in mathematics: identifying the holes of knowledge,” Applied Network Science 3, 1–23 (2018).
  * Sizemore  _et al._ (2018) A. E. Sizemore, E. A. Karuza, C. Giusti, and D. S. Bassett, “Knowledge gaps in the early growth of semantic feature networks,” Nature Human Behaviour 2, 682–692 (2018).
  * Edelsbrunner, Harer  _et al._ (2008) H. Edelsbrunner, J. Harer, _et al._ , “Persistent homology-a survey,” Contemporary Mathematics 453, 257–282 (2008).
  * Edelsbrunner (2013) H. Edelsbrunner, “Persistent homology: theory and practice,” (2013).
  * Aldridge and Teichmann (2020) S. Aldridge and S. A. Teichmann, “Single cell transcriptomics comes of age,” Nature Communications 11, 1–4 (2020).
  * Lähnemann  _et al._ (2020) D. Lähnemann, J. Köster, E. Szczurek, D. J. McCarthy, S. C. Hicks, M. D. Robinson, C. A. Vallejos, K. R. Campbell, N. Beerenwinkel, A. Mahfouz, _et al._ , “Eleven grand challenges in single-cell data science,” Genome Biology 21, 1–35 (2020).
  * Moon  _et al._ (2018) K. R. Moon, J. S. Stanley III, D. Burkhardt, D. van Dijk, G. Wolf, and S. Krishnaswamy, “Manifold learning-based methods for analyzing single-cell RNA-sequencing data,” Current Opinion in Systems Biology 7, 36–46 (2018).
  * Barthélemy (2011) M. Barthélemy, “Spatial networks,” Physics Reports 499, 1–101 (2011).
  * Granovetter (1973) M. S. Granovetter, “The strength of weak ties,” American Journal of Sociology 78, 1360–1380 (1973).
  * Watts (2011) D. J. Watts, “A simple model of global cascades on random networks,” in  _The Structure and Dynamics of Networks_ (Princeton University Press, 2011) pp. 497–502.
  * Porter and Gleeson (2016) M. A. Porter and J. P. Gleeson, “Dynamical systems on networks,” Frontiers in Applied Dynamical Systems: Reviews and Tutorials 4 (2016).
  * Tralie, Saul, and Bar-On (2018) C. Tralie, N. Saul, and R. Bar-On, “Ripser. py: A lean persistent homology library for python,” Journal of Open Source Software 3, 925 (2018).
  * Gu  _et al._ (2019) C. Gu, S. Liu, Q. Wu, L. Zhang, and F. Guo, “Integrative single-cell analysis of transcriptome, dna methylome and chromatin accessibility in mouse oocytes,” Cell Research 29, 110–123 (2019).
  * Zhang  _et al._ (2021) H. Zhang, M. Esposito, M. G. Pezet, J. Aryaman, W. Wei, F. Klimm, C. Calabrese, S. P. Burr, C. H. Macabelli, C. Viscomi, _et al._ , “Mitochondrial DNA heteroplasmy is modulated during oocyte development propagating mutation transmission,” Science Advances 7, eabi5657 (2021).
  * Wolf, Angerer, and Theis (2018) F. A. Wolf, P. Angerer, and F. J. Theis, “Scanpy: large-scale single-cell gene expression data analysis,” Genome biology 19, 1–5 (2018).
  * Luecken and Theis (2019) M. D. Luecken and F. J. Theis, “Current best practices in single-cell RNA-seq analysis: a tutorial,” Molecular systems biology 15, e8746 (2019).
  * Klamt, Haus, and Theis (2009) S. Klamt, U.-U. Haus, and F. Theis, “Hypergraphs and cellular networks,” PLoS Computational Biology 5, e1000385 (2009).
  * Klimm, Deane, and Reinert (2021) F. Klimm, C. M. Deane, and G. Reinert, “Hypergraphs for predicting essential genes using multiprotein complex data,” Journal of Complex Networks 9, cnaa028 (2021).
  * Vassilvitskii and Arthur (2006) S. Vassilvitskii and D. Arthur, “k-means++: The advantages of careful seeding,” in  _Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms_ (2006) pp. 1027–1035.
  * De Silva and Tenenbaum (2002) V. De Silva and J. B. Tenenbaum, “Global versus local methods in nonlinear dimensionality reduction.” in  _NIPS_ , Vol. 15 (2002) pp. 705–712.