Aurora Rossi

Samuel Deslauriers-Gauthier

Emanuele Natale

On null models for temporal small-worldness in brain dynamics

Keywords: Brain Networks, Temporal Networks, Small-worldness, Hyperbolic Graph, Null Model, fMRI

Brain dynamics can be modeled as a temporal brain network starting from the activity of different brain regions in functional magnetic resonance imaging (fMRI) signals. When validating hypotheses about temporal networks, it is important to use an appropriate statistical null model that shares some features with the treated empirical data. The purpose of this work is to contribute to the theory of temporal null models for brain networks by introducing the Random Temporal Hyperbolic Graph Model (RTH), an extension of the Random Hyperbolic Graph (RH), known in the study of complex networks for its ability to reproduce crucial properties of real-world networks. We focus on temporal small-worldness which, in the static case, has been extensively studied in real-world complex networks and has been linked to the ability of brain networks to efficiently exchange information. We compare the RTH Graph Model with standard null models for temporal networks and show it is the null model that best reproduces the small-worldness of resting brain activity. This ability to reproduce fundamental features of real brain networks, while adding only a single parameter compared to classical models, suggests that the RTH Graph Model is a promising tool for validating hypotheses about temporal brain networks.

INTRODUCTION

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 [START_REF] Karahanoglu | Dynamics of large-scale fMRI networks: Deconstruct brain activity to build better models of brain function[END_REF][START_REF] Maltbie | Comparison of Resting-State Functional MRI Methods for Characterizing Brain Dynamics[END_REF][START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF][START_REF] Thompson | From static to temporal network theory: Applications to functional brain connectivity[END_REF]. 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-Oxygenation-Level Dependent (BOLD) signal from the scans.

To test certain hypotheses about observed functional connectivity, after modeling brain activity as a complex temporal network [START_REF] Bahrami | A mixed-modeling framework for whole-brain dynamic network analysis[END_REF][START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF][START_REF] Thompson | From static to temporal network theory: Applications to functional brain connectivity[END_REF], it is important to compare it with an appropriate null model [START_REF] Lurie | Questions and controversies in the study of time-varying functional connectivity in resting fMRI[END_REF][START_REF] Váša | Null models in network neuroscience[END_REF]. 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 [START_REF] Costa | Analyzing and modeling real-world phenomena with complex networks: a survey of applications[END_REF]. They are graphs with non-trivial topological properties, two main properties that have come to define them are scale-freeness and small-worldness [START_REF] Bassett | Small-World Brain Networks Revisited[END_REF]Chung & Lu, 2006;[START_REF] Tomasi | Brain Network Dynamics Adhere to a Power Law[END_REF].

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 [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF], as it replaces the temporal version clustering coefficient with the temporal correlation coefficient T C, which measures the change between two consecutive time steps. In the Measures Section we will discuss the reasons behind the adoption of this new definition.

Until recently, neuroscientists mainly analyzed the static version of brain networks [START_REF] Bassett | Small-World Brain Networks Revisited[END_REF][START_REF] Fornito | Fundamentals of brain network analysis[END_REF][START_REF] Liao | Small-world human brain networks: Perspectives and challenges[END_REF][START_REF] Váša | Null models in network neuroscience[END_REF], 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 [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF]. For consistency with the other models used in this work, we will refer to them as the Random Temporal Egde Graph (RTE) Model and the Random Temporal Permuted Times Graph (RTPT) 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 [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF].

Lately, in the research community of real-world complex networks, the Random Hyperbolic Graph Model (RH) has become famous because it can exhibit both a high-tail degree distribution and small-worldness (Bode, Fountoulakis, & Muller, 2014;[START_REF] Fountoulakis | On the evolution of random graphs on spaces of negative curvature[END_REF][START_REF] Fountoulakis | Law of large numbers for the largest component in a hyperbolic model of complex networks[END_REF][START_REF] Krioukov | Hyperbolic geometry of complex networks[END_REF]. It is a geometric model whose nodes are points that lie in the hyperbolic space. The latter cannot be embedded in Euclidean space, due to the exponential growth of its volume. It also happens in the case of functional brain networks, that are best represented in the hyperbolic disc as shown by [START_REF] Whi | Hyperbolic disc embedding of functional human brain connectomes using resting-state fMRI[END_REF]. A temporal version of the RH Graph Model has recently been considered [START_REF] Hartle | Dynamic Hidden-Variable Network Models[END_REF][START_REF] Von Looz | Updating Dynamic Random Hyperbolic Graphs in Sublinear Time[END_REF].

We also consider two other temporal geometric graphs, the Random Temporal Square Graph (RTS) Model and the Random Temporal Torus Graph (RTT) Model as possible good candidates for null models.

We compare them to the state of the art and to the networks extracted from the signals contained in the empirical data of 1047 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 [START_REF] Rossi | Temporal Brain Networks[END_REF]. 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 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 [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF].

METHODS

Extracting temporal brain networks

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 [START_REF] Rossi | Temporal Brain Networks[END_REF].

Data.

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 1047 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 [START_REF] Smith | Resting-state fMRI in the Human Connectome Project[END_REF] and [START_REF] Glasser | The minimal preprocessing pipelines for the Human Connectome Project[END_REF], and had artifacts removed using ICA-FIX [START_REF] Griffanti | ICA-based artefact removal and accelerated fMRI acquisition for improved resting state network imaging[END_REF][START_REF] Salimi-Khorshidi | Automatic denoising of functional MRI data: Combining independent component analysis and hierarchical fusion of classifiers[END_REF]. 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 [START_REF] Liao | Small-world human brain networks: Perspectives and challenges[END_REF].

The rfMRI data were acquired in four runs of 14 minutes and 33 seconds 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 1200 and the interval between two time steps is 0.72 seconds. 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 902629. The dimension of a voxel is 2 × 2 × 2 millimeters.

Processing.

The pipeline we implemented takes the BOLD signals of the rfMRI scan for each voxel of the brain and transforms them into temporal graphs [START_REF] Van Dijk | Intrinsic Functional Connectivity As a Tool For Human Connectomics: Theory, Properties, and Optimization[END_REF]. We release them in the publicly available dataset [START_REF] Rossi | Temporal Brain Networks[END_REF]. 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-minute procedure. The data also contains 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 [START_REF] Van Dijk | Intrinsic Functional Connectivity As a Tool For Human Connectomics: Theory, Properties, and Optimization[END_REF]. 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 [START_REF] Gallardo | Groupwise structural parcellation of the whole cortex: A logistic random effects model based approach[END_REF], [START_REF] Glasser | A multi-modal parcellation of human cerebral cortex[END_REF] or [START_REF] Schaefer | Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity mri[END_REF]. 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 seconds with an overlap of 30 seconds [START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF], each window being equivalent to a 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 [START_REF] Fornito | Fundamentals of brain network analysis[END_REF]. 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 [START_REF] Fornito | Fundamentals of brain network analysis[END_REF]. 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.

Random graph models

In this paper, we consider five temporal random graph models. The first two models are the Random Temporal Permuted Times Graph Model (RTPT) and the Random Temporal Edges Graph Model (RTE) [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF], 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 Graph Model (RTS), the Random Temporal Torus Graph Model (RTT), and the Random Temporal Hyperbolic Graph Model (RTH), 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 timestep 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 to the temporal brain networks extracted from data.

Temporal graph.

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 ]}. Random Temporal Permuted Times Graph Model (RTPT(G)).
The Random Temporal Permuted Times Graph Model [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF] 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, E T ) from RTPT(G) model by sampling a permutation π uniformly at random and defining the set

