# Linking Through Time:
Memory-Enhanced Community Discovery in Temporal Networks

Giulio Virginio Clemente<EMAIL_ADDRESS>IMT School for Advanced
Studies, Piazza San Francesco 19, 55100 Lucca (Italy)  Diego Garlaschelli  IMT
School for Advanced Studies, Piazza S. Francesco 19, 55100 Lucca (Italy)
Lorentz Institute for Theoretical Physics, Niels Bohrweg 2, 2333 CA Leiden
(The Netherlands)  INdAM-GNAMPA Istituto Nazionale di Alta Matematica (Italy)

###### Abstract

Temporal Networks, and more specifically, Markovian Temporal Networks, present
a unique challenge regarding the community discovery task. The inherent
dynamism of these systems requires an intricate understanding of memory
effects and structural heterogeneity, which are often key drivers of network
evolution. In this study, we address these aspects by introducing an
innovative approach to community detection, centered around a novel modularity
function. We focus on demonstrating the improvements our new approach brings
to a fundamental aspect of community detection: the detectability threshold
problem. We show that by associating memory directly with nodes’ memberships
and considering it in the expression of the modularity, the detectability
threshold can be lowered with respect to cases where memory is not considered,
thereby enhancing the quality of the communities discovered. To validate our
approach, we carry out extensive numerical simulations, assessing the
effectiveness of our method in a controlled setting. Furthermore, we apply our
method to real-world data to underscore its practicality and robustness. This
application not only demonstrates the method’s effectiveness but also reveals
its capacity to indirectly tackle additional challenges, such as determining
the optimal time window for aggregating data in dynamic graphs. This
illustrates the method’s versatility in addressing complex aspects of temporal
network analysis.

††preprint: APS/123-QED

##  I Introduction

Networks often possess intricate mesostructures that conceal crucial
information about the system’s functionality. Among these mesostructures, the
presence of communities or groups has garnered significant interest within the
scientific community. This interest is primarily due to their profound
connection with various characterizing aspects of the system, such as the
diffusion of processes within it [1], its stability [2, 3], and the ability to
predict link formation [4].

When considering the network’s topological structure, communities emerge as
clusters of nodes exhibiting a significantly higher link density than
expected. Based on this intuition, numerous methodologies have been developed
to tackle community detection, ranging from the use of objective functions to
quantify the quality of a partition, such as the modularity function [5], to
approaches based on generative models like the Stochastic Block Model [6], and
flow models that leverage network dynamics to cluster nodes based on
persistent flow patterns [7].

Fundamental is the determination of the ”expected”, a concept at the heart of
the primary tool used in this work: the modularity function. This function
quantifies the deviation between a given benchmark, which defines what is
’expected,’ and the observed number of links within the network. This
benchmark, often referred to as a null model, plays a crucial role in
obtaining accurate results, a role that we are going to remark in the rest of
the paper.

While many of the above-mentioned methodologies have been well-developed for
static networks, they encounter challenges when it comes to generalizing them
for time-varying graphs. This difficulty arises primarily due to the
additional degree of freedom introduced by time, which needs to be considered.
Nevertheless, community detection for temporal networks is a rapidly expanding
field, and several studies have endeavored to adapt static methods to the
dynamic context [8, 9, 10].

When talking about modularity function, it is the absence of appropriate null
models that constitutes a significant barrier to generalizing it to a temporal
framework. In this paper, we fill this gap by demonstrating how the maximum
entropy approach [11] can be utilized to construct various forms of the
modularity function and leverage them to account for different aspects of
network structure. Specifically, in the subsequent sections, we introduce a
novel formulation of the modularity function based on a recent Markovian null
model for temporal networks [12]. This formulation allows simultaneous control
over anomalies related to link persistence between nodes and link formation,
thus allowing an accurate generalization of the modularity function to time-
varying graphs. Additionally, we illustrate how incorporating memory between
blocks can generally lead to a dramatic improvements in algorithmic
performance.

To demonstrate the advancements achieved by leveraging these properties, we
investigate the problem of the Detectability Threshold, which serves as a
theoretical limit defining the boundary for the ability of any algorithm to
detect the presence of communities [13]. Initially, we introduce the threshold
for the static case, and subsequently, by incorporating the insights presented
in [14], we compute the threshold while considering the presence of memory
between nodes within the same blocks.

We supported our theoretical results with numerical simulations by conducting
a series of analysis on synthetic temporal networks, wherein both the number
of nodes and memberships remain fixed throughout the network’s evolution. Our
findings reveal that when the membership of a group of nodes influences both
the probability of link formation and the probability of link persistence,
algorithm performance consistently improves when considering also information
regarding the latter property.

Furthermore, we also conducted empirical analysis by utilizing a real-world
dataset of proximity relations in a primary school [15]. We illustrate how our
proposed methodology, can help with the selection of an optimal time window
for data aggregation, and can be effectively used for community detection when
memory is driven by the nodes’ membership.

Our synthetic and empirical analyses aim to underline the importance of
considering temporal characteristics and memory in network structures, and to
highlight how these considerations can enhance community detection methodology
in temporal networks.

##  II Community detection, Null models
and Modularity Function

Consider a graph $G=(N,E)$, where $N$ represents a set of nodes and
$E\subseteq N\times N$ denotes the set of edges. A community structure within
$G$ is defined as a partition $g=\\{g_{1},g_{2},...,g_{q}\\}$ of $N$, where
each $g_{i}\subset N$ is a cluster of nodes such that the density of edges
within each $C_{i}$ is significantly greater than the density of edges between
$g_{i}$ and $N\setminus g_{i}$. In this context, we focus on non-overlapping
communities, implying that each node belongs to exactly one $g_{i}$ within
$g$.

Discovering the optimal subdivision for a network with an inherent block
structure is a challenging task that has captured the attention of numerous
researchers [16]. To address this issue, Newman and Girvan introduced a
quality function in 2004, known as modularity [5], to evaluate the
effectiveness of network partitions. The original formulation of this function
is expressed as follows:

$Q(g)=\sum_{i,j\in N}\left[a_{ij}-\langle
a_{ij}\rangle\right]\delta_{g(i)g(j)}$ (1)

Where $a_{ij}$ denotes the observed element of the adjacency matrix $A$ for
the graph $G$, indicating the presence of a link between node $i$ and node
$j$. On the other hand, $\langle a_{ij}\rangle$ represents the expected value
of the same link in a null model that assumes the absence of any community
structure. The vector $g$ indicates the membership of each node within the
community structure.

The modularity function measures the deviation between the observed and
expected number of intra-community edges as defined by the null model.
Therefore, to identify the optimal partition, our objective is to maximize the
modularity function with respect to the assigned node membership, effectively
seeking the community structure with the most significant deviation with
respect to the specified null model.

When the optimal partition is such that the two quantities are similar, the
modularity function approaches zero, indicating a lack of strong community
structure. Conversely, when a robust community structure exists, the
modularity value tends to approach 1. As already said, it is important to note
that the cornerstone of this measure lies in the definition of the null model
(a pattern-generating model) that assumes specific properties remain constant
while allowing other characteristics to vary stochastically [17]. The
randomization process aims to recreate certain aspects of the data that would
be expected in the absence of any particular mechanisms, such as the presence
of a community structure.

It is precisely by leveraging this characteristic that the modularity function
accomplishes its objective.

The choice of the null model is subjective and depends on the properties we
aim to reproduce in our random network. However, much attention has been given
to null models that accurately replicate the degree sequence of undirected and
unweighted networks, primarily to study the impact of local inhomogeneity
[18]. In this paper, since our main contribution relies on null models that
precisely reproduce the degree sequence of a given network, we will
exclusively consider null models that achieve this objective.

When constraining the degree sequence, for binary and undirected networks, the
widely used specification for $\langle A_{ij}\rangle$ is the well-known Chung-
Lu model: $\langle a_{ij}\rangle=\frac{k_{i}k_{j}}{2m}$, where $k_{i}$ and
$k_{j}$ represent the degrees of nodes $i$ and $j$, respectively, and $2m$ is
the total number of possible links. However, it is important to note that the
Chung-Lu model is an approximation that holds under stringent conditions [19],
which many real networks do not satisfy.

An important example of such a deviation is highlighted in the study conducted
by Maslov et al. [20], where they demonstrate that applying the Chung-Lu model
to compute the probability of connecting two hubs in a snapshot of the
internet network on January 2, 2000 would result in a probability of 43.5.
Clearly, this value cannot be interpreted as a probability. Utilizing such a
quantity in the modularity function would introduce bias, ultimately impacting
the performance of the method.

In this work, we propose a modularity function derived from null models
developed through a Maximum Entropy approach [11]. This approach overcomes the
limitations associated with the Chung-Lu model and avoids the aforementioned
bias.

###  II.1 Maximum Entropy Null Models

The maximum entropy approach is a methodology that enables the definition of a
probability distribution over a set of elements, such as graphs, in a manner
that ensures specific expected values are reproduced while all other aspects
remain maximally random. This means that if certain characteristics emerge
from the model without being explicitly imposed, it is indicative of a genuine
correlation between these properties and the constraints utilized in
constructing the model.

Applying the maximum entropy formalism to real-world problems involves a two-
step process. The first step addresses a graph G with certain distinctive
properties C(G). In this step, we seek to obtain the most unbiased probability
distribution P(G) over the set of all possible graphs, ensuring that the
expected value of C(G) aligns with the imposed constraints. To achieve this
distribution, we leverage the Shannon-Gibbs entropy:

$S(P(G))=-\sum_{G}P(G)ln(P(G))$ (2)

Imposing that the probability distribution reproduces on average

$\langle C\rangle=\sum_{G}C(G)P(G)=C^{*}$ (3)

Where G is an ensemble element, $C^{*}$ is the measured value of the
constraint, and $\langle\cdot\rangle$ indicates the expected value over the
ensemble.

It is crucial to emphasize that, as we have defined it, this class of models
allows for the imposition of constraints where the expected value, rather than
the exact value for each network realization, is reproduced. This has
implications not only on computational aspects [11] but also on the model’s
ability to replicate metrics not predetermined by constraints [21].

The constrained maximization problem is approached by employing the method of
Lagrange multipliers ($\theta$). The resulting probability distribution has
the following functional form [22]:

$P(G|\vec{\theta})=\frac{e^{-H(G|\vec{\theta})}}{Z(\vec{\theta})}$ (4)

Where $H(G|\vec{\theta})=\sum_{i}\theta_{i}C_{i}$ is called graph Hamiltonian,
and $Z(\vec{\theta})=\sum_{G}e^{-H(G|\vec{\theta})}$ is the normalization
constant. It is important to note that the analytical expression for $Z$ is
often unknown, making it challenging to obtain an exact expression for $P(G)$.
However, in our case, we will consider only situations where solutions have
been found.
The second step involves determining the numeric values of the parameters
$\vec{\theta}$ that reproduce the measured quantities on a real graph. To
tackle this problem, we can employ the maximum likelihood principle, which
entails setting the parameters to the specific values that maximize the
likelihood function $P(G^{*}|\vec{\theta})$ given the model [23].

The first maximum entropy model we introduce is the Binary Configuration Model
[22], which serves as the counterpart to the Chung-Lu model constructed using
the maximum entropy formalism. As previously mentioned, this represents the
primary model employed in the classical definition of modularity. The Binary
Configuration Model is obtained by imposing constraints on the degree
sequence, leading to the following form of the Hamiltonian:

$H(\mathbf{G}|\vec{\theta})\equiv\sum_{i=1}^{N}\theta_{i}k_{i}$ (5)

Where $k_{i}$ is the degree of node i. Hence for this specific problem, we
have that

$\langle a_{ij}\rangle\equiv\frac{x_{i}x_{j}}{1+x_{i}x_{j}}.$ (6)

where $x_{i}\equiv e^{-\theta_{i}}$ for each node $i$ [23, 11].

In the subsequent sections, we will utilize the maximum entropy formalism to
employ specific null models for defining new modularity expressions and
characterizing their properties.

###  II.2 Detectability Threshold In Static Graphs

In order to explore the impact of temporal information on community detection,
we investigate the behavior of different modularity functions in relation to
the detectability threshold problem. This phenomenon has garnered significant
attention from researchers in recent years [24, 25]. It has been demonstrated
that there exists a fundamental limitation to the ability of algorithms to
detect community structures in networks, characterized by a phase transition
known as the detectability threshold (DT) [25].

To introduce this issue, we consider networks generated using the well-known
Stochastic Block Model (SBM) [26]. The SBM is a generative model that produces
networks with a block structure, making it a suitable framework for studying
the performances of community detection algorithms in a controlled
environment. By analyzing networks generated by the SBM, we can gain insights
into the behavior of community detection algorithms and investigate how
temporal information can potentially enhance their performance.

Given a set of nodes N and q groups denoted by
$g=\\{g_{1},g_{2},...,g_{q}\\}$, the SBM in its basic formulation defines a
network structure in which unweighted edges are placed between node i and node
j with probability $p_{g_{i}g_{j}}$ depending on the groups which they belong.
Hence the probability of observing a graph G is:

$P(G|g)=\prod_{(ij)\in
N}p_{g_{i}g_{j}}^{a_{ij}}\left(1-p_{g_{i}g_{j}}\right)^{1-a_{ij}}$ (7)

Where $a_{ij}$ is the element of the adjacency matrix A corresponding the
graph $G$.
In the following, we will refer to a particular class of block model defined
by clusters having the same dimension: Planted Partition Model, specifically
all the theoretical results we show are related to a Planted Partition model
with two communities that will be characterized by a probability of connection
between nodes of different clusters $p_{out}$ and within nodes of the same
cluster $p_{in}$.
For this specific case, it has been shown rigorously [27] that given a sparse
graph G, i.e. with $p_{ij}=O(\frac{1}{N})$, with two communities, an average
degree k, an average number of connections between groups
$k_{out}=\frac{N}{2}p_{out}$, and within groups $k_{in}=\frac{N}{2}p_{in}$, it
exists a value $k_{in}-k_{out}=\Delta^{*}$ such that, above this threshold
(DT), no algorithm can capture the community structure better than by chance:

$\Delta^{*}=2\sqrt{k}$ (8)

This expression of the threshold represents a result that will be used
subsequently as a reference to demonstrate the benefit of using multiple
realizations of the same network, rather than the single network on which the
previous expression was obtained.

In the following, to measure the quality of the partition computed by our
models, we use the the Adjusted Rand Index (ARI) [28] defined in the
following:

$ARI(g_{t},g_{e})=\frac{\sum_{ij}\binom{n_{ij}}{2}-\left[\sum_{i}\binom{a_{i}}{2}\sum_{j}\binom{b_{j}}{2}\right]/\binom{N}{2}}{\frac{1}{2}\left[\sum_{i}\binom{a_{i}}{2}+\sum_{j}\binom{b_{j}}{2}\right]-\left[\binom{a_{i}}{2}\sum_{j}\binom{b_{j}}{2}\right]/\binom{N}{2}}$
(9)

Where, index i goes over all the sets in $g_{t}$, while index j goes over all
sets in $g_{e}$.
$n_{ij}$ defines the number of elements in common between the sets in
$g_{t}(i)$ and $g_{e}(j)$, $a_{i}=\sum_{j}n_{ij}$ is the sum of the elements
in common between set i and all sets in j, and $b_{j}\sum_{i}n_{ij}$ is the
sum of the elements in common between the set j and all the possible sets j.
Hence, the value of ARI varies between 0, when none of the elements is well
classified, and 1, when there is a perfect matching between the real
membership and the one inferred. In this way we can compare the communities
obtained by the modularity maximization with the one given as ground-truth.

As a first experiment, we use the modularity function defined with the null
model in eq. (6), denoted as:

$Q(g)=\sum_{i,j}\left[a_{ij}-\frac{x_{i}x_{j}}{1+x_{i}x_{j}}\right]\delta_{g(i)g(j)}$
(10)

we can examine how the community detection algorithm behaves when varying the
relationship between $p_{in}$ and $p_{out}$. Since the expected degree is the
same for each node, a suitable null model for our case involves a single
parameter, in contrast to the binary configuration model (CM) that requires N
parameters. However, the choice of the null model does not significantly
affect the final results, so we employ the binary configuration model for
illustrative purposes. Figure (1) displays the transition between a regime
where the modularity function-based algorithm successfully detects the
community structure and a regime where it fails to do so.

Figure 1:  In figure are shown the performances using the modularity defined
in eq. (10) on a static network with k = 8 and 200 nodes. Each point is
averaged over 30 iteration.

The displayed plot results from the examination of graphs, each composed of
two communities, in which the relationship between $p_{in}$ and $p_{out}$
varies as previously described. These artificially generated graphs maintain a
constant average node degree. The plot substantiates the theoretical forecast
that a distinct transition point, denoted as $\Delta=\Delta^{*}$, exists. The
experiments confirms the existence of a point that demarcates the boundary
between regions where community structure is detectable and where it remains
undetectable.

##  III Modularity for Temporal Networks with Stationary Data Generating
Process

The transition from static to temporal graphs brings an added dimension that
can significantly impact tasks like community detection. The network’s
evolution over time, as noted in [29], can be driven by alterations in the
community structure.

To deal with these evolving systems, it is necessary to acknowledge that
temporal networks can be represented through various methods. Thus, selecting
an approach that best aligns with our proposed algorithm is crucial. For our
study, we face the community detection problem working with a snapshot
representation, namely a sequence of T static sparse graphs, each denoted as
$G_{i}$. In our context, every graph in the sequence maintains a constant
number of nodes, N, and is characterized as unweighted and undirected. We
refer to this series as a ’graph trajectory’, which is denoted as follows:

${\mathcal{G}}\equiv\\{\mathbf{G}_{1},\dots,\mathbf{G}_{T}\\}$ (11)

As said, communities in temporal networks can experience a variety of
transformations [30]. Our focus, however, is on persistent communities that
actively shape the network’s evolution. Though this represents a very simple
case of study, it is worth investigating because, as we’ll show, it helps to
demostrate how the right null model can greatly boost community detection
accuracy, and, in this specific case, to show how memory can be enlightening
in community discovery.

Figure 2: Schematic Representation of Temporal Networks with Fixed
Communities. Each node in the network is assigned a color to indicate its
membership that remains fixed throughout the evolution of the network.

Our study postulates that networks evolve with stable community memberships
and parameters, as illustrated in Fig. 2. This implies that communities are
assumed to maintain consistency across the graph trajectory, generated by a
stationary process. We particularly explore two different data generating
processes: initially, a time-varying graph with no correlation between
snapshots, subsequently evolving into a temporal network with inherent memory.

###  III.1 Temporal Networks with independent observations

