On null models for temporal small-worldness in brain dynamics

The brain is always active, even when it is not performing a cognitive task. Therefore, it is interesting to investigate its fluctuations starting from this case: the resting-state condition (Karahanoğlu & Van De Ville, 2017; Maltbie, Yousefi, Zhang, Kashyap, & Keilholz, 2022; Preti, Bolton, & Van De Ville, 2017; Thompson, Brantefors, & Fransson, 2017). In this context, functional magnetic resonance imaging (fMRI) techniques allow us to study the underlying functional architecture of the brain and its temporal evolution. We can represent cerebral activity by extracting and manipulating the blood oxygen level–dependent (BOLD) signal from the scans.

To test certain hypotheses about observed functional connectivity, after modeling brain activity as a complex temporal network (Bahrami, Laurienti, Shappell, Dagenbach, & Simpson, 2022; Sizemore & Bassett, 2018; Thompson et al., 2017), it is important to compare it with an appropriate null model (Lurie et al., 2020; Váša & Mišić, 2022). A null model is a random statistical object that has certain properties in common with the empirical data under consideration and is used to evaluate whether the latter has noteworthy features or properties that cannot be attributed to randomness or other constraints.

In the complex temporal networks of brain activity just mentioned, each node corresponds to a brain region, and edges represent the presence of interactions greater than a certain threshold between two connected regions. As the graph evolves over time, edges may appear and disappear, meaning that the interaction between two brain regions jumps over and below a certain value. Complex networks can be found across different fields such as biology, sociology, epidemiology, and brain dynamics (da Fontoura Costa et al., 2011). They are graphs with nontrivial topological properties; two main properties that have come to define them are scale-freeness and small-worldness (Bassett & Bullmore, 2017; Chung & Lu, 2006; Tomasi, Shokri-Kojori, & Volkow, 2017).

In this paper, the common property of the null model and the real data we focus on is small-worldness, which expresses the efficiency of information exchange between nodes. It is defined as the ratio between two properties of the network: the clustering coefficient and the average path length. The clustering coefficient measures how many connections exist among a node’s neighbors, which reflects the local connectivity of the network. The average path length measures the average number of steps it takes to travel from one node to another, which reflects the global connectivity of the network. Given the fact that we are considering temporal brain networks, the temporal version of small-worldness can be interpreted as a measure of efficient communication among brain regions. We propose a new definition of temporal small-worldness S = C/L as the ratio of the temporal clustering coefficient C and the temporal path length L. The temporal clustering coefficient C is the time average of the clustering coefficients computed at each time step. The temporal path length L is the average of the fastest temporal paths between all pairs of nodes. S differs from S_SB, the temporal small-worldness definition of Sizemore and Bassett (2018), as it replaces the temporal version clustering coefficient with the temporal correlation coefficient TC, which measures the change between two consecutive time steps. In the Measures subsection we will discuss the reasons behind the adoption of this new definition.

Until recently, neuroscientists mainly analyzed the static version of brain networks (Bassett & Bullmore, 2017; Fornito, Zalesky, & Bullmore, 2016; Liao, Vasilakos, & He, 2017; Váša & Mišić, 2022), as the temporal extension of graph models is a recent research topic even in computer science. Currently, the temporal null models that have been considered in the literature are the randomized edge model and the randomly permuted times model, which are based on empirical data (Sizemore & Bassett, 2018). For consistency with the other models used in this work, we will refer to them as the random temporal edge (RTE) graph model and the random temporal permuted times (RTPT) graph model. In the RTE graph, each edge is rewired at each time step by changing one of its endpoints in the original graph. In the RTPT graph, the temporal structure is deteriorated by randomly permuting the time at which edges occur. However, as we show in this work, these random models do not reflect the typical structure of small-world networks, which exhibit short average path lengths and low clustering coefficients (Watts & Strogatz, 1998).