E T = {(u, v, π(t)) : (u, v, t) ∈ E T }. Random Temporal Edges Graph Model (RTE(G)).
The Random Temporal Edges Graph Model [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF] 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, E T ) 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.

Random Square Graph Model (RS(r)).

A RS(r) graph G(V, E) [START_REF] Penrose | Random Geometric Graphs[END_REF] is obtained by distributing nodes independently and uniformly at random on a unit square [0, 1] 2 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 Random Torus Graph Model (RT(r)).

E = {(u, v) : (u, v) ∈ V × V, ||p u -p v || 2 ≤
A RT(r) graph G(V, E) is obtained by distributing nodes independently and uniformly at random on a torus [0, 1] 2 and connecting each pair of nodes that are at distance less than radius r (see nodes contained in orange balls in Figure 3b). On the torus [0, 1] 2 , the distance between two nodes p A , p B with coordinates, respectively (x A , y A ), (x B , y B ), is defined as

d(p A , p B ) = 1 2 -1 2 -|x A -x B | 2 + 1 2 -1 2 -|y A -y B | 2 .
Temporal version of Random Square and Torus Models.

To obtain the temporal version of both the RS and RT graphs 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). v is called the speed parameter. In the case of 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. likely to be near the border of the disk, otherwise, they are more likely to be near the center (see Figure 4a). Two nodes u and v are connected if they are at a hyperbolic distance d H less than R, where, considering the polar coordinate of u = (r u , θ u ) and v = (r v , θ v ), d H is defined as

d H ((r u , θ u ), (r v , θ v )) = 1 ζ acosh(cosh ζr u cosh ζr v -sinh ζr u sinh ζr v cos(π -|π -|θ u -θ v ||)).
Thus, the parameter R affects the average degree of the graph. The higher the value of R, the larger the distances between the nodes and the lower the average degree (see Figure 4b). 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] [START_REF] Hartle | Dynamic Hidden-Variable Network Models[END_REF]. 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).

Measures

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. [START_REF] Humphries | Network 'Small-World-Ness': A Quantitative Method for Determining Canonical Network Equivalence[END_REF] combined the previous measures into a quantifiable ratio: L rand . 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 [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF], which we will define in the next paragraph. The clustering coefficient at a given time t is denoted by C(t) and is required to compute the temporal clustering coefficient.

C(t) = i,j,k A ij (t)A jk (t)A ki (t) i k i (k i -1) (1)
where k i = j A ij (t) and A ij (t) is the value of the binary adjacency tensor of the network at time t corresponding to nodes i and j. The formula (1) is equivalent to the one introduced by Watt and Strogatz: the clustering coefficient is the average over all nodes of the local clustering coefficient

C i = 2|{e jk : v j , v k ∈ N i , e jk ∈ E}| k i (k i -1) ,
C i measures the tendency of neighbors of a node to form a clique. Some examples are shown in Figure 6.

The temporal clustering coefficient C is the average of the clustering coefficient computed at each time step: Temporal path length.

C = 1 T t C(t).
The temporal path length is defined as the average of the fastest paths between all pairs of nodes. Paths are taken across time. Only one edge can be used at each time step, and pauses are allowed at nodes (i.e., an edge can be not selected to move forward in a time step). The temporal path length L is thus defined as follows

L = 1 N (N -1) i =j l ij ,
where l ij is the minimum time needed to get from i to j starting at time 0. If there is no path between two nodes, the distance is set to the maximum distance, which is the number of time steps. [START_REF] Fornito | Fundamentals of brain network analysis[END_REF] suggest other ways to address the latter case, such as considering only paths that exist or choosing a threshold that yields only connected graphs.

Temporal correlation coefficient.

The temporal correlation coefficient quantifies how much the temporal structure changes over time. In particular, it shows the tendency of a node's neighbors to remain connected to it over successive time steps. The higher the coefficient, the smaller the changes over time [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF]. The temporal correlation coefficient of a node i is defined as

T C i = 1 T -1 T -1 t=1 j A ij (t)A ij (t + 1) j A ij (t) j A ij (t + 1)
.

We can thus write the temporal correlation coefficient for the whole temporal network as

T C = 1 N i T C i .
Temporal small-worldness.

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 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 4 time steps. In the first, shown in Figure 7a, we have a 5-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 6-node temporal graph in which a 3-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:

(a) A 5-node star graph. Its temporal correlation coefficient is 1, while its temporal clustering coefficient is 0.

(b) A 6-node temporal graph. Its temporal correlation coefficient is 0, while its temporal clustering coefficient is 0.5. The temporal small-worldness is the ratio between the clustering coefficient and the temporal path length, without normalization; formally

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

S SB = T C L .

EXPERIMENT SETUP

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

Computing Platform

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 seconds, on a 64-core node with the following technical characteristics:

AMD dual-EPYC 7542 processors, 1024 GB of RAM.

Pipeline

The pipeline is implemented in Python. As explained in the Method Section, it transforms the BOLD signal into a temporal graph. The code, which is available at the following link https://github.com/aurorarossi/TemporalBrainNetworksCode, takes various parameters as input, which are described in the next sections.

Atlas.

Regarding the type of atlas, it is possible to choose among Gallardo (2018), [START_REF] Schaefer | Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity mri[END_REF] or [START_REF] Glasser | A multi-modal parcellation of human cerebral cortex[END_REF]. 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 1002 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.

Windows.

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 seconds, respectively. Among the different types of windows, we chose to work with rectangular windows.

Correlation.

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

Thresholds.

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.

Generators and measures

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 gigabytes of data generated by the pipeline starting from 975 gigabytes of raw data.

Optimization: fitting parameters

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, while 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.

RESULTS

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 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 1047 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 [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF].

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 grey and pink.

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 (a) (b)

Figure 8: Result of the first optimization method: comparison between empirical and random networks with 302 nodes at each average degree for the two definitions of small-worldness.

small-worldness S SB , which is based on the temporal correlation coefficient T C. 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 due to the fact that the variance is very small. The quantitative results describing the plots 9a and 9b are given in Table 2 andTable 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 16861 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 is 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.

Model

DISCUSSION AND CONCLUSIONS

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 is collected in time and processed using the window sliding 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.