The first Data Generating Process (DGP) that we analyze generates a graph
trajectory characterized by a sequence of independent and sparse snapshots,
each of them built according to a Planted Partition Model, the very same model
used for the static case. In this way the average connectivity for each block
in each snapshot is equal to k, implying that the average connectivity in
time, $\sum_{t}\frac{1}{T}k(t)=\overline{k}=k$. Hence, the model so defined
has the following expression, representing the probability of observing a
given graph trajectory $\mathcal{G}$:

$P(\mathcal{G}|g)=\prod_{t=1}^{T}\prod_{(ij)}p_{g_{i}g_{j}}^{a_{g_{i}g_{j}}(t)}\left(1-p_{g_{i}g_{j}}\right)^{1-a_{g_{i}g_{j}}(t)}$
(12)

We refer to this model as a Dynamic Stochastic Block Model (DSBM).

To address community detection within this framework, we employ a modularity
function that utilizes as a null model a naive temporal generalization of the
Binary Configuration Model. This model is constructed by substituting the
degree sequence of a single snapshot with its temporal average as a
constraint. Consequently, the related Hamiltonian is expressed as follows:

$H(\mathcal{G}|\vec{\alpha})\equiv\sum_{i=1}^{N}\alpha_{i}\sum_{t=1}^{T}\frac{k_{i}(t)}{T}$
(13)

Formally, the resulting expression for the expected value of the link
existence between two nodes is the same as the non-temporal binary
configuration model, but this time it refers to its average over time:

$\langle\overline{a_{ij}}\rangle=\frac{x_{i}x_{j}}{1+x_{i}x_{j}}$ (14)

Where, in this case, $x_{i}=e^{-\alpha_{i}/T}$.
Therefore, the related modularity can be defined as:

$Q(g)=\sum_{i,j}\left[\ \overline{a_{ij}}-\langle\overline{a_{ij}}\rangle\
\right]\delta_{g(i)g(j)}$ (15)

With $\langle\overline{a_{ij}}\rangle$ as in eq. (14) and
$\overline{a_{ij}}=\frac{1}{T}\sum_{T}a_{ij}(t)$.

Intuitively, when the community structure does not change in time, as the case
presented in this section, using eq. (15) allows to filter out the spurious
links by observing several realizations of the same network, hence directly
putting more evidence on the hidden block structure.

However, by employing this modularity function, we implicitly assume that each
snapshot represents an independent realization of the same network. This
approach overlooks the temporal order in the graph’s trajectory, thereby
discarding relevant information when it is present.

###  III.2 Community detection in Markovian
Time-Varying Graphs

We now aim to introduce a refined version of the modularity function, which
incorporates node features associated with the nodes’ capacity to maintain
links. These ideas are derived from the findings presented in [12], where the
authors demonstrate the importance of capturing node-level memory and how a
model incorporating this aspect may be crucial for analyzing real-world
systems. Based on this consideration, here ahead, we want to demonstrate how,
by leveraging these features, we can improve the performances of the
modularity function on the community detection task. To achieve the
aforementioned objectives, we consider a DGP with a Markovian structure, which
allows us to exert control over the memory embedded in the system.

The DGP employed can be seen as a Dynamic Stochastic Block Model that differs
from others of its kind found in the literature [9, 14, 31]. Indeed, contrary
to other DGPs where memory is an external parameter, independent of the
network’s characteristics such as degree distribution, in our model, we
incorporate a memory effect. This effect is embodied in a Markovian evolution
of the system, governed by node-specific features that simultaneously control
link formation and persistence, interdependently.

Therefore, to introduce the model, we begin by identifying the features to be
encoded in the maximum entropy model, which will serve as the foundation for
both the DGP and the definition of modularity. In this context, we focus on
what is termed the ’persisting degree’, the portion of the degree that is
maintained between one snapshot and the next:

$h_{i}(t,1)\equiv\sum_{j\neq i}a_{ij}(t)a_{ij}(t+1)=\sum_{j\neq
i}b_{ij}(t,1).$ (16)

We will refer to the quantity $a_{ij}(t)a_{ij}(t+1)$ as the one-lagged
persisting link.

Hence, by imposing this as further constraints at the case described by the
Hamiltonian in eq. (13), we generate the following Hamiltonian:

$H(\mathcal{G}|\vec{\alpha},\vec{\beta})\equiv\sum_{i=1}^{N}\left[\alpha_{i}\sum_{t=1}^{T}\frac{k_{i}(t)}{T}+\beta_{i}\sum_{t=1}^{T}\frac{h_{i}(t)}{T}\right]$
(17)

By solving the related problem [12], we can obtain the analytical expression
for the expected value of the average link and the average persisting-link:

$p_{ij}=\langle\overline{a_{ij}(t)}\rangle$ (18)

and

$q_{ij}=\langle\overline{a_{ij}(t)a_{ij}(t+1)}\rangle$ (19)

Both the expressions as function of $\vec{\alpha}$ and $\vec{\beta}$ (see the
supplementary information for the details).
Where $p_{ij}$ represents the marginal probability of having a link between
node i and node j at time t, and

$P[a_{ij}(t+1)=1\ \&\ a_{ij}(t)=1]=q_{ij}$ (20)

is the probability that two nodes have a connection at time t and $t+1$.

Once obtained $p_{ij}$ and $q_{ij}$, we can compute the following stochastic
matrix:

$P_{ij}(q_{ij},p_{ij})=\begin{bmatrix}\frac{q_{ij}}{p_{ij}}&\frac{p_{ij}-q_{ij}}{p_{ij}}\\\
\frac{p_{ij}-q_{ij}}{1-p_{ij}}&\frac{1-2p_{ij}+q_{ij}}{1-p_{ij}}\end{bmatrix}$
(21)

Which defines the transition matrix associated with the evolution of each
link.

This stochastic matrix essentially defines the system’s evolution, ensuring
that the expected values for the persisting degree and the degree precisely
match those constrained in eq. (17).

Operationally, our version of the Dynamic Stochastic Block Model begins with
the generation of an initial network snapshot using a static Planted Partition
Model. Subsequently, we allow the link presence to evolve according to a
Stochastic Matrix in eq. (21), creating the full trajectory.

Similar to the other cases, in our experiments we generate graph trajectories
by specifying our DGP with two distinct probabilities: $p_{out}$ and $p_{in}$
for connection, and $q_{out}$ and $q_{in}$ for preserving links. These
probabilities automatically induce an average degree, k, and an average
persisting degree, h, respectively. Consequently, the membership influences
both the link density within blocks and the persisting link density. All
probabilities are set to ensure the network’s sparsity is maintained, being of
the order $O(\frac{1}{N})$, which is a crucial condition for defining the
expression for the Detectability Threshold.

Once the graph sequence $\mathcal{G}$ is generated, we can obtain two
sequences of matrices:

$\mathcal{A}=\left(A(1),A(2)...A(T)\right)$ (22)

And

$\mathcal{B}=\left(B(1),B(2)...B(T)\right)$ (23)

Where $A(t)=\left(a_{ij}(t)\right)$ and
$B(t)=\left(a_{ij}(t+1)a_{ij}(t)\right)$.

From which we can define a new matrix

$C(t)=A(t)+B(t)$ (24)

Such that each element of $C_{ij}(t)=a_{ij}(t)+b_{ij}(t)$.

By computing the average in time of A(t) and B(t), we obtain the matrix
$\overline{\mathcal{C}}$, whose elements are used to define a new modularity
in the following way:

$Q(g)=\sum_{ij\in
N}\left[\overline{c_{ij}(t)}-\langle\overline{c_{ij}}\rangle\right]\delta(g_{i},g_{j})$
(25)

Where $\langle\overline{c_{ij}}\rangle=p_{ij}+q_{ij}$.

Through this approach, we account for two contributions that are, in
principle, interdependent yet capable of providing distinct insights into the
DGP.

##  IV Detectability Threshold for Stationary Time-Varying Graphs

This section illustrates how the DT varies based on the assumptions about the
DGP and the model applied for data analysis. As previously stated, the DGP
remains constant over time, indicating that our data are samples from a
stationary process. Under this premise, it is intuitive that an increase in
the number of samples, i.e., the length of the trajectory, reduces the
influence of noise on the observed block structures. This reduction occurs as
noise is averaged out over more samples, consequently enhancing the DT, and
leading to improved outcomes.

Nevertheless, this intuition is not the whole story. Indeed, as we are going
to show, using the right features and a suitable model to capture those
features, can positively change the value of the DT and the results of the
algorithm. Building upon findings in the literature, we demonstrate how the
implementation of the above-mentioned Maximum Entropy null model for temporal
networks alters the expression for the DT and the regions where the algorithm
works well increase.

###  IV.1 Detectability Threshold for non-correlated Time-Varying Graphs

We now recall the expression for the DT in a first and well known case, where
we consider a DGP in which every snapshot is sparse and independent, as in eq.
(12).

It is known that under the condition of sparsity and independence, an
expression for the DT has been obtained by [14] for $T\rightarrow\infty$ and
conjectured for finite T, a conjecture that has then been demonstrated in
[31]. The focus of the computation relies on the sparsity of each snapshot
that implies a locally tree-like structure for the graph trajectory
$\mathcal{G}$. Hence the authors in [14, 31], interpreted the community
detection task in temporal networks as a label recovery problem, showing that,
for our particular case in which communities remain stable in time, community
detection is feasible only if:

$\Delta_{nm}^{*}=2\sqrt{\frac{\overline{k}}{T}}$ (26)