Lately, in the research community of real-world complex networks, the random hyperbolic (RH) graph model has become famous because it can exhibit both a high-taildegree distribution and small-worldness (Bode, Fountoulakis, & Muller, 2016; Fountoulakis, 2012; Fountoulakis & Müller, 2016; Krioukov, Papadopoulos, Kitsak, Vahdat, & Boguñá, 2010). It is a geometric model whose nodes are points that lie in the hyperbolic space. The latter cannot be embedded in Euclidean space, owing to the exponential growth of its volume. It also happens in the case of functional brain networks, which are best represented in the hyperbolic disc as shown by Whi, Ha, Kang, and Lee (2022). A temporal version of the RH graph model has recently been considered (Hartle, Papadopoulos, & Krioukov, 2021; von Looz & Meyerhenke, 2018).

We also consider two other temporal geometric graphs, the random temporal square (RTS) graph model and the random temporal torus (RTT) graph model as possible good candidates for null models. We compare them with the state of the art and with the networks extracted from the signals contained in the empirical data of 1,047 subjects from the WU-Minn Human Connectome Project (WU-Minn HCP 1200 Subjects Data Release Reference Manual, 2018), which we collected in a dataset that we made publicly available (Rossi, Deslauriers-Gauthier, & Natale, 2023). We test whether our random temporal models have small-worldness values close to those of empirical data at different connectivity thresholds. In particular, we show that the random temporal hyperbolic (RTH) graph model is a suitable null model with respect to the small-worldness property, especially when compared to the previous RTPT and RTE graph null models considered by Sizemore and Bassett (2018).

In this section, we describe the fMRI preprocessed data we started from and the process we implemented to extract the temporal brain networks from them (Rossi et al., 2023).

The data we used were taken from the Human Connectome Project (WU-Minn HCP 1200 Subjects Data Release Reference Manual, 2018), specifically, the 3 Tesla Siemens fMRI dataset of 1,047 selected resting-state subjects (rfMRI) scanned at the University of Washington or the University of Minnesota. The data are preprocessed and minimally preprocessed according, respectively, to Smith et al. (2013) and Glasser et al. (2013), and had artifacts removed using ICA-FIX (Griffanti et al., 2014; Salimi-Khorshidi et al., 2014). During the acquisition of the brain images, the patients would lie quietly in a darkened room with their eyes open, looking at a fixed bright point on a dark background, without performing any tasks. This is arguably the simplest case one can imagine, and it is natural to start the investigation from this point because even when the subject does nothing, brain regions are activated and the signal fluctuates. Understanding this scenario is important to proceed with the comparison between this state and more complicated situations, such as when the patient has some neural disease or is performing cognitive tasks involving memory and attention (Liao et al., 2017).

The rfMRI data were acquired in four runs of 14 min and 33 s each. The data were captured in both right-to-left (RL) and left-to-right (LR) phase encoding, with two scans per direction. We used only one of the two left-to-right scans. Regarding the temporal resolution of the data, the total number of time steps is 1,200 and the interval between two time steps is 0.72 s. Regarding the spatial resolution, there are 91 voxels in the x- and z-axis and 109 in the y-axis, making the total number of voxels 902,629. The dimension of a voxel is 2 × 2 × 2 mm.