Figure 1 :

 1 Figure 1: Steps of the pipeline to extract temporal brain networks from an fMRI signal for each voxel of the brain. The first step is linear regression and bandpass filtering of the data. The second step is the division of the brain into regions according to an atlas. The third step is the sliding window method, which individuates rectangular windows within which, in the fourth step, the Pearson correlation coefficient is computed between each pair of brain region time series. The last step is the thresholding of the data.

  (a) A temporal graph. (b) A RTPT graph. (c) A RTE graph.

Figure 2 :

 2 Figure 2: Example of 4-time steps illustrating how the RTPT and RTE graphs are obtained from a given temporal graph. The RTPT graph is obtained by randomizing the time each edge appears over different time steps. The RTE graph is obtained by reassigning one endpoint of an edge to a different node within each window time step.

  (a) A RS graph. (b) A RT graph.

Figure 3 :

 3 Figure 3: Showing the main difference between the RS and the RT graphs. In both cases, the nodes are connected if they are at a distance less than a radius r. In the RS graph the nodes are distributed on a square, hence the node on the bottom right and the node on the top left are not connected; on the torus, on the contrary, they are connected.

  Figure 3a). The set of edges is thusE = {(u, v) : (u, v) ∈ V × V, ||p u -p v || 2 ≤r} where p u and p v are the coordinates of the positions of the points u and v.

Random

  Hyperbolic Graph Model (RH(ζ, α, R)). Krioukov's method for generating a Random Hyperbolic graph is to distribute nodes quasi-uniformly within a disk of radius R centered on the upper half of a hyperboloid of given negative curvature K = -ζ 2 (Krioukov et al., 2010). Writing the node positions in polar coordinates (r, θ), where r ∈ [0, R] and θ ∈ [0, 2π], we have that the angular density is ρ(θ) = 1 2π and the radial coordinate density is ρ(r) = α sinh(αr) cosh(αR)-1 . The α parameter controls the spread of the point positions. If α = ζ, the points follow a uniform distribution. If α > ζ, the points are more

  As for the curvature parameter ζ, we set it to ζ = 1, since (Lemma 1.1 of Bode et al. (2014)) if two RH graphs have the same ratio ζ α and the other parameters equal then they produce the same distribution on graphs.

  (a) Representations of the RTH graph on the Poincaré disk for different values of the parameter α. (b) Example showing how the parameter R affects the average degree of the RTH graph when the other parameters are fixed α = 1.2, ζ = 1 and v = 0.9. (c) Example showing how the parameter v affects the percentage of edges that remain the same between two consecutive time steps when the other parameters are fixed α = 0.5, R = 4.5 and ζ = 1.

Figure 4 :

 4 Figure 4: Behavior of the RTH Graph Model parameters.

  Temporal clustering coefficient. 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.

Figure 5 :

 5 Figure 5: Open and closed triplets. The open triplet in this graph is the green path, instead, the orange triangle is a closed triplet.

Figure 6 :

 6 Figure 6: Local clustering of the orange node. In (a) the clustering coefficient is equal to zero because there are zero connections among the neighbors of the orange node, over three possible ones. In (b) there is one connection over three. In (c) there are all possible connections, hence the clustering coefficient has value 1, which is the maximum.

Figure 7 :

 7 Figure 7: Small examples of temporal graphs where the temporal correlation coefficient differs from the temporal clustering coefficient in a significant way.

Figure 9 :

 9 Figure 9: Results of the second optimization method: comparison between empirical and random networks with 302 nodes for different measures.-23-

Table 1 :

 1 Summary of the random temporal graph models described in the Random graph section.

	Random Temporal Permuted Times	RTPT	Data-based	-	G	Sizemore and Bassett (2018)
	Random Temporal Edges	RTE	Data-based	-	G	Sizemore and Bassett (2018)
	Random Temporal Square	RTS	Synthetic	Square [0, 1] 2	r, v	Penrose (2003)
	Random Temporal Torus	RTT	Synthetic	Torus [0, 1] 2	r, v	Penrose (2003)
	Random Temporal Hyperbolic	RTH	Synthetic	Hyperboloid with ζ < 0 ζ, α, R, v Krioukov et al. (2010),
						von Looz and Meyerhenke (2018)
	clustering coefficient divided by average path length, normalized by the corresponding measure
	Lg calculated on a random graph, formally Cg	C rand

Table 2 :

 2 Area Minimum Maximum Mean Absolute value Mean Euclidean distance Distances between the median of the empirical data and the one of the model for the temporal small-worldness S, as discussed in Result Section. The best results are highlighted in bold and the second best is underlined.

	RTH α = 0.8, v = 0.9	2.84 0.0000014	0.067	0.017	0.00018
	RTH α = 0.65, v = 0.6 2.73 0.0000011	0.040	0.016	0.00015
	RTS v = 0.1	8.02 0.0052702	0.166	0.048	0.00053
	RTT v = 0.2	9.88 0.0000023	0.169	0.059	0.00057
	RTE	33.92 0.0115314	0.248	0.201	0.00160
	RTPT	21.42 0.0128656	0.166	0.127	0.00102

Table 3 :

 3 Distances between the median of the empirical data and the one of the model for the temporal small-worldness S SB , as discussed in Result Section. The best results are highlighted in bold and the second best is underlined.

	Model	Area Minimum Maximum Mean Absolute value Mean Euclidean distance
	RTH α = 0.8, v = 0.9	4.46 0.0000006	0.067	0.026	0.00024
	RTH α = 0.65, v = 0.6 9.60 0.0093862	0.096	0.057	0.00050
	RTS v = 0.1	37.23 0.0000018	0.266	0.221	0.00174
	RTT v = 0.2	17.91 0.0136169	0.130	0.106	0.00084
	RTE	20.56 0.0396405	0.142	0.122	0.00095
	RTPT	7.25 0.0000087	0.060	0.043	0.00035

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.

-24-Authors: Aurora Rossi, Samuel Deslauriers-Gauthier and Emanuele Natale ACKNOWLEDGMENTS 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.

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 if 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 non-resting 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, transform them in temporal networks and then apply the same analysis as in this work. This would allow to compare the results of the two different approaches [START_REF] Liégeois | Interpreting temporal fluctuations in resting-state functional connectivity MRI[END_REF][START_REF] Liégeois | Interpreting null models of resting-state functional MRI dynamics: not throwing the model out with the hypothesis[END_REF]. Node degree: number of edges incident to a node.

High tail distribution: a distribution that presents a significant number of extreme events.

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.

Time-respecting path: a finite or infinite sequence of distinct edges, ordered in time, connecting a set of distinct nodes.

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

Binary temporal graph: a temporal graph whose edges have weights of 1 or 0.


Claude Pasquier
email: claude.pasquier@unice.fr

Frédéric Flouvat

Jérémy Sanhes

Nazha Selmaoui-Folcher

Attributed graph mining in the presence of automorphism

