Mapping Time Series to Networks Using Visibility Graph

A time series can be converted into a network by mapping the connection between timepoints using visibility criteria (Lacasa et al., 2008; Zhang et al., 2022). A visibility graph (VG) for a time series xi is defined as a graph G = (V, E) such that each time stamp, t, is a node (V) in the graph and the edge (E) between nodes ??, and ??, is a line of visibility between the signal amplitudes xi and xj.

Any two data points (ti, xi) and (tj, xj) will be connected nodes of the graph if any other data points (tk, xk) which lie between them meet the following criteria.

$$_< _+ _-_) \frac_-_) }_-_) }$$

This natural visibility can also be illustrated by using vertical bars to represent time-series data (Fig. 1b). These VGs are connected graphs, in which each node has at least one connection (i.e., neighbouring timepoints are connected to each other). This network is invariant under affine transformations and vector translations of the data. Weighted VGs (Silva et al., 2022) can be obtained by using Euclidian distance between timepoints as weighting factor such that:

$$_= \frac_- _) }^+_- _) }^ }}$$

VG graph of fMRI time-series represent a simplified mathematical construct to generate graph format representation of timeseries. Each node in the graph represents a time-point and the binary edges represent visibility of signal at a timepoint to another. Thus, edges represent temporal relationship in the data. VG networks are useful for generating unique features, which can reveal certain patterns or relationships within the time series data. However, it cannot provide insights into the specific neural connections, synapses, or the underlying spatial architecture of the brain (Silva et al., 2022).

An important aspect of VG analysis is that there is a natural correspondence of nodes (i.e., timepoints in VG) across different brain regions and subjects. Hence, any local features of nodes such as degree-sequence can be compared across various brain regions. One previous work has used this natural alignment within multiplex visibility framework as compact way of extracting at once both the local temporal structure and the global connectivity pattern (Sannino et al., 2017). Our current work aimed to go back one step and investigate whether basic local and global properties of visibility graphs are reliable and robust.

Graph Theoretical Features of Visibility Graphs

The variability of time-series data can be analysed using graph-theoretical methods. These analyses can extract features such as centrality, distance, community structure, and connectivity, which are important for understanding the characteristics of a graph. In this work, we use five global properties of the graph - average weighted degree, average path length, global clustering coefficient, number of communities, and modularity - to characterise the time-series graphs. A brief description of these measures is provided below.

The average weighted degree is the arithmetic mean of the weighted degrees of all nodes in a network. The weighted degree is a measure of the strength of connectivity for each node and characterises the intensity of connectivity in the node’s neighbourhood. The average path length is the mean of the shortest paths between all pairs of nodes. This is a measure of the flow of information in the network. The global clustering coefficient is a measure of the extent to which the nodes of a graph tend to cluster. It measures the probability that two nodes connected to a given node are also connected. The number of communities is the number of clustered groups of nodes in the network. The modularity, Q, measures how well the graph can be divided into communities. A high modularity indicates a graph with a dense internal community structure and sparse connections between nodes of different communities.

The Walktrap Community-finding algorithm was utilized for community detection. This algorithm employs a random walks approach to identify regions of the network where nodes are more likely to be interconnected, indicating the presence of a community. Furthermore, to calculate the modularity of a graph, we determined the degree of separation of nodes belonging to different communities using the approach described in a previous work (Silva et al., 2022). We used the NetF toolbox to generate VG features from fMRI timeseries (Silva et al., 2022). The toolbox is available via github (https://github.com/vanessa-silva/NetF).

Network Connectivity Using Degree Synchrony

The VG networks can also be represented as a system of multilayer network in which each brain region forms a layer (Fig. 1) (Ahmadlou & Adeli, 2012). Such a multilayer network has one-to-one correspondence between nodes of each layer such that node i in one layer corresponds to the same node in other layers, making the multilayer graph a multiplex network. Graph theoretical features of multiplex graph can be used to measure similarity between layers.

The degree synchronization is a measure of similarity between the series of connectivity degree of each layer (Ahmadlou & Adeli, 2012). Since each node in each layer represent a time-point, the series of connectivity degree of visibility graphs at layers is a representation of fluctuations in degree over time. Thus, any correlations in fluctuations in visibility graph degrees over time across multiple layers represents an alternative measure of synchronization. Thus, degree synchronization is given as:

$$S \left[Dx, Dy\right]= \frac) \left(Dy-\overline \right)}$$

where Dx and Dy are degree sequence of visibility graphs representing timeseries x and y.

While correlation-based connectivity mapping of fMRI time-series provides valuable static insights into functional connectivity, degree synchrony-based connectivity mapping introduces a unique dimension, emphasizing how fluctuations in the amplitude of time signals synchronize across different regions over time. Although more work is needed to better understand the biological underpinnings of degree-synchrony-based connectivity mapping, it represents an important step towards a computationally efficient method with the potential to capture multi-scale dynamics in fMRI data.

Fig. 1

A framework for visibility graph analysis of fMRI data. a fMRI data is parcellated into different brain regions using an atlas b Average time series is extracted from all brain regions. This time-series data can be represented as a temporal landscape such that visibility between each data point can be identified. Visibility between two time points exists (i.e., an edge) if any other time point between them has a corresponding intensity below the line connecting the two time points (i.e., there is a direct line-of-sight between the peaks of time points) c The graphs generated using the visibility criteria has number of nodes equal to number of time points in the data. The graphs can then be processed using standard graph theoretical analysis methods to generate graph features. d The visibility graph degree vectors from each ROI can be correlated to generate a Degree Connectivity network, providing a new measure of functional connectivity