Where we define $\Delta_{nm}=\overline{k_{in}}-\overline{k_{out}}$111nm stands
for no-memory.. For this kind of temporal network, all the information related
to the block structure is contained in the degree sequence, so the threshold
has an expression that is a direct function of this quantity. Interestingly,
as anticipated, the threshold is also a function of T and is such that as T
increases $\Delta_{nm}^{*}$ goes to zero, so, if T is big enough, we will
always spot the block structure.

Given this framework, here we can employ a modularity function that uses only
information on the degree sequence to spot the community structure.

###  IV.2 Detectability Threshold for Markovian Time-Varying Graphs

Transitioning to Markovian temporal networks introduces correlations between
snapshots, effectively embedding intrinsic memory. In our model, this memory
is assimilated by enabling link persistence, a feature previously noted for
its undesirable effects on community detection task performance [9].

Introducing snapshot correlations complicates the observed matrix sequence
$\mathcal{A}$, with links influenced more by memory than indicative of
community structure.

Furthermore, the lack of independence between consecutive snapshots leads to
the loss of prerequisites for computing the DT, as described in [14, 31].
However, it is worth noting that within our DGP, we can precisely determine
the necessary waiting period to consider two consecutive matrices
uncorrelated, thereby restoring the conditions required to calculate the
threshold.

This can be done by acknowledging that each link in our graph evolves
according to the stochastic matrix defined in eq.(21). This matrix is
characterized by two eigenvalues: the first, by definition, is equal to 1, and
the second is given by:

$\lambda_{ij}^{(2)}=\frac{q_{ij}-p_{ij}^{2}}{p_{ij}-p_{ij}^{2}}$ (27)

This eigenvalue is directly linked to the system’s memory. Indeed, by
calculating $\lambda_{ij}^{(2)}$, as detailed in [12], we can determine the
correlation length and the rate at which the correlation diminishes between
two consecutive snapshots. Consequently, we obtain $\tau^{*}$, which is
defined as:

$\tau_{ij}^{*}=\frac{1}{ln\left(\frac{1}{\lambda_{ij}^{(2)}}\right)}$ (28)

As schematically represented in Fig. (3), waiting for a duration equal to the
maximum of all $\tau_{ij}^{*}$, denoted as $\tau^{*}$, between one snapshot
and the next ensures that the snapshots can be considered independent.

Figure 3: In the figure, a schematic representation of the behavior for the
memory in time for system evolving according to eq. (21). $\tau^{*}$
represents the relaxation time.

At this point, by discounting the entries of the matrix
$\overline{\mathcal{A}}$ by the maximum correlation length $\tau^{*}$, we
reproduce the condition for the calculation of the DT for finite T, resulting
in the following expression for the DT:

$\Delta_{M,nm}^{*}=2\sqrt{\frac{\tau^{*}}{T}\overline{k}}$ (29)

Where M stands for ”Markovian”.

Moving ahead, we can state that in scenarios where the process exhibits
memory, relying only on the degree sequence to identify communities results in
the neglect of crucial information that can provide insights into the block
structure. In the following, we will demonstrate how incorporating the
persisting degree enhances the algorithm’s performance and changes the
expression for DT.

In this regard, our innovative modularity function, as defined in eq. (25),
integrates information from two distinct matrices, $\overline{\mathcal{B}}$
and $\overline{\mathcal{A}}$. This approach allows us to present how the
method proposed by Ghasemian et al. [14] to compute the DT can be re-adapted
to our context, thereby implying a new and better threshold.

We have previously demonstrated that, by adjusting $\overline{\mathcal{A}}$
with $\tau^{*}$, its elements can be interpreted as those realized by an
equivalent DSBM with independent snapshots [31, 9]. To extend this rationale
to the entire matrix $\overline{\mathcal{C}}$, it is necessary to justify how
the same discounting factor can be applied to the matrix $\mathcal{B}$. This
is achieved by explicitly calculating the autocorrelation between elements at
time $t$ and $t+1$ for matrix $\mathcal{B}$.

As detailed in the Supplementary Information, by adjusting for a time
$\tau^{*}$, even in this scenario, the elements of matrix $\mathcal{B}$ can be
considered independent. This adjustment allows for the application of the
previously discussed methodology in calculating the DT. Thus, when considering
the matrix resulting from the sum of $\mathcal{A}$ and $\mathcal{B}$, which
remains sparse, we re-establish the conditions necessary for the analytical
computation of the DT:

$\Delta_{M,m}^{*}=2\sqrt{\frac{\tau^{*}}{T}\overline{c}}$ (30)

Where $c=\overline{k}+\overline{h}$, and $\overline{h}<\overline{k}$ by
definition. This result shows that, if the block membership guides the link
persistence, this property can be used to improve the quality of the
discovered community structure.

##  V Numerical Simulations

To validate the analytical results presented in the previous sections, we have
conducted a series of numerical simulations. These simulations demonstrate the
consistency between the predicted and observed thresholds. Additionally, we
characterize the algorithm’s performance across various combinations of
parameters, showing how our novel modularity function works better in
different scenarios.

To achieve this, we utilize the Data Generating Processes previously
introduced for temporal networks, along with their associated modularity
functions. For all experiments, we optimize the modularity function using the
algorithm detailed in [32], adapting the expression for the null model
specific to each case. We generate the trajectories by varying the values of
$\Delta k$ and $\Delta h$, as defined in:

$\Delta k=k_{in}-k_{out}\quad\textit{and}\quad\Delta h=h_{in}-h_{out}$ (31)

with

$\displaystyle k_{in}=\frac{N}{2}p_{in},\quad h_{in}=\frac{N}{2}q_{in}$ (32)
$\displaystyle k_{out}=\frac{N}{2}p_{out},\quad h_{out}=\frac{N}{2}q_{out}$
(33)

Here, $N$ denotes the number of nodes within a single snapshot. In general,
each trajectory is generated to maintain both the average degree and the
average persisting degree across the entire temporal network, when those are
specified, or otherwise only the average degree.

In this setting, $\Delta k$ signifies the disparity between the average
(temporal) number of links formed within the communities and those outside,
whereas $\Delta h$ denotes the difference in the persistence of links within
these communities compared to external links.

In all the experiments, we have two communities whose elements are known and
play the role of ground-truth communities.

The initial series of experiments involves generating graph trajectories with
independent observations, as in eq. (12), by varying $k_{in}$ and $k_{out}$,
thus altering the value of $\Delta k$ while maintaining a constant average
degree $k$. By specifying the number of time steps for each trajectory, we
consequently establish the threshold $\Delta_{nm}^{*}$.

We generate 50 trajectories for each $\Delta K$ value and various $T$ values,
then plot the average performance across these 50 instances for the Community
Detection task.

Figure 4: The figure represents three different profiles for the performances
of the community detection task on graph trajectories of different lengths.
The vertical line blu represents the DT for the temporal network with T=15,
while the red vertical line is for a single static snapshot.

In Figure (4), are represented the results obtained using the modularity
function described in eq. (15). The red vertical line corresponds to the
static DT eq.(8) as computed on the first snapshot of the trajectory (see Fig.
1). It is interesting to note how even with few snapshots, the gain of
performance quality is evident. As expected, in all the cases, the DT
predicted and the ones computed numerically are the same.

In our next analysis, we explore graph trajectories that are generated through
a Markovian data-generating process, as detailed in eq. (21). Similar to our
previous experiments, we generate 50 trajectories, each for different values
of $T$ and $\Delta K$. However, in this updated scenario, we incorporate a
memory effect by keeping the values of $q_{in}$ and $q_{out}$ constant, but
allowing $p_{in}$ and $p_{out}$ to vary. This approach, as discussed in
Section (IV.2), is expected to significantly influence the performance of the
modularity, previously defined in eq. (15).

Figure 5: In figure is displayed the average ARI values for different values
of $\Delta K$ and T. The vertical lines represent the computed thresholds for
the static case (in red) and the dynamic case with T=15 (in blue).

Figure (5) illustrates the adverse effects of memory on the quality of the
detected communities. This diminishes the benefits of using additional
snapshots compared to the static case, thereby confirming the observations in
eq. (29). Notably, the curve patterns in Fig. (5) differ from those observed
in Fig. (4). Intrinsically, the value of $\tau^{*}$, which reflects the memory
level in our process, is influenced not only by $q_{in}$ and $q_{out}$ but
also by $\Delta K$, as expressed in eq. (27). This dependency impacts the
memory and so, it alters the curve shapes.

Having demonstrated how memory can impact on the ability to discover
communities in temporal networks, we now explore how it can be leveraged to
recover some of the lost information about community structure. We proceed to
the most comprehensive scenario addressed in this study, generating graph
trajectories while varying $\Delta k$ and $\Delta h$ for different values of
$T$. It’s important to note that the average value of $\overline{h}$ should
remain lower than $\overline{k}$ during trajectory generation in order to
create realistic trajectories, which, of course, it influences the value of
$\tau^{*}$. In fact, as we adjust these two parameters, while maintaining the
average values of $\overline{k}$ and $\overline{h}$ constant, the resultant
changes in $\tau^{*}$ significantly affect both the curve shapes and the value
of the DT as already mentioned and showen in Fig. (5).

To realize this task we use the Dynamic Stochastic Block model described in
Sec. (III.2).

In the experiments we choose all the parameters such that we keep the
threshold $\Delta_{m}*$ fixed, but varying the relationships between
$k_{in},k_{out},h_{in}$ and $h_{out}$, so to explore the algorithmic
performances and see whether the threshold predicted analytically is validated
numerically.