Keywords: Attributed graph, Frequent pattern mining, Automorphism, Structure mining, Itemset mining

Attributed directed graphs are directed graphs in which nodes are associated with sets of attributes. Many data from the real world can be naturally represented by this type of structure, but few algorithms are able to directly handle these complex graphs. Mining attributed graphs is a difficult task because it requires combining the exploration of the graph structure with the identification of frequent itemsets. In addition, due to the combinatorics on itemsets, subgraph isomorphisms (which have a significant impact on performances) are much more numerous than in labeled graphs. In this paper, we present a new data mining method that can extract frequent patterns from one or more directed attributed graphs. We show how to reduce the combinatorial explosion induced by subgraph isomorphisms thanks to an appropriate processing of automorphic patterns.

Introduction

Directed graphs are well suited to model complex structures present in the real world. Gene regulatory networks, for example, are directed graphs in which nodes represent genes and edges represent regulatory influences. The World Wide Web is well modeled by a directed graph in which nodes are pages and edges are hyperlinks. Email communications; social graphs in which an individual can follow the activities of other people; citation networks in which an article cites other articles are also ideally represented by directed graphs.

Because of these numerous applications, graphs are extensively studied in graph theory, and, more recently, in the context of data mining. Most researches focus either on unlabeled graphs or on graphs whose nodes are associated with a unique label. In a social graph, for example, labels can be the ID of individuals; in a gene regulatory network, they can represent genes' names; in a Web graph, nodes are often labeled with pages' URL, etc.

However, in many applications, objects (represented by nodes) can be associated with many features. Individuals in social graphs have many characteristics, as genes in regulatory networks. In citation networks, articles are associated with many data like keywords, authors, publication dates and patents. In these data, each characteristic associated with an object is an attribute of the corresponding node. Graphs in which nodes are annotated with sets of attributes (or itemsets) are named attributed graphs, and up to now, only few studies are devoted to their analysis. 1One of the classical tasks when analyzing graph data is the discovery of frequent subgraphs. The interest of such patterns is to highlight how node labels are often organized. Mining attributed graphs allows the identification of structural patterns, but also highlights the relationship between nodes' attributes.

Two main reasons make the mining of attributed graphs very difficult. First, it is necessary to combine the exploration of the graph structure with the identification of frequent itemsets associated with nodes. The second reason is that, as for labeled graphs, the performance of the mining process is highly affected by the cost of subgraph isomorphism tests [START_REF] Yan | gspan: graph-based substructure pattern mining[END_REF]. As noted by [START_REF] Yan | gspan: graph-based substructure pattern mining[END_REF], in sparse labeled graphs with diverse labels, the number of subgraph isomorphisms remains pretty small. However, it is no longer the case in attributed graphs where the combinatorics of attributes often implies a large number of subgraph isomorphism tests.

The number of subgraphs with automorphisms is also much higher. As highlighted by [START_REF] Kuramochi | An efficient algorithm for discovering frequent subgraphs[END_REF], the presence of automorphism is problematic for all labeled graph mining algorithms. Using algorithms based on pattern extension, each automorphism in a subgraph of size k can lead to a different candidate pattern of size k + 1. In the worst case, when the processed subgraph is an unlabeled clique (or a clique with the same label attached to all nodes), the number of automorphisms is k!.

However, up to now, little attention is given to this problem. Indeed, in labeled graphs, nodes are associated with unique labels, which are, in general, manifold. Therefore, it is not very frequent to have the same label associated with many nodes in a subgraph. This is different in attributed graphs because entities represented by nodes are characterized by multiple attributes, and some of them are very common. Thus, it greatly increases the number of possible automorphisms. In a social graph, for example, a given individual is often connected with a group of other individuals sharing several characteristics (age group, hobbies, social class, musical taste, etc.).

The key contributions of our work are the following: [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF] We present the problem of mining frequent attributed subgraphs in one or more directed attributed graphs, and (2) we show how the identification of frequent patterns can be achieved by exploring the spanning trees of the graphs. (3) We define a new canonical form adapted to the exploration of the search space, and (4) we show how to extend such patterns to generate attributed subgraphs. [START_REF] Borgelt | Mining molecular fragments: finding relevant substructures of molecules[END_REF] We analyze the particular case of automorphic patterns and show how the combinatorial explosion caused by subgraph isomorphisms can be significantly limited. Before concluding, [START_REF] Bringmann | What is frequent in a single graph?[END_REF] we present the results obtained by analyzing several artificial and real datasets with the program AADAGE (Automorphism Aware Directed Attributed Graph Explorer) that we have developed.

Related works

Many algorithms have been proposed for mining labeled graphs (with one item associated with each node). Almost all proposed solutions are adaptations of frequent itemset mining. They explore the solution space by enlarging candidate subgraphs and eliminating infrequent patterns [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF][START_REF] Jiang | A survey of frequent subgraph mining algorithms[END_REF]. These methods differ in their way to explore the search space and in their pruning strategies to avoid the generation of redundant solutions. Indeed, unlike itemset mining where it is easy to generate all possible solutions in a unique way, in graph mining, it is very common to build the same pattern several times.

One of the privileged ways to avoid repeatedly exploring the same part of the search space is to define a systematic and orderly manner to add edges and nodes to candidate patterns based on a canonical representation of subgraphs. The best-known methods for generating all candidate subgraphs are depth-first search (gSpan [START_REF] Yan | gspan: graph-based substructure pattern mining[END_REF], Closegraph [START_REF] Yan | CloseGraph: mining closed frequent graph patterns[END_REF]) and breadth-first search (MoSS/MoFa [START_REF] Borgelt | Mining molecular fragments: finding relevant substructures of molecules[END_REF]).

The underlying idea behind the use of canonical forms is to extend only one solution in each isomorphism class of subgraphs. In graph theory, there are many ways to define the canonical form of a graph. A commonly used method is to represent a graph as an adjacency matrix, to swap rows and columns in order to generate all isomorphic graphs and to choose one of the matrices as the canonical form. In data mining, one avoids as much as possible to perform complete (and costly) isomorphism tests for each new candidate. The canonical representation is chosen such that it is easy (and effective) to test whether a pattern has already been generated. The most commonly used solution is to create a sequence of all edges contained in a subgraph in the order in which they were added [START_REF] Borgelt | Canonical forms for frequent graph mining[END_REF][START_REF] Yan | gspan: graph-based substructure pattern mining[END_REF].