The pipeline we implemented takes the BOLD signals of the rfMRI scan for each voxel of the brain and transforms them into temporal graphs (Van Dijk et al., 2010). We release them in the publicly available dataset, Rossi et al. (2023). The first step of the pipeline is to perform a linear regression on the movement parameters of the data, which reduces the contribution of the subject’s movements during the 14-min procedure. The data also contain some noise caused by the patient’s involuntary respiratory and cardiac rates. This problem is mitigated by applying a bandpass filter that isolates the BOLD fluctuations within a range of 0.01–0.08 Hz (Van Dijk et al., 2010). To obtain the division of the brain into regions that correspond to the nodes of our temporal networks, we use a brain atlas, such as the Gallardo et al. (2018), Glasser et al. (2016), or Schaefer et al. (2018). Using an atlas, we extract groups of voxels from gray matter areas that differ by anatomical and functional criteria. By averaging their signal, a time series is created for each brain region. Using the sliding window method, we individuate rectangular windows of length 60 s with an overlap of 30 s (Preti et al., 2017), each window being equivalent to one time step of the network. The Pearson correlation coefficient is calculated between each pair of brain region time series within each window; the resulting correlation value is the weight of the edge between the respective brain region nodes. In this way, we obtain an adjacency matrix for each window. The rows and columns of the matrix are ordered so that regions of the left hemisphere are grouped in the first block and regions of the right hemisphere are grouped in the second block (see the adjacency matrix in Figure 1). The two more correlated blocks (corresponding to the lighter pixels in the figure), upper left and lower right, represent strong intrahemispheric connectivity, while the two extra-diagonal blocks correspond to weaker interhemispheric connections. Three diagonals can be distinguished, the middle one representing the correlation of a region with itself. The other two sub-diagonals correspond to the connectivity of each region with its counterpart in the opposite hemisphere (Fornito et al., 2016). The values range between –1 and 1. If they are high (close to 1), the interpretation is clear: The two areas are synchronously activated, thus more likely to communicate. On the other hand, there is no consensus on the interpretation of negative indices. Given that many measures are defined for unweighted graphs, it is also common to threshold the data (Fornito et al., 2016). Currently, there is no consensus on the most appropriate threshold, so we analyze the data according to different positive thresholds applied to the Pearson correlation coefficient. The sequence of the obtained matrices defines the temporal brain networks.

In this paper, we consider five temporal random graph models. The first two models are the random temporal permuted times (RTPT) graph model and the random temporal edges (RTE) graph model (Sizemore & Bassett, 2018), which are constructed by randomizing empirical data; in particular, the RTPT graph is obtained by randomizing the time at which each edge appears, while the RTE graph is obtained by randomly reassigning an edge endpoint to a different node within each time step. The other three models are the random temporal square (RTS) graph model, the random temporal torus (RTT) graph model, and the random temporal hyperbolic (RTH) graph model, which do not depend on empirical data. They are synthetic graphs obtained by randomly placing nodes on different spaces: the square, the torus, and the hyperboloid, respectively. At each time step, the temporal graph is constructed by connecting the nodes depending on their position. From one time step to the next, the node positions are updated according to a movement model. Both the adjacency and their movement are computed according to the geometry of the underlying space. A summary of the five models and their parameters can be found in Table 1. Next, we provide their formal mathematical definitions. In the next section, we will compare them with the temporal brain networks extracted from data.

A temporal graph, denoted as G(V, E_T), is given by a set of nodes V and a set of temporal edges E_T. Each temporal edge is a triple (u, v, t) where u, v ∈ V are the endpoints of the edge and t ∈ [1, …, T] is the time at which the edge appears thus E_T = {(u, v, t) : u, v ∈ V, t ∈ [1, …, T]}.

The random temporal permuted times graph model (Sizemore & Bassett, 2018) is obtained by taking the temporal graph G(V, E_T) of a subject (see Figure 2A) and randomizing the time at which each edge appears (see Figure 2B). We construct a G(V, $ET′$⁠) from RTPT(G) model by sampling a permutation π uniformly at random and defining the set $ET′$ = {(u, v, π(t)) : (u, v, t) ∈ E_T}.

The random temporal edges graph model (Sizemore & Bassett, 2018) is obtained from the temporal graph G(V, E_T) of a subject (see Figure 2A) by reassigning one endpoint of an edge to another node in each window time step (see Figure 2C). We construct a G(V, $ET′$⁠) from RTE(G) by sampling a permutation σ uniformly at random of the set V and by setting (u, σ(v), t) ∈ E′ for each contact (u, v, t) ∈ E.