Given the process used to generate the trajectories, we ensure that the
evolution does not change the relation between $p_{in}$ and $p_{out}$,
facilitating the direct comparison between the modularity built on the model
with memory and the one built on the memory-less one.

Figure 6: *

(a)

Figure 7: *

(b)

Figure 8: *

(c)

Figure 9: *

(d)

Figure 10: *

(e)

Figure 11: *

(g)

Figure 12: In these figures, we present three heatmaps for each of the
following modularity functions: the no-memory modularity function (Figures a,
b, and c at the top) and the memory modularity function (Figures d, e, and f
at the bottom). These heatmaps illustrate the performance on three Markovian
networks of increasing lengths (T = 15, T = 50, and T = 100, respectively). In
each plot, the analytical DT is represented by a red line.

The experimental outcomes, as depicted in Fig. (12), align with expectations.
The no-memory approach fails to recognize the varying propensity of nodes
within different groups to preserve links, inadvertently discarding valuable
information.

As the number of samples generated by the same DGP increases, the accuracy of
parameter estimation improves. Consequently, a greater number of snapshots
leads to more precise inference of the community structure.

It is evident that $\Delta_{nm}^{*}>\Delta_{m}^{*}$, while there are zones
where, specifically for $\Delta h(0)$, the distinction between modularity
incorporating the persisting degree and that without it becomes negligible.
This is attributed to the neglectable disparity at this point between the
persisting links within and outside communities.

However, as illustrated in Fig. (12), for all other scenarios, incorporating
this additional information consistently enhances the quality of the detected
communities, showing how our novel modularity performs better.

##  VI Proximity Networks

We now proceed to examine a real-world application, aiming to test the
implementation of the previously defined modularities on a dataset where
ground-truth community information is accessible.

The dataset provides information about proximity relations between students
and teachers in a primary school [15]. Electronic proximity detectors worn by
the participants were used to collect the data, capturing the presence of
other nearby detectors at 20-second intervals. The data was collected over a
period of two days, specifically from October 1st to October 2nd, 2009. For
simplicity, in this analysis, we focus only on the data from the first day.
Based on this data, we aim to construct a temporal network, which will be
characterized by a series of snapshots created by aggregating data over a
specific time window. This implies that within each window, two nodes will be
connected if they have come into contact at least once during that time
interval. In our study, we leverage verified information about existing
communities. The underlying hypothesis is that the network dynamics are
primarily driven by the community memberships of the nodes. The aim is to
investigate whether, by examining these dynamics, it’s possible to discern the
node memberships. It’s important to note that these dynamics are far from
trivial, as even during class hours, students change classrooms and have the
freedom to move around.

Working with snapshot representations of real temporal networks requires to
take into account how the choice of time windows for aggregating information
can impact the quality of results, as noted by Krings et al [33] and Clauset
et al [34]. Before suggesting a solution to this challenge for our specific
scenario, we begin analyzing the real network by initially adopting a
10-minute time aggregation window. This initial aggregation will be used for
preliminary analyses, which will lay the foundation for developing a strategy
to determine the optimal time window for our task.

The proposed approach for community identification is based on the assumption
that the process that generated the data is stationary. Therefore, it becomes
crucial to identify temporal segments for which it makes sense to consider the
network as generated by the same process. To solve this task and identify
potential breakpoints in the system’s evolution, we apply the approach
introduced in [12], by using the model that takes into account memory at node
level as the reference model.

We have to emphasize that in choosing periods using this model, we are
implicitly assuming that the observed stationary period is one during which
the memory among links remains relatively stable.

We assume that the network is generated using a Dynamic Stochastic Block Model
using the ground truth communities, and we evaluate the performances of our
algorithm. We then observe the behavior of the algorithm and the quality of
the communities found in the various identified segments.

Figure 13: The figure displays the ARI values for 10 minutes temporal windows
at various identified breakpoints. Additionally, the values for different
types of null models are indicated.

In Fig. (13) are presented the performances of the algorithm when executed
over a 10-minute time window. We identify several changing points that reflect
diverse dynamic behaviors. In this particular scenario, these dynamics are
understood to result from the students’ class memberships and their movements
across various classrooms in the school. It is noteworthy that our approach
finds that influence of class-based membership diminishes significantly during
periods in the playground. This result allows us to deduce that the supposed
data generating process, the Dynamic Stochastic Block Model with which the
trajectory has formed, can be observed during class hours. However, during
playground time, the decrease in the quality of the observed communities
compared to the ground truth leads us to believe that the assumed data
generating process may be incorrect during these hours, and that the dynamics
within that interval are not well described by this model.

We now take a further step in our investigation by examining how the quality
of the communities identified changes with different time windows. Our
objective is to determine if there is a natural scale where the effects of
membership on link density and persistent link density are more evident. To
this end, we distinguish between class hours, where an higher similarity
between the discovered communities and the ground truth was noted, and
playground hours. Our algorithm is applied across various time windows. For
each time window, we identify distinct stationary segments. In each segment,
we assess the Adjusted Rand Index between the communities detected using our
modularities and the actual communities. Additionally, we characterize each
stationary segment using $\Delta k$ and $\Delta h$.

In this context, $\Delta k$ is defined as the difference between the average
(temporal) number of links formed within the ground-truth communities and
those formed outside, while $\Delta h$ represents the difference in the number
of links that persist within the ground truth communities compared to those
outside, hence the terms in eq. (31) becomes:

$\displaystyle
k_{in}=\sum_{ij}\overline{a}_{ij}^{*}\delta(g_{i}^{*},g_{j}^{*}),\quad
h_{in}=\sum_{ij}\overline{b}_{ij}^{*}\delta(g_{i}^{*},g_{j}^{*})$ (34)
$\displaystyle
k_{out}=\sum_{ij}\overline{a}_{ij}^{*}\textit{I}_{(g_{i}^{*}\neq
g_{j}^{*})},\quad
h_{out}=\sum_{ij}\overline{b}_{ij}^{*}\textit{I}_{(g_{i}^{*}\neq g_{j}^{*})}$
(35)

with, I as indicator function

$\textit{I}_{(g_{i}\neq g_{j})}=\begin{cases}1&\text{if }g_{i}\neq g_{j},\\\
0&\text{otherwise}.\end{cases}$ (36)

Figure 14: *

(a)

Figure 15: *

(b)

Figure 16: Figure displays the performance of community detection on a real
network. Each point corresponds to a stationary segment identified using
structural break detection. To each point, a value of $\Delta k$ and $\Delta
h$ is attributed. The color of the points indicates the ARI value, calculated
by comparing the communities detected with the modularity with memory against
the ground truths. In plot (a), the values during class hours are shown, while
plot (b) depicts the values during playground.The label associated with each
point indicates the size of the temporal window used to create the temporal
network associated with that segment. Only labels corresponding to ARI values
greater than 0.8 are displayed.

Figure (16) presents the outcomes of this analysis. It reveals that during
class hours, an increase in both $\Delta k$ and $\Delta h$ is associated with
a corresponding rise in the ARI. However, during playground hours, as
previously observed in Fig. (13), this correlation appears to be absent. It’s
important to emphasize that in this case, $\Delta k$ and $\Delta h$ do not
assume the same meaning as in cases involving only two partitions. Here, what
we observe is the cumulative sum computed across all possible clusters, which
amounts is greater then two in our dataset. Hence, it might be possible to
observe an increase in the difference between links due to the contribution of
one of the communities for which there was a reference and not in an
homogeneous way for all the clusters.

So our results show that for the non-playground part the performances,
regardless of the aggregation window, seem to follow the expectations for a
scenario where the evolution is genuinely driven by membership. As the
temporal window increases, there is an increase in $\Delta k$ and $\Delta h$
for certain segments, which the model can leverage to improve the quality of
the communities. During class hours, memory and community structures appear to
be related.

The situation differs during playground hours. In this latter case, it’s
probable that other effects are mixed with the existing ones, rendering the
presumed Dynamic Stochastic Block Model with the class membership less
suitable for accurately describing the system.

To strengthen the relationship between memory and the quality of the detected
communities, we propose an additional analysis. We examine if there is a
relationship between the presence of memory and the quality of the observed
communities: the hypothesis is that if the system’s evolution follows the
membership of the elements, then there may exist a natural temporal scale at
which both the quality of the communities, measured with the ARI, and the
memory are at their peak. This natural time scale would result in the
preferred time window for aggregating data in such a way that the communities
become more distinct. We measure memory with normalized autocovariance. The
normalization is done in such a way to allow us to control for the effect of
degree. In general, the average normalized autocovariance will be higher if
the degree is larger. By dividing by the degree of each node, we discount this
effect.

Hence, as a metric to this end, we employ the average autocovariance, defined
as follows:

$\displaystyle\overline{\textbf{A}}(W,\tau=1)=\frac{1}{N}\sum_{i}^{N}\frac{A_{i}(\tau=1)}{k_{i}}$
(37)

Here, the average is computed over all nodes N in the network, and

$A_{i}(\tau)=\frac{1}{T}\sum_{t=1}^{T}\big{[}h_{i}(t,\tau)-\langle
h_{i}(t,\tau)\rangle_{nm}\big{]}$ (38)