The use of canonical forms, if it avoids exploring several times the same subgraph, does not obviate the need to identify all the possible ways to generate patterns. Thus, if the nodes are ordered in the lexicographical order of their labels, the pattern in Fig. 1b is canonical, and it appears 4 times in the initial graph presented in Fig. 1a. The patterns of Fig. 1c,d, which are also canonical, can be generated in 12 different ways (in the initial graph, there are 12 possible ways to choose 2 b nodes out of the 4 children labeled b or a). It is necessary to consider all these forms because the chosen configuration influences the future extensions of the pattern; for example, whether or not to extend each of the nodes labeled b with a node labeled c. The patterns with many subgraph isomorphisms in input data are challenging for all existing algorithms because of the problem of subgraph isomorphism, which is itself NPcomplete. Patterns presented in Fig. 1c,d, which have automorphisms, generate the biggest problems and contribute significantly to performance degradation of the algorithms, because the number of subgraph isomorphisms with the initial graph is increased.

(d)

In the case of sparse labeled graphs, this number remains relatively small. In attributed graphs, however, since many attributes are associated with nodes, the probability of observing subgraphs sharing one or more identical attributes is greatly increased. One way to circumvent it is to remove from the graph the most frequent items, but we are depriving ourselves of potentially interesting patterns.

Some studies deal specifically with the analysis of graphs in which the nodes are associated with itemsets as in the works of [START_REF] Fukuzaki | Finding itemset-sharing patterns in a large itemsetassociated graph[END_REF] or [START_REF] Miyoshi | Frequent pattern discovery from a single graph with quantitative itemsets[END_REF]. The issues addressed in these studies are different from ours in the sense that they deal with the analysis of undirected graphs and search for subgraphs sharing the same itemsets.

In our work, we seek to identify the recurrent presence of certain attributes associated with connected nodes. In its objectives, our work is closer to sequence mining [START_REF] Agrawal | Mining sequential patterns[END_REF][START_REF] Ayres | Sequential pattern mining using a bitmap representation[END_REF] or attributed tree mining [START_REF] Pasquier | Frequent pattern mining in attributed trees: algorithms and applications[END_REF]. In our previous work [START_REF] Pasquier | Frequent pattern mining in attributed trees: algorithms and applications[END_REF], we proposed a method to mine substructures reflecting the evolution of itemsets. However, our algorithm was dedicated to the analysis of datasets that can be represented as attributed trees.

Basic concepts and problem statement

In this section, we give basic definitions and concepts, and then we introduce the problem of attributed graph mining.

Preliminaries

Let I = {i 1 , i 2 , .., i n } be a set of items. An itemset is a subset of I. The items belonging to an itemset are sorted by lexicographical order and are accessed by their index (e.g., I 2 to access the second item of I).

A directed attributed graph G = (V, E, λ) on a set of items I consists of a set of vertices or nodes V , a set of edges E ⊆ {(u, v) ∈ V 2 |u = v} and a labeling function λ(v) : V → P(I) which associates a set of labels i ∈ I to each vertex v ∈ V . P(I) denotes the power set of I. For an edge (u, v) ∈ E, u is the parent of v and v is the child of u. There is a cycle in the directed graph if a path can be found from a vertex to itself. A directed labeled graph is rooted if one can find a vertex v such that there exists a path between v and any other vertex of the graph.

Two attributed graphs

G 1 = (V 1 , E 1 , λ 1 ) and G 2 = (V 2 , E 2 , λ 2 ) are isomorphic iff there exists a bijective function ϕ : V 1 → V 2 such that: 1. ∀v ∈ V 1 : λ 1 (v) = λ 2 (ϕ(v)), 2. ∀(u, v) ∈ E 1 ⇔ (ϕ(u), ϕ(v)) ∈ E 2 .
An automorphism is an isomorphism where G 1 = G 2 . Attributed graphs can be seen as itemsets organized in a graph structure. As such, the notion of attributed subgraph can be defined w.r.t. itemsets inclusion and structural inclusion. For itemset inclusion, we say that an attributed graph G 1 is contained in another attributed graph G if both attributed graphs have the same structure, and for each vertex of G 1 , the associated itemset is contained in the itemset of the corresponding vertex in G.

Definition 1 (Itemset/content inclusion

) G 1 = (V 1 , E 1 , λ 1 ) is contained in G = (V, E, λ) and is denoted by G 1 I G, iff there exists a function ϕ : V 1 → V such that: 1. ∀v ∈ V 1 : λ 1 (v) ⊆ λ 2 (ϕ(v)), 2. ∀(u, v) ∈ E 1 ⇔ (ϕ(u), ϕ(v)) ∈ E 2 .
Structural inclusion is represented by the classical concept of subgraph [START_REF] Inokuchi | An apriori-based algorithm for mining frequent substructures from graph data[END_REF][START_REF] Kuramochi | Frequent subgraph discovery[END_REF][START_REF] Yan | Graph indexing: a frequent structure-based approach[END_REF]. In labeled graph mining, the labels on vertices are preserved by the mapping function. In the case of attributed graphs, the mapping function preserves the itemsets associated with vertices.

Definition 2 (Structural inclusion) G

1 = (V 1 , E 1 , λ 1 ) is included in G = (V, E, λ) and is denoted by G 1 S G, iff there exists a function ϕ : V 1 → V such that: 1. ∀v ∈ V 1 : ϕ(v) ∈ V , 2. ∀v ∈ V 1 : λ 1 (v) = λ(ϕ(v)), 3. ∀(u, v) ∈ E 1 : (ϕ(u), ϕ(v)) ∈ E.
An attributed subgraph is defined in the following way to include both content and structural inclusion.

Definition 3 (Attributed subgraph) An attributed graph G 1 = (V 1 , E 1 , λ 1 ) is an attributed subgraph of G = (V, E, λ) and is denoted by G 1 G, iff there exists a function ϕ : V 1 → V such that: 1. ∀v ∈ V 1 : ϕ(v) ∈ V , 2. ∀v ∈ V 1 : λ 1 (v) ⊆ λ(ϕ(v)), and 3. ∀(u, v) ∈ E 1 : (ϕ(u), ϕ(v)) ∈ E. If G 1 is an attributed subgraph of G, we say that G is an attributed supergraph of G 1 . Each attributed graph G 1 that is isomorphic to G defines one embedding of G 1 in G.
Two different embeddings may refer to the same nodes and edges, simply by mapping the nodes in a permuted way.

Problem statement

Frequent attributed subgraph mining

Our work considers two settings w.r.t. input data: the transactional setting and the single graph setting. In the first case, input data are organized as a transactional database B in which each transaction is a graph. In such case, it is possible to calculate, for each generated pattern P, a per-transaction frequency, which is given by the number of graphs in B for which P is a subgraph. In the second case, input data are only composed of a single graph (not necessarily connected). In such case, we use the frequency measure defined by [START_REF] Bringmann | What is frequent in a single graph?[END_REF]. Using this definition, the absolute frequency of a graph pattern P is equal to the minimum number of unique nodes in G that a node of P is mapped to. Based on this measure, we can calculate a relative frequency, which is given by the absolute frequency divided by the number of nodes in the graph.

These two methods of frequency calculation are used the same way. In this paper, we will refer to these measures by the generic term "frequency."