Application of Visibility Graph Analysis in rs-fMRI DataParticipants and Data

The dataset used in the present study was obtained from the Human Connectome Project (HCP) S1200 release (Van Essen et al., 2013). Participants (n = 1113) were young healthy adults aged 22 to 37 years. Each participant took part in two sessions (conducted on two consecutive days) of resting-state fMRI scans acquired over two runs (right-to-left and left-to-right phase encoding) of 14m33s each. The MRI data were acquired on a customized Siemens Skyra 3 Tesla MR scanner using a multiband echo planar imaging sequence (TR = 720ms, TE = 33.1ms, voxel dimension: 2 × 2 × 2 mm3) (Smith et al., 2013). The current study uses data from the participants who completed all four runs with a final sample of 1010 individuals (mean age 29 ± 4 years, 453 males).

fMRI Pre-processing

We used fMRI data processed using HCP pipeline and available publicly. The detailed pre-processing steps used in HCP pipeline are described elsewhere (Glasser et al., 2016; Smith et al., 2013). Briefly, the fMRI pre-processing steps within the HCP data included (1) removal of spatial artifacts and distortions, (2) correction of head motion, (3) spatial registration to the MNI (Montreal Neurological Institute) standard space and (4) removal of motion-related and structured physiological noise artefacts using ICA-FIX (Salimi-Khorshidi et al., 2014). Data were analysed in CIFTI (Connectivity Informatics Technology Initiative) format, in which cortical surface time series and subcortical volume time series were structured in a grayordinate dense time series file. Head motion was quantified using framewise displacement (Power et al., 2015).

Estimation of Visibility Graph Features

The pre-processed fMRI time-series data from each session (REST1 and REST2) were temporally concatenated across the two runs within the session. This generated approximately 30 min of rs-fMRI data per session. The fMRI data were parcellated using Yeo-17 functional atlas (Yeo et al., 2011) for obtaining average time series across functionally-defined brain regions (114 brain regions). The time-series data obtained from the parcellation were then converted into VGs using natural visibility graph algorithms in R (Silva et al., 2022). To characterise the degree properties of the VG obtained from fMRI time-series, log-log degree distribution was estimated for each subject. The degree distributions were obtained from several regions of interests within the Yeo-17 atlas for demonstration purpose. The degree distribution plots were pooled across participants and plotted. Power-law fit was obtained by using powerLaw package in R. The VG for each region was analysed using graph theoretical analysis tools in implemented in R by a previous work (Silva et al., 2022) to obtain the five features: average weighted degrees, average path length, global clustering coefficient, number of communities and modularity. In order to characterise degree distribution of VG associated with fMRI time-series, log-log plot of degree distribution were generated for a region from each canonical functional network. Power law fit was approximated for each degree distribution plot using PoweRlaw package in R.

Impact of Motion and Test-Retest Reliability

To identify the impact of motion intrusions on VG features and identify the threshold required for removing any motion related impacts, we characterised the relationship between proportions of motion corrupted data and VG features across the participants. First, to estimate the impact of motion on the features, we calculated the correlation (Pearson’s) between each feature and the percentage of fMRI data points associated with motion across the participants. This percentage was defined as the percentage of fMRI data points associated with frame-wise displacement (FD) greater than 0.2 mm (i.e., data points with greater than 0.2 mm FD divided by total data points). Next, to determine the allowable percentage of data points with motion without compromising the VG features, we ran the correlation analysis using various subsets of the data, with each subset exclusively comprising subjects exhibiting less than x% of frames characterized by frame displacement (FD) exceeding 0.2 mm. The value of x ranged from 10 to 40% in the steps of 0.5%. This allowed us to identify appropriate motion thresholds (in terms of tolerable percentage of motion corrupted data points) necessary for reliable VG features.

To estimate test-retest reliability of VG features between two sessions acquired on different days (REST1 and REST2), we computed the intraclass correlation coefficient (ICC). A two-way mixed effects model with absolute agreement as a reliability measure was used, as per the recommendation from previous work (Koo & Li, 2016). This analysis was applied only to the low motion dataset identified using the motion threshold established in the previous section.

Functional Network Estimation Using Degree Synchrony

Mapping of functional network using VG degree synchrony is a novel approach. Only one previous study has attempted to do this in a small dataset (Gao et al., 2020). We used the low motion fMRI dataset and mapped the functional networks in the brain using the degree synchrony. We generated degree synchrony based functional connectivity maps for all individuals, which were averaged to obtain a mean connectivity map.

The primary motivation for this analysis was to introduce methodological approach for generating these networks, paving the way for future in-depth analyses. This analysis should be taken as an initial step in hypothesis generation, opening the door for subsequent detailed analyses. The detailed analysis of network properties and the impact of motion on these analyses is outside the scope of this manuscript.

To capture the degree synchrony, a pairwise correlation analysis was performed on the VG degree time series of each region, yielding a matrix of inter-regional correlations. The construction of the functional network was achieved by thresholding the degree synchrony matrix to retain connections with significant synchrony.