represents the autocovariance specific to node i. In this context, $\langle
h_{i}(t,\tau)\rangle_{nm}$ denotes the expected node degree persistence
according to a network randomization model that precludes memory. To
accomplish this, we utilize a null model specifically designed to not
contemplate memory effects, as defined in eq. (14). This approach retains any
genuine temporal correlations between snapshots, enabling a more accurate
investigation of memory effects.

Having identified class hours as the temporal segments where dynamics are
driven by membership, we focus on these to study the behavior of normalized
autocovariance. Our strategy involves varying the time window within a range
of 1 to 20 minutes:

$W\in\\{1^{\prime},2^{\prime},3^{\prime},\ldots,20^{\prime}\\}$ (39)

We identify all the stationary segments that fall within class hours using the
same approach as before, corresponding to the different time windows. For each
of these, we calculate the normalized autocovariance for $\tau=1$. We use the
average of the normalized autocovariance across all identified stationary
segments as a measure of memory. Thus, defining brks as the set of segments
where we can assume the memory to be constant, the memory associated with each
time window will be defined as:

$\mu_{A}(W)=\frac{1}{|\textit{brks}|}\sum_{s\in\textit{brks}}\overline{\textbf{A}_{s}}(W,\tau=1)$
(40)

Hence, once for each segment, is computed the expected value of the persisting
degree relative to a memoryless model ($\langle
h_{i}(t,\tau)\rangle_{\text{nm}}$), we can compute $\mu_{A}$.

At this stage, we need to establish a criterion for selecting the optimal time
window, which will be utilized for evaluating the efficacy of our community
detection modularity.

We define the ”best window” $W_{best}$ as the index $W$ for which $\mu_{A}(W)$
is maximum:

$W_{best}=\arg\max\mu_{A}(W)$ (41)

Figure 17: The figure represents the normalized autocovariance measurement
across three time segments for different time window.

in Fig. (17) is presented the measurement of the normalized autocovariance as
function of the time window and the average ARI computed over the segments in
class hours. Each point on the plot is designed for a non-overlapping time
window. The red points represent the maximum value taken between all the
windows. Between these two quantities we measure the Pearson Correlation,
funding a strong and significative correlation with value 0.73, demostrating a
relationship between the presence of memory and the quality of the communities
found.

In our data, $W_{best}=16\ \textit{minutes}$, which also corresponds to the
maximum on the ARI.

In support of this assessment, we have included controlled experiments
reported in the supplementary information (A.5) where the system evolves
exactly following these the hypothesis tested above.

##  VII Discussion

Effective community detection in temporal networks necessitates the
development of new, appropriate null models [35]. Such models may provide a
basis for comparison and enable the identification of dynamic structures in
networks enhancing the quality of the discovered communities. In this work, we
make a contribution in this direction by demonstrating how the use of the
right null model can lead to the discovery of hidden structures with higher
precision.

We show the potential of modularity-based methods in community detection by
utilizing Maximum Entropy null models to tailor their formulation to the
specific challenges encountered. We propose a new modularity function that
acknowledges the persistence of connections among nodes within the same
community. By focusing on the issue of the detectability threshold, our
findings reveal that selecting an appropriate null model can substantially
elevate the quality of the identified communities. This highlights the
critical role of accounting for memory effects, when present, and explains how
incorporating memory changes the detectability threshold for community
identification.

In real-world applications, we demonstrate the importance of considering
memory in community discovery, and how this feature can assist in indirectly
related challenges, such as determining the ideal time window for data
aggregation to develop a time-varying graph. We propose a data-driven
methodology to identify the most suitable time window that reveals a more
distinct delineation of community structures.

The results indicate that if memory is present but not adequately incorporated
into the model, it can deteriorate the quality of the identified communities.
However, if properly integrated into the algorithm used for community
detection, it can become a fundamental advantage.

In numerous real-world systems, the influence of community membership on
memory is a critical factor that cannot be ignored. Under these circumstances,
a function that integrates memory effects, such as the one we propose in this
study, is essential. It holds the potential to uncover new patterns that were
previously hidden. By acknowledging and incorporating memory, our approach
opens the door to revealing deeper insights into the dynamics of these
systems.

##  VIII acknowledgements