A pattern is frequent if its frequency is greater than or equal to a minimum threshold value that is called support. The problem consists in enumerating all frequent attributed subgraphs in a given dataset.

Mining closed attributed subgraphs

The problem with frequent attributed graph mining is that the number of frequent patterns is often huge. In real applications, generating all solutions can be very expensive or even impossible. Moreover, many of these frequent attributed subgraphs contain redundant information.

Huge efforts have been made to design condensed representations that are able to summarize solutions in smaller sets (see, e.g., the works of [START_REF] Mannila | Multiple uses of frequent sets and condensed representations[END_REF] and [START_REF] Pasquier | Discovering frequent closed itemsets for association rules[END_REF]). The set of closed patterns is an example of such a condensed representation [START_REF] Pasquier | Discovering frequent closed itemsets for association rules[END_REF]. We say that an attributed graph G is a closed attributed graph if none of its proper attributed supergraphs has the same support as G. However, mining frequent closed attributed graphs could also be very expensive. We propose to relax the closure property on attributed graphs and to mine only attributed graphs which are closed w.r.t. their content (i.e., itemsets associated with vertices). This condensed representation, called the set of c-closed attributed graphs (content closed), is a superset of the set of closed attributed graphs (and a subset of the set of attributed graphs). We say that an attributed graph G is a c-closed attributed graph if there is no pattern with the same structure, the same support and supersets in the corresponding vertices.

Canonical form of directed attributed graphs

We propose a canonical form based on the spanning tree of the pattern (a subgraph). Our canonical form is the sequence of nodes obtained by performing a depth-first traversal of the spanning tree (note that other approaches use a sequence of edges).

In order to unambiguously identify each spanning tree, it is sufficient to identify each node by its associated attributes and its depth. As we also need to represent reentrant links, we use a third attribute that contains the index of the destination node if it is already present in the pattern.

Ordering of itemsets

In our application, nodes are associated with itemsets and one should define a total order on itemsets. Given two itemsets I and J (I = J ), we say that

I < J iff ∃k ∈ [1, min(|I|, |J |)] such that (1) ∀i < k : I i = J i and (2) I k < J k or k = |J |.

Ordering of nodes

The code of a node is a triple (d, Q, p) where d is the depth of the node in the spanning tree, Q is the set of items associated with the node, and p is 0 if the node is used for the first time or p is equal to the position of a node in the sequence (starting position is 0) if the destination node is already in the sequence. To compare two nodes, it is sufficient to compare their codes as triples.

Given two triples

T 1 = (d 1 , Q 1 , p 1 ) and T 2 = (d 2 , Q 2 , p 2 )
, we say that T 1 < T 2 iff one of the following assertions is true:

1. d 1 > d 2 , 2. d 1 = d 2 and Q 1 < Q 2 , 3. d 1 = d 2 and Q 1 = Q 2 and p 1 < p 2 .
In Fig. 2, for example, the three children of node a in pattern P 1 are represented by the following triples: (1, cd, 0), (1, cd, 3) and (1, cd, 0). All nodes are located at a depth of 1, they are all associated with the same itemset cd, but the second node references the fourth node present in the sequence (whose index is 3).

Ordering of attributed subgraphs

Let P be a pattern with a root node r (P), we note code(r (P)) the code of this pattern. It is constructed by concatenating the codes of the nodes of the corresponding spanning tree explored in a depth-first fashion. The smallest code represents the canonical form.

Given two patterns P 1 and P 2 , respectively, encoded by

C 1 =< c 1 1 , ..., c m 1 > and C 2 =< c 1 2 , ..., c n 2 >, we say that C 1 < C 2 iff one of the following assertions is true: 1. ∃k ∈ [1, min(m, n)] such that ∀i < k : c i 1 = c i 2 and c k 1 < c k 2 2. ∀i < min(m, n) : c i 1 = c i 2 and m > n.
Figure 2 shows four isomorphic patterns of the same directed attributed graph. Their codes are represented by the following sequences:

1. code(r (P 1 )) = < (0, a, 0)(1, cd, 0)(2, cd, 0)(2, cd, 0)(1, cd, 3)(1, cd, 0)(2, cd, 0) > 2. code(r (P 2 )) = < (0, a, 0)(1, cd, 0)(2, cd, 0)(2, cd, 0)(1, cd, 2)(1, cd, 0)(2, cd, 0) > 3. code(r (P 3 )) = < (0, a, 0)(1, cd, 0)(2, cd, 0)(2, cd, 0)(1, cd, 0)(2, cd, 0)(1, cd, 3) > 4. code(r (P 4 )) = < (0, a, 0)(1, cd, 0)(2, cd, 0)(2, cd, 0)(1, cd, 0)(2, cd, 0)(1, cd, 2) >
By using the ordering of triples previously defined, we found code(r (P 4 )) < code(r (P 3 )) < code(r (P 2 )) < code(r (P 1 )); then, the pattern P 4 is the canonical form.

Enumeration of patterns

Subgraph patterns are generated by performing a depth-first traversal of the search space based on their spanning tree. All items present in the input graph(s) constitute the initial patterns that will be progressively extended. Each pattern extension is built using the canonical form sequence generated in the previous step.

The use of the strategy defined by [START_REF] Pasquier | Frequent pattern mining in attributed trees: algorithms and applications[END_REF] allows enumerating all spanning trees. The generation of new patterns is performed by using the rightmost path extension method [START_REF] Chi | Cmtreeminer: mining both closed and maximal frequent subtrees[END_REF]. Two types of extensions are possible: The itemset extension adds a new item to the itemset associated with the rightmost node of the spanning tree, and the structural extension adds a new child to one of the nodes composing the right path of the spanning tree. The method is complete but might generate redundant patterns in the form of isomorphic graphs. Duplicate candidates are detected and discarded before the candidate extension process by performing a canonical check. The limitation of this approach is that it only allows finding rooted patterns.

Identification of canonical and automorphic patterns

The code for the graphs defined above has, like most of the codes used by graph mining algorithms, the following property: Each prefix of a canonical code is itself canonical [START_REF] Huan | Efficient mining of frequent subgraphs in the presence of isomorphism[END_REF][START_REF] Yan | gspan: graph-based substructure pattern mining[END_REF]. The extension of a pattern does not change the order of the triples in the case of the generation of a new canonical pattern. Thus, each canonical pattern can only be obtained by adding a new triple or adding an item to the itemset associated with the last triple.

By using the rightmost path expansion strategy, it is possible to check the canonicity of a pattern by testing only the nodes belonging to the right path. Let l and p be two functions that return the last and penultimate child of a node, respectively. It is possible to determine the canonicity of a pattern in the following way: For all the nodes n belonging to the right path, if ∃ p(n) and code( p(n)) ≤ code(l(n)), then the pattern is in canonical form. It is, in the same way, fairly easy to identify the extensions that produce patterns with automorphisms. If at least one of the nodes belonging to the right path has more than one child and for all the nodes n belonging to the right path, if ∃ p(n) and code( p(n)) = code(l(n)), then the pattern has automorphisms. On such a pattern, we call split node the node from the right path that is closest to the root and that has several children with the same label. Thus, the patterns shown in Fig. 1c, d have automorphisms. On these two figures, the split node is the node a.