An RS(r) graph G(V, E) (Penrose, 2003) is obtained by distributing nodes independently and uniformly at random on a unit square [0, 1]² and connecting each pair of nodes whose Euclidean distance is less than radius r (see nodes contained in orange balls in Figure 3A). The set of edges is thus E = {(u, v) : (u, v) ∈ V × V, ‖p_u − p_v‖₂ ≤ r}, where p_u and p_v are the coordinates of the positions of the points u and v.

To obtain the temporal version of both the RS and RT graph models, which we will refer to them as RTS(r, v) and RTT(r, v) graph models, we can define a movement model according to which the positions of the points are updated and the distances are recomputed. The number of updates is set equal to the number of time steps of the real data. The new position of a point is determined by applying a displacement vector to it, whose direction is chosen uniformly in [0, 2π) and whose length is chosen uniformly in (0, v). The v is called the speed parameter. In the case of an RT graph, there is no problem of updating the point position, since there are no boundaries. In the case of the unit square, boundaries need to be taken into account. If a point happens to move outside the boundary of the unit square, its trajectory is reflected inwards, meaning that its position is updated as if it had bounced off the boundary.

To obtain a temporal version, which we will refer to it as the RTH(ζ, α, R, v) graph model, we update the point position as many times as the number of time steps of the real data we want to compare with. The point position update is chosen in such a way that guarantees that the marginal distribution of each time step is the same as the initial one. Specifically, the polar coordinate is updated by adding to the previous value a number chosen according to a uniform distribution over (0, v), and then computing the modulo 2π. For the radial coordinate, we add to the previous value a number chosen according to a uniform distribution over (−v, v), and we reflect the result in the interval [0, 1] to keep the distribution uniform in [0, 1] (Hartle et al., 2021). Thus, this parameter controls the movement of the points, and the higher the value of v, the more the points move and the lower the percentage of edges that remain the same between two consecutive time steps (see Figure 4C).

In this section, we introduce some measures to define a temporal version of small-worldness. The first static qualitative definition of small-worldness was introduced by Watts and Strogatz (1998) as the combination of a high clustering coefficient, as in lattice graphs, and a short average path length, as in random graphs. Humphries and Gurney (2008) combined the previous measures into a quantifiable ratio: clustering coefficient divided by average path length, normalized by the corresponding measure calculated on a random graph, formally $CgLg/CrandLrand$⁠. If the ratio is greater than 1, it indicates a significant small-world property of the network. The extension to temporal complex networks requires adapting the previous measures to capture the time evolution. The path length can naturally be extended to the average length of the fastest time-respecting path between each pair of nodes. Tang, Scellato, Musolesi, Mascolo, and Latora (2010) replaced the local clustering coefficient with the temporal correlation coefficient, which was also considered by Sizemore and Bassett (2018), which we will define in the next paragraph.

The clustering coefficient is defined as the ratio between the total number of closed triplets and the total number of triplets (open and closed). More specifically, closed triplets are sub-triangles in the graph that are counted three times, and open triplets are paths of length two, as illustrated in Figure 5.

In the static case, which corresponds to a graph that remains unchanged at each time step, the average path length is retrieved from the temporal path length. On the contrary, in the latter case, the temporal correlation coefficient does not coincide with the clustering coefficient. For this reason, the temporal correlation coefficient does not seem to properly generalize the clustering coefficient by Watts and Strogatz (1998) to the temporal setting, while our definition appears to be consistent with the latter. To illustrate the differences between the two measures, consider the following simple examples of temporal graphs with four time steps. In the first, shown in Figure 7A, we have a five-node star graph that does not change over time in which the temporal correlation coefficient is 1, while the temporal clustering coefficient is 0 because there are no closed triplets. In the second case, we have a six-node temporal graph in which a three-node triangle subgraph alternates between two configurations, as shown in Figure 7B. In this case, the temporal correlation coefficient is 0, while the temporal clustering coefficient is 0.5. As a consequence, in our experiments, we consider and compare two versions of small-worldness:

• The temporal small-worldness is the ratio between the clustering coefficient and the temporal path length, without normalization; formally,

  $S=CL.$
• The Sizemore and Bassett temporal small-worldness is the ratio between the temporal correlation coefficient and the temporal path length, without normalization; in formulas,

  $SSB=TCL.$

In this section, we provide some technical information, such as the details of the machine on which we run the experiments and the inputs of the pipeline, and we describe the optimization methods we used to choose the model parameters.

We run the pipeline on the Nef Cluster Computing Platform of Inria Sophia Antipolis. By setting the parameters described below, the average execution time of our pipeline for a subject is 250 s, on a 64-core node with the following technical characteristics:

• AMD dual-EPYC 7542 processors,

• 1024 GB of RAM.

The pipeline is implemented in Python. As explained in the Methods section, it transforms the BOLD signal into a temporal graph. The code, which is available at https://github.com/aurorarossi/TemporalBrainNetworksCode (Rossi, 2024), takes various parameters as input, which are described in the next sections.

Regarding the type of atlas, it is possible to choose among Gallardo et al. (2018), Schaefer et al. (2018), or Glasser et al. (2016). For the first two, it is possible to select the number of regions into which the brain is divided, starting from 102 nodes up to a maximum of 1,002 nodes. The Glasser atlas instead has 362 regions. In our experiment, we chose the Schaefer atlas with 302 nodes. Using an atlas with a higher number of regions could potentially introduce errors due to the small size of the regions, which could result in regions not matching across different subjects.

Several parameters describe the windows, such as the seconds corresponding to the window length and the window overlap, which we chose to be 60 and 30 s, respectively. Among the different types of windows, we chose to work with rectangular windows.

As correlation measure we chose the classical Pearson’s correlation.

Since there was no consensus on the most appropriate threshold, we consider different values. Specifically, we consider all the thresholds from 0.2 to 0.9 with a step of 0.05 and from 0.9 to 0.98 with a step of 0.02. The output is a binary temporal network. The latter is sparse for high threshold values and dense for low values.

All temporal graph generators and measures are written in the Julia programming language. In addition to optimizing the execution of serial code, we applied parallel programming techniques such as multithreading, which were exploited to make the scripts efficient. These scripts take as input 532 GB of data generated by the pipeline starting from 975 GB of raw data.

After creating the temporal graphs with the pipeline, our goal is to compare them with the random temporal models, which all depend on some parameters. The latter should be chosen to share as much as possible some common properties with the real data. In particular, we aim to have the same number of nodes, average degree, and similar temporal small-worldness values. We thus choose the set of parameters for which the random temporal models best match the aforementioned properties in the empirical data. In the RTPT and RTE graph models, the number of nodes and the average degree are inherited from the real data by construction and they do not depend on other parameters. For the RTS and RTT graph models, the number of nodes can be set directly and, by changing the connection parameter r, it is possible to obtain any average degree; instead, the velocity parameter v has to be optimized. For the RTH graph model, the number of nodes has to be set, the average degree is obtained by changing R, and the parameters to be optimized are α and v, while ζ is set to 1, because, as explained in the Methods section, we can fix ζ and vary just α. We used two different optimization methods to find the best values of the parameters.

The first finds a set of parameters for the three random temporal models that would fit the small-worldness values of the real data for some average degree values in the interval 0–170 (see Figure 8). To obtain different average degrees for the real data networks we use different threshold values; the higher the threshold value, the lower the average degree. For the other models, we mention how to set it in the previous paragraph.