This work is supported by the European Union - NextGenerationEU - National
Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR),
project ‘SoBigData.it - Strengthening the Italian RI for Social Mining and Big
Data Analytics’ - Grant IR0000013 (n. 3264, 28/12/2021)
(https://pnrr.sobigdata.it/). D.G. acknowledges support from the Dutch
Econophysics Foundation (Stichting Econophysics, Leiden, the Netherlands) and
from the project NetRes - ‘Network analysis of economic and financial
resilience’, Italian DM n. 289, 25-03-2021 (PRO3 Scuole), CUP D67G22000130001
(https://netres.imtlucca.it).

## References

  * [1] Ide, Kiyotaka; Namatame, Akira; Ponnambalam, Loganathan; Fu, Xiuju; Siow, Rick. A Study on Relationship between Modularity and Diffusion Dynamics in Networks from Spectral Analysis Perspective, International Journal of Advanced Computer Science and Applications, 5, (2014).
  * [2] May, Robert M. Will a large complex system be stable?, Nature, 238, 5364, 413–414 (1972).
  * [3] de Castro, Francisco; Adl, Sina; Allesina, Stefano; Bardgett, Richard; Bolger, Thomas; Dalzell, Johnathan; Emmerson, Mark; Fleming, Thomas; Garlaschelli, Diego; Grilli, Jacopo; Hannula, Silja; Vries, Franciska; Lindo, Zoë; Maule, Aaron; Opik, Maarja; Rillig, Matthias; Veresoglou, Stavros; Wall, Diana; Caruso, Tancredi. Local stability properties of complex, species-rich soil food webs with functional block structure. Ecology and Evolution, 2021, 1-12 (2021).
  * [4] Kumari, Anisha; Behera, Ranjan; Sahoo, Bibudatta; Sahoo, Satya. Prediction of link evolution using community detection in social network, Computing, 104, (2022).
  * [5] Newman, Mark; Girvan, Michelle. Finding and Evaluating Community Structure in Networks, Physical review. E, Statistical, nonlinear, and soft matter physics, 69, 026113 (2004).
  * [6] Decelle, Aurelien; Krzakala, Florent; Moore, Cristopher; Zdeborova, Lenka. Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications, Physical review. E, Statistical, nonlinear, and soft matter physics, 84, 066106 (2011).
  * [7] Rosvall, Martin; Bergstrom, Carl. Maps of Random Walks on Complex Networks Reveal Community Structure, Proceedings of the National Academy of Sciences of the United States of America, 105, 1118-23 (2008).
  * [8] Mucha, Peter; Richardson, Thomas; Macon, Kevin; Porter, Mason; Onnela, Jukka-Pekka. Community Structure in Time-Dependent, Multiscale, and Multiplex Networks, Science (New York, N.Y.), 328, 876-8 (2010).
  * [9] Barucca, Paolo; Lillo, Fabrizio; Mazzarisi, Piero; Tantari, Daniele. Detectability thresholds in networks with dynamic link and community structure. arXiv preprint arXiv:1701.05804 (2017).
  * [10] Zhang, Xiao; Moore, Cristopher; Newman, Mark. Random graph models for dynamic networks, The European Physical Journal B, 90, (2016).
  * [11] Squartini, Tiziano; Garlaschelli, Diego. Maximum-Entropy Networks: Pattern Detection, Network Reconstruction and Graph Combinatorics, (2017).
  * [12] Clemente, Giulio Virginio; Tessone, Claudio J.; Garlaschelli, Diego. Temporal networks with node-specific memory: unbiased inference of transition probabilities, relaxation times and structural breaks. arXiv:2311.16981 (2023).
  * [13] Decelle, Aurelien; Krzakala, Florent; Moore, Cristopher; Zdeborova, Lenka. Inference and Phase Transitions in the Detection of Modules in Sparse Networks, Physical review letters, 107, 065701 (2011).
  * [14] Ghasemian, Amir; Zhang, Pan; Clauset, Aaron; Moore, Cristopher; Peel, Leto. Detectability Thresholds and Optimal Algorithms for Community Structure in Dynamic Networks, Physical Review X, 6, (2016).
  * [15] Stehlé, Juliette; Voirin, Nicolas; Barrat, Alain; Cattuto, Ciro; Isella, Lorenzo; Pinton, Jean-François; Quaggiotto, Marco; Van den Broeck, Wouter; Régis, Corinne; Lina, Bruno; Vanhems, Philippe. High-Resolution Measurements of Face-to-Face Contact Patterns in a Primary School, PLOS ONE, 6, e23176 (2011).
  * [16] Fortunato, Santo. Community Detection in Graphs, Physics Reports, 486, (2009).
  * [17] Gotelli, Nicholas; Graves, Gary R. Null Models in Ecology. Ecology, 78, (1996).
  * [18] Newman, M. Finding Community Structure in Networks Using the Eigenvectors of Matrices. Physical review. E, Statistical, nonlinear, and soft matter physics, 74, 036104 (2006).
  * [19] Squartini, Tiziano; Mastrandrea, Rossana; Garlaschelli, Diego. Unbiased sampling of network ensembles, New Journal of Physics, 17, (2015).
  * [20] Maslov, Sergei; Sneppen, Kim; Zaliznyak, Alexei. Detection of topological patterns in complex networks: Correlation profile of the internet, Physica A: Statistical Mechanics and its Applications, 333, 529-540 (2004).
  * [21] Caruso, Tancredi; Clemente, Giulio; Rillig, Matthias; Garlaschelli, Diego. Fluctuating ecological networks: A synthesis of maximum-entropy approaches for pattern detection and process inference. Methods in Ecology and Evolution, 13, (2022).
  * [22] Park, Juyong; Newman, M. E. J. Statistical mechanics of networks, Phys. Rev. E, 70, 066117 (2004).
  * [23] Garlaschelli, Diego; Loffredo, Maria. Maximum likelihood: Extracting unbiased information from complex networks, Physical review. E, Statistical, nonlinear, and soft matter physics, 78, 015101 (2008).
  * [24] Reichardt, Jörg; Leone, Michele. (Un)detectable Cluster Structure in Sparse Networks, Physical review letters, 101, 078701 (2008).
  * [25] Nadakuditi, Raj; Newman, M. Graph Spectra and the Detectability of Community Structure in Networks, Physical review letters, 108, 188701 (2012).
  * [26] Holland, Paul; Laskey, Kathryn; Leinhardt, Samuel. Stochastic Blockmodels: First Steps, Social Networks - SOC NETWORKS, 5, 109-137 (1983).
  * [27] Mossel, Elchanan; Neeman, Joe; Sly, Allan. Reconstruction and estimation in the planted partition model, Probability Theory and Related Fields, 162, (2014).
  * [28] Hubert, Lawrence; Arabie, Phipps. Comparing Partitions, Journal of Classification, 2, 193-218 (1985).
  * [29] Rossetti, Giulio; Cazabet, Remy. Community Discovery in Dynamic Networks: A Survey, ACM Computing Surveys, 51, (2017).
  * [30] Palla, Gergely; Barabasi, Albert-Laszlo; Vicsek, Tamás. Quantifying social group evolution, Nature, 446, 664-7 (2007).
  * [31] Dall’Amico, Lorenzo; Couillet, Romain; Tremblay, Nicolas. Community detection in sparse time-evolving graphs with a dynamical Bethe-Hessian, Advances in Neural Information Processing Systems, 33, (2020).
  * [32] Clauset, Aaron; Newman, M; Moore, Cristopher. Finding community structure in very large networks, Physical review. E, Statistical, nonlinear, and soft matter physics, 70, 066111 (2005).
  * [33] Krings, Gautier; Karsai, Márton; Bernharsson, Sebastian; Blondel, Vincent; Saramäki, Jari. Effects of time window size and placement on the structure of aggregated networks, EPJ Data Science, 1, (2012).
  * [34] Clauset, Aaron; Eagle, Nathan. Persistence and periodicity in a dynamic proximity network, arXiv, 12117343, (2012).
  * [35] Bassett, Danielle; Porter, Mason; Wymbs, Nicholas; Grafton, Scott; Carlson, Jean; Mucha, Peter. Robust Detection of Dynamic Community Structure in Networks. Chaos (Woodbury, N.Y.), 23, 013142 (2013).
  * [36] Bohrnstedt, George; Goldberger, Arthur. On the Exact Covariance of Products of Random Variables, Journal of the American Statistical Association, 64, 1439-1442 (1969).
  * [37] Akaike, H. A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716-723 (1974).

##  Appendix A Supplementary Information

###  A.1 Explicit Expressions for $p_{ij}$ and $q_{ij}$

The expressions for $p_{ij}$ and $q_{ij}$ are derived from the solution to the
problem defined by Hamiltonian in eq. (17), here recalled:

$H(\mathcal{G}|\vec{\alpha},\vec{\beta})\equiv\sum_{i=1}^{N}\left[\alpha_{i}\sum_{t=1}^{T}\frac{k_{i}(t)}{T}+\beta_{i}\sum_{t=1}^{T}\frac{h_{i}(t)}{T}\right]$
(42)

In [12] it is provided an exact formulation for these quantities.
Specifically, through appropriate reparametrization, the model outlined by the
Hamiltonian can be precisely mapped to a combination of non-interacting one-
dimensional Ising models. These models are then solvable analytically, with
their solutions expressible in terms of a constant connection probability.

$\displaystyle
p_{ij}=\left(\frac{x_{i}x_{j}y_{ij}-1}{2\sqrt{4x_{i}x_{j}+(x_{i}x_{j}y_{ij}-1)^{2}}}\right)\left(\frac{(\lambda_{ij}^{+})^{T}-(\lambda_{ij}^{-})^{T}}{(\lambda_{ij}^{+})^{T}+(\lambda_{ij}^{-})^{T}}\right)+\frac{1}{2}$
(43) $\displaystyle
q_{ij}(\tau)=p_{ij}^{2}+p_{ij}(1-p_{ij})\left(\frac{(\lambda_{ij}^{-})^{\tau}(\lambda_{ij}^{+})^{T-\tau}+(\lambda_{ij}^{+})^{\tau}(\lambda_{ij}^{-})^{T-\tau}}{(\lambda_{ij}^{+})^{T}+(\lambda_{ij}^{-})^{T}}\right)$

With

$\displaystyle\lambda_{ij}^{\pm}=e^{J_{ij}}\cosh
B_{ij}\pm\sqrt{e^{2J_{ij}}\sinh^{2}B_{ij}+e^{-2J_{ij}}}$ (45)

and

$\displaystyle B_{ij}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\left[\ln(x_{i})+\ln(x_{j})+\ln(y_{i})+\ln(y_{i})\right]$
(46) $\displaystyle J_{ij}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{4}\left[\ln(y_{i})+\ln(y_{j})\right]$ (47)

With, $x_{i}=e^{-\alpha_{i}/T}$, while $y_{i}=e^{-\beta_{i}/T}$.

Based on these terms, we can determine the probability of observing or not
observing a link at time $t+\tau$ given its presence or absence at time $t$.
This concept is then translated in the stochastic matrix defined in eq. (21).

###  A.2 Auto-covariance for the matrix B

Our main argument for the computation of the DT in the case in which
modularity using the persisting degree is employed exploits a heuristic, based
on the assumption that after $\tau^{*}$ define in equation (28) all the
matrices being in the sequence $\mathcal{B}$ and $\mathcal{A}$ can be
considered independent.

Using the argument in [36], we compute:

$\displaystyle Cov(B_{ij}(t),B_{ij}(t+\tau))=$ (48) $\displaystyle
Cov\left(a_{ij}(t)a_{ij}(t+\tau),a_{ij}(t+\tau)a_{ij}(t+2\tau)\right)$

Where, by conducting the calculations, we have the result as follows:

$\displaystyle Cov(B_{ij}(t),B_{ij}(t+\tau))=$ $\displaystyle
E[a_{ij}(t)]^{2}\big{[}Cov(a_{ij}(t),a_{ij}(t+\tau))^{2}+Var(a_{ij}(t))+Cov(a_{ij}(t),a_{ij}(t+2\tau))\big{]}+$
$\displaystyle+Cov(a_{ij}(t+\tau))^{2}+Var(a_{ij}(t))Cov(a_{ij}(t),a_{ij}(t+2\tau))$
(49)

Where we use the fact that
$Cov(a_{ij}(t),a_{ij}(t+\tau)=Cov(a_{ij}(t+\tau),a_{ij}(t+2\tau))$

That, after substituting the values of the covariances as computed in [12] and
the expected values, becomes:

$\displaystyle Cov(B_{ij}(t),B_{ij}(t+\tau))=$
$\displaystyle=p_{ij}^{2}\big{[}p_{ij}(1-p_{ij})+2g_{ij}(\tau)+g_{ij}(2\tau)\big{]}+g(\tau)^{2}+p_{ij}(1-p_{ij})g_{ij}(2\tau)=$
$\displaystyle=p_{ij}^{3}(1-p_{ij})+p_{ij}^{2}g_{ij}(2\tau)+2p_{ij}^{2}g(\tau)+g_{ij}(\tau)^{2}+p_{ij}(1-p_{ij})g_{ij}(2\tau)$
(50)

Where $g(\tau)$ is a function of $\tau$, defined in [12], whose value decays
with $\tau$:

$g(\tau)=p_{ij}(1-p_{ij})(\mu)^{\tau}$ (51)

with $\mu<1$. Considering 50, we note that after $\tau^{*}$, and considering
the term $p_{ij}^{3}(1-p_{ij})$ small enough, it is equal to 0.

###  A.3 Generating Community-Preserving Graph Trajectories via Markov Chains

Operationally, we begin by defining a network generated according to a
stochastic block model characterized by an internal cluster connection
probability, $p_{in}$, and an external connection probability, $p_{out}$. At
this point, two different stochastic matrices are used to define the evolution
of connections in this initial network. These matrices are characterized by
$p_{in}$ and $q_{in}$ for internal cluster connections, and $p_{out}$ and
$q_{out}$ for external connections. By varying the ratio between the densities
of link and link persistence, we are able to study the behavior of our
methodology.

Trajectories are generated in a way that preserves both the network’s average
total degree $\bar{k}^{*}$ and the average total persistent degree
$\bar{h}^{*}$.

The concept revolves around altering the ratio between the average
contribution from internal cluster connections and that from external links.
This goal is achieved by adjusting $p_{in}$ and $p_{out}$ as well as $q_{in}$
and $q_{out}$. To maintain $\bar{k}^{*}$ and $\bar{h}^{*}$ as mentioned, we
use the following formulas.

$\displaystyle
p_{in}\frac{N}{n_{groups}}+p_{out}\frac{N(n_{groups}-1)}{n_{groups}}=\overline{k}^{*}$
(52) $\displaystyle
q_{in}\frac{N}{n_{groups}}+q_{out}\frac{N(n_{groups}-1)}{n_{groups}}=\overline{h}^{*}$
(53)

Where $n_{groups}$ represents the number of clusters, $N$ is the total number
of nodes in the network.

Hence by choosing $p_{out}$ and $q_{out}$ as function of $p_{in}$ and
$q_{in}$:

$\displaystyle p_{out}=\epsilon p_{in}$ (54) $\displaystyle q_{out}=\eta
q_{in}$ (55)

Where $\eta$ and $\epsilon$ are two parameters that define the internal and
external link/persisting link density within the clusters.

We can vary the difference between the internal cluster and the external one,
by working on these two parameters. Using this further conditions we obtain
the following formulas:

$\displaystyle\left\\{\begin{array}[]{l}p_{out}=\epsilon p_{in}\\\
p_{in}=\frac{n_{groups}\overline{k}^{*}}{\left[1+\epsilon(n_{groups}-1)\right]N}\end{array}\right.$
(58)

$\displaystyle\left\\{\begin{array}[]{l}q_{out}=\eta q_{in}\\\
q_{in}=\frac{n_{groups}\overline{h}^{*}}{\left[1+\eta(n_{groups}-1)\right]N}\end{array}\right.$
(61)

By using equations (58) and (61), we ensure that the expected value on the
trajectory of $\left<\bar{h}\right>=\bar{h}^{*}$ and
$\left<\bar{k}\right>=\bar{k}^{*}$.

By changing $\eta$ and $\epsilon$, we can increase or decrease the disparity
between the internal and external link density. Specifically these parameters
can assume values in the following intervals:

$\epsilon\in\left\\{kv\mid k\in\mathbb{N},\,nv\leq 1\ \&\ nv>0\right\\}$ (62)

$\eta\in\left\\{r\epsilon\mid r\in(0,1)\right\\}$ (63)

In this way, the network’s evolution will preserve the condition whereby the
persistent degree within the group and outside the group will always be lower,
respectively, than the internal degree to the group and the external one.

In the numerical simulations presented in Section V, we focused on a scenario
where $n_{groups}=2$.

###  A.4 A Schematic Presentation of the Algorithm

Here we explain the steps followed to optimize the modularity function. The
algorithm aims to maximize the modularity function by changing the labels
associated with each element of the network. We exploit the method proposed by
Clauset and Moore in [32], by readapting it to use the Maximum Entropy null
models. In this sense we rewrite the modularity function in a general form, so
that the function for which we want to find the maximum is:

$Q(g)=\sum_{i,j}\left[C_{ij}-P_{ij}^{(null)}\right]\delta_{g(i)g(j)}$ (64)

Where depending on the features considered, $C_{ij}$ is:

$C_{ij}=\begin{cases}a_{ij}\quad&\text{Static Case}\\\
\overline{a_{ij}}\quad&\text{Dynamic Case No Memory}\\\
\overline{a_{ij}}+\overline{b_{ij}}\quad&\text{Dynamic Case With
Memory}\end{cases}$ (65)

While $P_{ij}^{(null)}$ assumes the form of the null models presented in the
main text. The implementation proposed in [32] works by defining a matrix
$\Delta Q$, whose elements are indicated with $\Delta Q_{vw}$ (defined below)
for each couple of communities. The size of this matrix is initially assumed
to be equal to the number of nodes, i.e. at each node is initially attached a
unique label. The size is then reduced according to the merging procedure, so
to reduce the number of communities. Hence, for our application, we define the
initial value of $\Delta Q_{vw}$ as:

$\Delta Q_{vw}=C_{vw}-P_{vw}^{(null)}\quad\text{for}\quad j,i\in N$ (66)

Then is selected the largest $\Delta Q_{vw}$, and if it is positive,
communities v and w are merged, the value of Q is incremented of $\Delta
Q_{vw}$ and the matrix $\Delta Q$ is updated according to the following rule:

$\Delta Q_{vk}^{new}=\Delta Q_{vk}+\Delta Q_{wk}$ (67)

The procedure is repeated until it is no more convenient to merge communities.

###  A.5 Analyzing the Evolution of Detected Communities and the Role of
Memory

The purpose of this section is to present additional analyses that support the
findings from the main sections. In the first subsection, we conduct numerical
simulations to demonstrate how selecting the appropriate time window for data
aggregation correlates with the observed memory in the data, as measured by
the normalized autocovariance. This selection is relevant in enhancing the
algorithm’s ability in detecting communities. In the second subsection, we
compare the performance of modularity with memory against modularity without
memory across various time windows, illustrating the distinct advantages of
each approach.

####  A.5.1 Numerical Simulations time windows

We carried out a series of numerical experiments to verify whether, in line
with the findings from Section VI, the point at which the normalized
autocovariance peaks also corresponds to a more discernible community
structure as detected by the memory-exploiting modularity. In our simulations,
we begin with a time-varying graph generated as described in Section A.3,
featuring a structure of 5 communities. From this graph, we create additional
structures by disaggregating and aggregating snapshots. The disaggregation
process involves dividing each snapshot from the original time-varying graph
into a number of sub-snapshots. Each sub-snapshot $G_{sub}$ comprises links
randomly selected from the original snapshot, ensuring that the sum of all
sub-graphs reconstructs the original, hence:

$\bigcup_{t=1}^{n}G_{sub}(t)=G^{original}$ (68)

Where the operator $\bigcup$ operates in such a way that if a link is present
in both snapshot $t_{1}$ and snapshot $t_{2}$, it is counted only once in the
aggregated graph.
Temporal order is maintained; for example, if the first snapshot is divided
into $n$ graphs, these will occupy the same position as the original snapshot
in the complete time-varying graph. The aggregation process functions by
summing multiple snapshots to form a new snapshot, which will contain all the
links from the combined structures. The procedure ensures that there will be
non-overlapping snapshots. Through this approach, we establish a controlled
setting with a natural scale at which memory is incorporated into the
generative process. This enables us to numerically check whether this scale
corresponds to the point where the normalized autocovariance is higher and
where the modularity, with known ground-truth communities, performs optimally.

Figure 18: Results from the numerical simulations described in Sec. (A.5.1)

In Fig (18) are shown the results from a series of numerical simulations
varying the values of $\Delta k$ and $\Delta h$ in a temporal network
featuring 5 ground truth communities. Each point pertains to an aggregation
state of the network, selecting progressively more aggregated temporal windows
as the value on the x-axis increases. The two charts display the Adjusted Rand
Index (ARI) value (top) for each aggregation window and the value of
normalized autocovariance (bottom). Points connected by the same line
represent the same network across different aggregation windows. The points
highlighted in red correspond to the peak aggregation window for the variable
of interest. The natural scale of aggregation is 3, which is the baseline from
which all other temporal networks were generated.

####  A.5.2 Real-Data results

Below, we present data related to the application of the model aggregated over
various temporal windows. This information facilitates a comparative analysis
of the results obtained using modularity with memory against those obtained
without memory.

If we were to use the no-memory model to identify potential structural breaks
and then conduct community detection, we would find that no breakpoints are
detected for any temporal window. In this regard, it’s possible to compare the
efficacy of the two models in explaining the data using the Akaike Information
Criterion (AIC) [37]. By doing so, considering the structural breaks
identified in the memory model and those detected by the no-memory model
(none), we can compare the two models by summing the AICs for each segment.
The result is a clear victory for the memory model, which consistently
explains the data better across all temporal windows, as indicated by lower
AIC values. When comparing the results of community detection while accounting
for the breaks identified by both the memory and no-memory models, we observe
that the memory model almost always outperforms the no-memory model.

Min | 1’ | 2’ | 3’ | 4’ | 5’ | 6’ | 7’ | 8’ | 9’ | 10’ | 11’ | 12’ | 13’ | 14’ | 15’ | 16’ | 17’ | 18’ | 19’ | 20’
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
ARI | 0.70 | 0.80 | 0.82 | 0.85 | 0.82 | 0.84 | 0.86 | 0.86 | 0.85 | 0.87 | 0.85 | 0.87 | 0.86 | 0.87 | 0.87 | 0.90 | 0.87 | 0.88 | 0.89 | 0.89
Table 1: The table displays the average value of ARI for the different
temporal aggregation windows for the modularity in eq. (25). The average is
computed by considering the value for each individual segment indicated by the
changing points.

In Table 1, we show the results of modularity with memory applied for
different time window. The ARI associated with modularity no-memory is always
0.76 because of the fact that with this model no breaks are found. In the
memory-less model, the data fed into the algorithm is simply a temporal
average of link presence. Evidently, in our scenario, this approach influences
community detection in a manner independent of the temporal window chosen,
further motivated by the absence of structural breaks.

This yields a value of ARI for the memory-less modularity where, only for
$W=1^{\prime}$, the outcome exceeds that of the counterpart with memory. This
suggests that at shorter temporal windows, there may be memory effects not
linked to membership. As noted, such effects can potentially compromise the
performance of the algorithm.

Figure 19: *

(a)

Figure 20: *

(b)

Figure 21: *

(c)

Figure 22: *

(d)

Figure 23: *

(e)

Figure 24: *

(g)

Figure 25: In these figures, we have a 3D representation of the average value
behavior over 20 instances for no-memory modularity (top) and memory
modularity (bottom) tested on graph trajectories created by varying $\Delta k$
and $\Delta h$. From left to right, we change the length of the trajectory
with values T = 15, T = 50, and T = 100.