Handling reentrant links

From the enumeration of spanning trees described above, it is possible to define a method for enumerating all rooted subgraphs (since a rooted subgraph is a tree with reentrant links). For this, we need to identify when a structural extension adds a node that is already present in the candidate pattern. When this happens, the node is added to the pattern as if it is a normal structural extension, but we indicate that it is a reuse of a node already present by storing the index of the pointed node. The exploration of the search space is stopped from this node. Indeed, because using this exploration strategy, we are sure that all the nodes belonging to the right path are extended in the current generation phase. All structures that are not part of the rightmost path have already been fully explored. The reuse of a node is only reflected in the code of the graph by adding a triple that refers to another triple. This operation generates no change in the properties of the code, the completeness of the method and the way to eliminate noncanonical patterns.

Handling cycles

The extension of a pattern can lead to the creation of a cycle if and only if the new added node points to one of the nodes already present in the rightmost path. When a cycle is created, we proceed in the same manner as other reentrant links: A new node is added with the index of the reused node, and the extensions from this node are stopped. The completeness of the method is not affected by this operation. However, if the reused node is the root, we can generate redundant solutions. In such case, many isomorphic patterns are generated by starting the exploration from the other nodes belonging to the cycle.

We solve this problem by only extending one of these patterns. First, we generate all the isomorphic patterns s.t. their root is one of the nodes belonging to the right path. Then, we keep the pattern with the smallest canonical form and continue the exploration only for this pattern. In Fig. 3, P 1 represents a pattern generated from the graph G. The node labeled with d is the last node added to the pattern. This extension creates a cycle composed of 3 nodes belonging to the right path (d, c and bc). Thus, we have 2 other isomorphic patterns (one for each other node of the cycle). From the node bc, we obtain the pattern P 2 . From the node c, we obtain the pattern P 3 . The pattern that has the smallest code is the one rooted in bc, so we stop here the exploration of the pattern P 1 .

Handling automorphic patterns

The use of canonical forms avoids the generation of redundant solutions. It is used in our algorithm to prune the search space. However, this does not solve the problem posed by the presence of automorphic patterns. For example, let us consider the directed attributed graph shown in Fig. 4a. This graph contains several automorphic substructures that make less efficient the pruning based on canonical forms. When the pattern shown in Fig. 4b is analyzed, it is accepted and extended since it is in canonical form. However, this pattern can be obtained in twelve different ways. The naive algorithm consists, for each of these forms, to extend the pattern by using, as defined above, the itemset extension or the structural extension.

Itemset extension generates the pattern shown in Fig. 4c. This pattern will be rejected because it is not in canonical form. Structural extensions from node b generate one of the patterns shown in Fig. 4e. In all cases, the resulting pattern is not in canonical form and, as a consequence, it will be discarded. The only operation that can generate a canonical pattern is a structural extension from the split node a (Fig. 4d).

When we identify a new pattern with automorphisms, we know that this pattern cannot be extended on the last node or on all the nodes located after the split node. In fact, for a split node n, we have code( p(n)) = code(l(n)), then any extension done on the pattern of root l(n) will lead to a pattern with a code lower than p(n), so this pattern will not be in canonical form. We can use this fact to avoid generating unnecessary patterns. This is a first optimization, but its impact is limited.

The second optimization is to remove the mappings of a pattern that are not able to generate new canonical patterns. In Fig. 1, for example, the pattern of Fig. 4d, which is obtained in twelve different ways, can be only extended by adding a third child to the split node annotated a. Whatever way the pattern of Fig. 4d has been obtained, the possibilities of expansion are unchanged. It is therefore possible, in this case, without omitting solution, to keep only one mapping of the pattern.

With the graph in Fig. 4a, the situation is more complicated because it is possible, from this pattern, to add a child c to the node labeled a. In this case, it will be possible in the next generation, to apply any type of extension for all the nodes belonging to the right path. In the example, we could, in particular, add e, f , g or h as a child of node c. In order that all possible patterns could be generated, we must keep the possibility of extending the pattern with 2 bc nodes from the 3 bc nodes with a child. This represents 3 2 different mappings. In a real graph, the situation is often less simpler and determining what mappings it is possible to remove is not an easy task. To avoid performing calculations on the graph structure that would affect the performance of the algorithm, we choose to keep only one mapping among those using the same nodes. This decreases the number of mappings from a number of permutations to a number of combinations. This solution is not optimal, but it is sufficient to limit the combinatorial explosion.

In a star graph in which a node is associated with n nodes with the same label, there exists has Automor phisms, implementing the method described in Sect. 5.1). Candidates with automorphisms are processed in line 8 to 9. In line 8, the mappings of the patterns that are not able to generate new canonical patterns are removed (see Sect. 5.4). In line 9, extension of the patterns (stored in the set X ) are only performed for nodes located before the split node (function extend Be f ore that takes as parameters the candidate node G and the splitNode determined by the function get Split N ode. see Sect. 5.4). In line 11, candidates that do not have automorphisms are extended by the regular method extend.

In lines 13 to 15, the set of extensions is added to the set of candidates only if the current candidate G is not contained in a pattern belonging to the set of solutions with exactly the same occurrences. We use for the test, the method om Frequency that computes the total number of occurrences of a pattern in B (the occurrence-match frequency).

In line 16 to 18, the candidate is added to the set of solutions only if it is content closed; i.e., G is not contained in a pattern belonging to the set of solutions with the same frequency.

The processing of a candidate finishes by removing it from the candidates' list (line 21).

Experimental results

The method described in this paper has been implemented in C++ using STL. Experiments were performed on a computer running Ubuntu 13.04 and based on a Intel©Core™i5-2400 @ 3.10 GHz with 12 GB main memory. The hard disk used for the experiment is a 500 GB capacity with a SATA 3 Gb/s interface, 32 MB cache and 5400RPM spin speed. All timings are based on total execution time, including all preprocessing and results output.

Synthetic datasets

A set of 10,000 directed graphs, each comprising 1,000 nodes and 5000 edges, were generated with the program of [START_REF] Johnsonbaugh | A graph generation software package[END_REF]. This set of graphs was used as the basis for the creation of five datasets of directed attributed graphs named A1 to A5 in which nodes were associated with itemsets of size ranging from 1 to 5. To simulate the fact that, as in many real datasets, most attributes are rare, but a few are very common, we assigned each item a random integer value distributed according to a power law (more specifically, we have used a Pareto distribution with parameter α = 0.05).

Google+ and Twitter datasets

We used the social graphs of Google+ and Twitter built by [START_REF] Mcauley | Learning to discover social circles in ego networks[END_REF]. In each dataset, a node represents an individual and the associated attributes are the characteristics of this person.