The second optimization method is to select a set of parameters for each random model that minimizes the area between the temporal small-worldness curves of the empirical data and the random temporal model over all different average degrees. The curves are obtained by linear interpolation of the temporal small-worldness values at different average degrees (see Figure 9). For the RTH graph model, two different sets of parameters are shown, one optimizing each definition of temporal small-worldness. For the RTS and RTT graph models, only one set of parameters is shown, corresponding to the result of optimization with respect to our definition of temporal small-worldness.

This section discusses the result of all the temporal measures described in this work, applied to the processed data from the Human Connectome Project (WU-Minn HCP 1200 Subjects Data Release Reference Manual, 2018) and to the random temporal models. We have calculated the values over different average degrees, namely in the range 0–170. Those of the empirical data vary by changing the threshold values. Consequently, this happens also for the RTPT and RTE graph models. As for the RTS and RTT graph models, this property depends on the radius of the connection: The higher the radius, the higher the average degree. In the hyperbolic case, as shown in the section above, it depends on R, which is the connection radius.

The plots in Figure 8 show the results of the first optimization method. The real data of all 1,047 patients are shown in blue, the RTH graph model in red, and the RTT and RTS graph models in orange and green, respectively. Figure 8A shows that for any average degree value we can find a set of parameters for the RTH graph model that can approximate the real data value of our definition of temporal small-worldness, while this is not the case for RTT and RTS graph models. Figure 8B shows that for any average degree value we can find a set of parameters for all the random models that can approximate the real data value of the temporal small-worldness value of Sizemore and Bassett (2018).

In the plots of Figure 9, as above, the real data are shown in blue. RTT and RTS graph models are shown in orange and green, respectively. RTPT and RTE graph models are shown in gray and pink, respectively. There are also two RTH graph models with different parameters. The one with α = 0.65, v = 0.6, and ζ = 1, illustrated in violet, better matches the small-worldness S value of the real data, while the one with α = 0.8, v = 0.9, and ζ = 1, illustrated in red, better fits the definition of small-worldness S_SB, which is based on the temporal correlation coefficient TC. For all models, the lines are the medians and the shadows correspond to the values between the first and third quartiles. Looking at the random models, the median of the RTPT graph model always has higher values with respect to those of the RTE graph model for all measures except the temporal path length (see Figure 9E). Note that for RTT, RTS, and RTH graph models the shadows are barely visible because the variance is very small. The quantitative results describing the plots 9A and 9B are given in Table 2 and Table 3, respectively. These tables illustrate the comparison between the median of the empirical data and that of the models. This comparison is made in two ways: by measuring the area between the curves, as explained in the Optimization subsection, and by using different types of norms, namely minimum, maximum, mean absolute, and mean Euclidean distance norms. The latter are computed across the differences between the linear interpolation of the empirical data and the model data values of S and S_SB. These differences are computed on 16,861 equispaced points in the range [0, 170]. The minimum captures the smallest difference between the curves, while the maximum captures the largest. The mean absolute value norm considers the absolute value of all differences at the selected points and then calculates their average. Analogously, the mean Euclidean distance norm calculates the Euclidean distance between the curves at the selected points and then calculates their average. The best results are highlighted in bold and the second best are underlined. In Table 2 we can see that the RTH graph model with α = 0.65, v = 0.6 is the one that best approximates the small-worldness S value of the real data, and the second best one is the other RTH graph model. In Table 3 we can see that the RTH graph model with α = 0.8, v = 0.9 is the one that best approximates the small-worldness S_SB value of the real data, while the second best one is the RTPT graph model.

In all plots, the situation is more noisy in the sparse regime. For the sake of completeness, we did not remove this part even though these regimes are rarely considered in the literature in neuroscience.

The main focus of this work is to investigate a novel appropriate null model for testing hypotheses about brain dynamics. A good null model must have similar properties to the empirical data under consideration. In this context, the two most important properties present in complex networks such as those modeling brain activity are small-worldness and scale-freeness. In this work, we focus on the small-worldness property.