There is an edge between two persons when one of them is following the activity of the other. The datasets are substantial (107,614 nodes and 13,673,453 edges for Google+, 81,306 nodes and 1,768,149 edges for Twitter).

PubMed Central citation network

PubMed Central is a bibliographic database of scientific articles in the field of life sciences. This database contains 322,526 Open Source papers but only 58,728 of them cite another paper which is itself Open Source. We built a citation graph where each node represents a paper, each edge, a citation link, and the attributes associated with nodes identify the papers' keywords. The graph is very sparse as it contains only 60,012 edges; however, some attributes are frequently associated with subsets of connected nodes.

Performances evaluation

The results of our algorithm are shown in Fig. 6. The comparison with existing algorithms is only possible with labeled graphs. The only dataset in which each node is associated with a single item is dataset A1 presented in Fig. 6a. For this dataset, the performance obtained with the implementations of Gspan carried out in the project ParSeMiS [START_REF] Wörlein | A quantitative comparison of the subgraph miners mofa, gspan, ffsm, and gaston[END_REF] (an implementation that allows the mining of directed graphs) is also presented. Gspan and AADAGE, which are set to mine directed, rooted and connected subgraphs, generate exactly the same results. For the smallest values of support, our strategy of mapping filtering allows making a difference with Gspan. Figure 6b shows the impact of itemset size (for the same graph) on algorithm performance. The execution times increase exponentially when the number of itemsets in each node increases. This highlights the difficulty of mining attributed graphs in comparison with labeled graphs. Indeed, when mining attributed graphs, we have to deal with the combinatorics of graphs and itemsets at the same time.

The Google+ and Twitter datasets contain many nodes, many attributes and especially some very common attributes (e.g., the attribute "gender = male," which is present in 52 % of the nodes) whose presence generates so many patterns that the mining fails. In Fig. 6c, the two gender attributes were removed from Google+ dataset. Twitter dataset has not been changed. Despite its size, Google+ dataset was successfully mined by setting a minimum support of 3 %. Twitter dataset, although smaller, has been treated up to a support of 5 %. The common feature of these two datasets is that they are treated fairly rapidly up to a With PubMed Central dataset, it was possible to use a very low support (Fig. 6d). This allowed identifying some interesting patterns. A pattern, for example, concerns a set of papers annotated with " p53" which is the name of an oncogene. These papers are cited by other papers annotated "rapamycin" (an immunosuppressant drug designed in 1975) that are themselves cited by papers annotated "aging" (the influences of rapamycin on the deceleration of cellular senescence were discovered in 2009) and "mT O R" (a gene that is inhibited by rapamycin). This pattern represents a highly summarized picture of the progress of some research on aging.

Conclusion and perspectives

We have presented in this paper a method to mine directed attributed graphs and shown its effectiveness in the analysis of multiple datasets. We focused specifically on graphs in which, as in real data, some attributes are shared by a large number of nodes. We have shown that, depending on the structure of the pattern being analyzed, it is not always necessary to extend all possible mappings of a pattern. We treated specifically the case of automorphic patterns, but there exists other patterns for which some optimizations can be performed. For example, when the pattern in Fig. 1b, which does not have automorphisms, is found 12 times on the graph of Fig. 4a, we can see quite easily that some mappings can be discarded. We are convinced that the detailed analysis of the patterns and their mappings may reveal other configurations where a reduction of the search space can be performed, which could make possible the mining of larger or denser datasets. It is a future work that we aim to explore.

Fig. 1

 1 Fig. 1 Example of a directed labeled graph (a) and 3 patterns (a, b, c and d) present in the graph (the numbers displayed on the edges show their creation order)

Fig. 2

 2 Fig. 2 Example of four isomorphic patterns (the numbers displayed on the edges show their creation order)

Fig. 3 Fig. 4

 34 Fig. 3 Illustration of redundant patterns related to cycle

Fig. 5

 5 Fig. 5 AADAGE Algorithm

Fig. 6

 6 Fig. 6 Performances of the algorithm on several datasets

Graphs with labeled edges can always be transformed into graphs with only labels on nodes (using, e.g., the method proposed by[START_REF] Inokuchi | Complete mining of frequent patterns from graphs: mining graph data[END_REF]). For this reason, we only consider attributed nodes in our study.

n i=1n P i different ways of generating all the included patterns. After the optimization, this number drops to n-1 i=0 (ni) n i .6 Mining algorithmMining all frequent rooted attributed subgraphs generates a huge number of patterns. In real applications, generating all solutions can be very expensive or even impossible. In addition, many of these frequent attributed subgraphs contain redundant information. On the other hand, mining frequent closed attributed graphs is very expensive in terms of processing time.Based on the experiments we conducted on attributed tree mining[START_REF] Pasquier | Frequent pattern mining in attributed trees: algorithms and applications[END_REF], we developed an algorithm called AADAGE (Automorphism Aware Directed Attributed Graph Explorer) that is designed to mine content closed attributed graphs. It has been shown (Pasquier et al.[START_REF] Pasquier | Frequent pattern mining in attributed trees: algorithms and applications[END_REF]) that mining content closed patterns is a good compromise between non redundancy of solutions and execution time.Figure5shows the high-level structure of the AADAGE algorithm. The algorithm takes in parameter a database (B) and a minimum support (minSup). First (line 1), a set C containing all attributed subgraphs of size 1 is built by scanning B. A subgraph of size 1 is composed of one node associated with one item. Therefore, C is the set of nodes associated with every item presents in B. In line 2, the set of solutions S is initialized with the empty set. Between line 3 and line 22, a loop allows to process every candidate in the set.The function Get First (line 4) returns the smallest candidate in the set according to the order of attributed graphs defined in Sect. 4.3. In line 5, the algorithm performs a canonical test (function isCanonical, implementing the method described in Sect. 5.1) and a frequency test (function f requency) and discards, from further processing, the candidates that do not pass these tests. In line 6, candidates for which the previous extension step leads to the creation of a cycle (detected with the function hasC ycle) are checked for canonicity with the method described in Sect. 5.3. Candidates that have a cycle and that are not canonical are discarded. In line 7, the algorithm performs a test for the presence of automorphisms (function

Acknowledgments This work was supported by the ANR Grant "FOSTER" ANR-2010-COSI-012-01.

Frédéric Flouvat is an Associate Professor at the University of New Caledonia (Nouméa, New Caledonia), where he is teaching Algorithmic and Databases. He is also a member of a multidisciplinary research team on material and environment. This laboratory brings together geologists, physicists and computer scientists to address both fundamental and applied questions related to the concepts of risk and sustainable development. His research interests are spatio-temporal data mining and its application to environmental sciences.

Jérémy Sanhes is a Ph.D. student at the University of New Caledonia and INSA Lyon (France). His thesis deals with spatio-temporal data mining issues.