To transform the data into temporal networks, the empirical data are collected in time and processed using the sliding window method. The main challenge is to adapt all definitions of measures and models to their temporal versions. In particular, we provide a new definition of temporal small-worldness that is consistent with its static version: If the time steps of the temporal network are all identical, the temporal small-worldness is equal to the static small-worldness.

We compare five different temporal models by looking at how close their temporal small-worldness values are to the ones of the data at different thresholds. In Figure 8 we compare the model at each threshold individually, while in Figure 9 we compare them across different thresholds. The two models introduced by Sizemore and Bassett (2018) perform worse than the other three models for the S measure. This may be not surprising given the fact that they do not have any parameters to optimize but depend only on the data. The RTS and RTT graph models perform better than the other two models, and the measures computed on them are consistently very close to each other. This similarity may be explained by the fact that the two graphs look locally the same for all points that are not close to the border of the square or the torus, since the connectivity only differs in the latter area (see Figure 3). The RTH graph model is the one that not only achieves the closest values of the temporal small-worldness, but its curve of values follows the same trend as that of empirical data across different average degrees. This result is consistent with the fact that RTH is a popular complex network model that exhibits both small-worldness and scale-freeness.

This analysis could be similarly extended to the more challenging investigation of a temporal version of the scale-freeness property. It would be interesting to see whether the RTH graph model is still preferable. In this work, we focus on the resting-state fMRI in healthy patients. A further development could be the comparison with patients with diseases that affect the functionality of the whole brain. For example, we could consider diseases such as Alzheimer’s and Parkinson’s but also ADHD, depression, and schizophrenia. Again, the main question is whether the best model to fit the data remains the same across different parameter ranges. Another condition that could be investigated is the case of nonresting state, where patients are engaged in cognitive processes generated by tasks involving memory, learning, and attention. On the null models side, as future work, one could consider the case in which one starts from phase randomization and autoregressive randomization null data of brain dynamics, transforms them in temporal networks, and then applies the same analysis as in this work. This would allow a comparison of the results of the two different approaches (Liégeois, Laumann, Snyder, Zhou, & Yeo, 2017; Liégeois, Yeo, & Van De Ville, 2021).

This work has been supported by the French government, through the UCA DS4H Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-17-EURE-0004. The authors are grateful to the OPAL infrastructure from Université Côte d’Azur for providing resources and support.

Aurora Rossi: Data curation; Formal analysis; Investigation; Methodology; Software; Validation; Visualization; Writing – original draft; Writing – review & editing. Samuel Deslauriers-Gauthier: Data curation; Formal analysis; Investigation; Methodology; Supervision; Validation; Writing – review & editing. Emanuele Natale: Conceptualization; Formal analysis; Funding acquisition; Investigation; Methodology; Project administration; Resources; Supervision; Writing – original draft; Writing – review & editing.

Aurora Rossi, Agence Nationale de la Recherche (https://dx.doi.org/10.13039/501100001665), Award ID: ANR-17-EURE-0004.

Null model:

A theoretical model used to test hypotheses about empirical data, sharing certain properties with them.

Threshold:

A criterion value that determines whether values above it are taken into account and values below it are not.

Node degree:

Number of edges incident to a node.

High-tail distribution:

A distribution that presents a significant number of extreme events.

Temporal graph:

A graph in which edges and nodes change over time. It can be represented as a sequence of graphs, each corresponding to a time step of the network at a given time.

Brain atlas:

Subdivision of the gray matter of the brain into smaller areas that share a common property.

Adjacency matrix:

A square matrix representing the graph structure, where the indices of rows and columns represent the nodes and the entries correspond to the edge weight between the respective nodes.

Lattice graph:

A graph that contains a regular pattern.

Time-respecting path:

A finite or infinite sequence of distinct edges, ordered in time, connecting a set of distinct nodes.

Binary temporal graph:

A temporal graph whose edges have weights of 1 or 0.

https://github.com/aurorarossi/TemporalBrainNetworksCode