Characterization of Cortical Networks and Corticocortical Functional Connectivity Mediating Arbitrary Visuomotor Mapping

Experimental conditions and behavioral tasks.

Eleven healthy participants accepted to take part in our study (all were right handed; the average age was ∼23 years; 4 were females and 7 males), gave written informed consent according to established institutional guidelines and local ethics committee, and received monetary compensation (€50). We used an associative visuomotor mapping task, where the relation between visual stimulus and motor response is arbitrary and deterministic (Wise and Murray, 2000). The domain of visuomotor control distinguishes two forms of guidance: standard and nonstandard (or arbitrary) mapping. In standard mapping, such as reaching toward an object, the spatial properties of the visual cue (i.e., target location) provide relevant information for visuomotor transformation. In arbitrary mapping, the visual cue and the action are not spatially congruent; and other properties of the object, such as its form or color, are transformed to plan actions (Petrides, 1987; Wise et al., 1996; Murray et al., 2000; Passingham et al., 2000; Wise and Murray, 2000). To investigate arbitrary visuomotor mapping, we asked participants to perform a finger movement associated to a digit number: digit 1 instructed the execution of the thumb, digit 2 for the index finger, digit 3 for the middle finger, and so on (Fig. 1). Maximal reaction time was 1 s. After a fixed delay of 1 s following the disappearance of the digit number, an outcome image was presented for 1 s and informed the subject whether the response was correct, incorrect, or too late (if the reaction time exceeded 1 s). Incorrect and late trials were excluded from the analysis because they were either absent or very rare (i.e., maximum 2 late trials per session). The next trial started after a variable delay ranging from 2 to 3 s (randomly drawn from a uniform distribution) with the presentation of another visual stimulus (Fig. 1b). Each participant performed two sessions of 60 trials each (total of 120 trials). Each session included three digits randomly presented in blocks of three trials.

Anatomical, functional, and behavioral data acquisition.

Anatomical MRI images were acquired for each participant using a 3-T whole-body imager equipped with a circular polarized head coil. High-resolution structural T1-weighted anatomical image (inversion-recovery sequence, 1 × 0.75 × 1.22 mm) parallel to the anterior commissure-posterior commissure plane, covering the whole brain, were acquired. MEG recordings were performed using a 248 magnetometers system (4D Neuroimaging, magnes 3600). Visual stimuli were projected using a video projection, and motor responses were acquired using a LUMItouch optical response keypad with five keys. Presentation software was used for stimulus delivery and experimental control during MEG acquisition. Reaction times were computed as the time difference between stimulus onset and motor response. Sampling rate was 2034.5 Hz. Location of the participant's head with respect to the MEG sensors was recorded both at the beginning and end of each session to potentially exclude sessions and/or participants with large head movements. However, none of the participants moved >3 mm during all sessions. Thus, all participants were considered for further analysis.

Preprocessing and spectral analysis of MEG signals.

MEG signals were down-sampled to 1 kHz, low-pass filtered to 250 Hz, and segmented into epochs aligned on finger movement (i.e., button press). Epoch segmentation was also performed on stimulus onset and the data from −0.5 and −0.1 s before stimulus presentation was taken as baseline activity for the calculation of the single-trial HGA. Artifact rejection was performed semiautomatically. For each movement-aligned epoch and channel, the MEG signal variance and z-value were computed over time and taken as relevant metrics for the identification of artifact epochs. All trials with a variance >1.5 × 10⁻²⁴ across channels were excluded from further analyses. Additional metrics, such as the z-score, absolute z-score, and range between the minimum and maximum values, were also inspected to detect artifact. Two MEG sensors were excluded from the analysis for all subjects.

Spectral density estimation was performed using a multitaper method based on discrete prolate spheroidal (slepian) sequences (Percival and Walden, 1993; Mitra and Pesaran, 1999). To define the frequency bandwidth of interest for the estimate of single-trial gamma-band activity, we performed time-frequency analyses of the MEG time series at the sensor level. MEG fields were transformed to a planar gradient configuration (i.e., by computing the gradient tangential to the scalp). Time-frequency power estimates were then computed for all sensors over the gamma band from 30 to 150 Hz (in steps of 5 Hz) using 4 orthogonal tapers 0.25 s in duration and 20 Hz of frequency resolution, each stepped every 0.01 s. Time-frequency power estimates were log-transformed and then normalized (i.e., z-transformed) with respect to the mean and SD computed over a baseline period from −0.5 to −0.1 s before stimulus onset. Time-frequency maps were averaged across virtual gradiometers and participants. Results showed increases in power with respect to baseline in the high-gamma range from ∼60 to 120 Hz (Fig. 2). Based on the spectral analysis at the sensors level, subsequent analyses focused on the high-gamma band from 60 to 120 Hz. MEG time series were multiplied by k orthogonal tapers (k = 8) (0.15 s in duration and 60 Hz of frequency resolution, each stepped every 0.005 s), centered at 90 Hz, and Fourier-transformed. Complex-valued estimates of spectral measures X_sensorⁿ (t, k), including cross-spectral density matrices, were computed at the sensor level for each trial n, time t, and taper k.

Source analysis and calculation of HGA.

Source analysis requires a physical forward model or leadfield, which describes the electromagnetic relation between sources and MEG sensors. The leadfield combines the geometric relation of sources (dipoles) and sensors with a model of the conductive medium (i.e., the headmodel). For each participant, we generated a headmodel using a single-shell model constructed from the segmentation of the cortical tissue obtained from individual MRI scans (Nolte, 2003). Leadfields were not normalized. As source location, a 3D grid with regular spacing between the dipole locations of 10 mm was generated for each participant. Individual MRI scans were then warped to the template MRI in MNI space, and the normalization parameters were applied to the dipole grid. Such procedure assures that individual subjects' grid points are located in equivalent brain areas across all subjects according to MNI space. The headmodel, source locations, and the information about MEG sensor position were combined to derive single-participant leadfields.

We used adaptive linear spatial filtering (Veen et al., 1997) to estimate the power at the source level. We used the Dynamical Imaging of Coherent Sources method, a beam-forming algorithm for the tomographic mapping in the frequency domain (Gross et al., 2001), which is well suited for the study of neural oscillatory responses based on single-trial source estimates of band-limited MEG signals (for review, see Hansen et al., 2010). At each source location, Dynamical Imaging of Coherent Sources uses a spatial filter that passes activity from this location with unit gain while maximally suppressing any other activity. The spatial filters were computed on all trials for each time point and session and then applied to single-trial MEG data. Dynamical Imaging of Coherent Sources allows the estimate of complex-value spectral measures at the source level, X_sourceⁿ (t, k) = A(t) X_sensorⁿ (t, k), where A(t) is the spatial filter that transforms the data from the sensor to source level (for a detailed description of a similar approach, see Hipp et al., 2011). The single-trial high-gamma power at each source location was estimated by multiplying the complex spectral estimates with their complex conjugate, and averaged over tapers k, P_sourceⁿ (t) = X_sourceⁿ (t, k) X_sourceⁿ (t, k)_k^*. Single-trial power estimates aligned on movement and stimulus onset were log-transformed to make the data approximate Gaussian and low-pass filtered at 50 Hz to reduce noise. Single-trial mean power and SD in a time window from −0.5 and −0.1 s before stimulus onset was computed for each source and trial, and used to z-transform single-trial movement-locked power time courses. Similarly, single-trial stimulus-locked power time courses were log-transformed and z-scored with respect to baseline period, so to produce HGAs for the prestimulus period from −1.6 to −0.1 s with respect to stimulation for subsequent FC analysis. Finally, the anatomical position of each source was labeled according to BA using the binary representation of the Talairach–Tournoux atlas (Talairach and Tournoux, 1988) digitized for the Talairach Daemon (Lancaster et al., 2000). Single-trial HGA at each BA was defined as the mean z-transformed power values averaged across all sources within the same BA and cerebral hemisphere (a total of 76 BAs covering both hemispheres). The preprocessing steps, artifact rejection, spectral analyses, and source analysis were performed using FieldTrip toolbox (Oostenveld et al., 2011).

Linear correlation.

A conventional undirected connectivity measure is linear correlation. We computed the linear correlation coefficient between the HGAs of pairs of BAs at the single-trial level over a time window ranging from −1 to 0.5 s around movement onset, thus including 300 time samples (i.e., HGA were computed every 0.005 s). Single-trial correlation coefficients were also computed among pairs of HGA that were computed during the prestimulus period, from −1.6 to −0.1 s before stimulus onset. The statistical analyses then searched for significant differences between the movement-related correlations coefficients and those computed in the prestimulus interval.

Covariance-based Granger causality measures.

Granger causality (Granger, 1980; Brovelli et al., 2004; Ding et al., 2006; Bressler and Seth, 2011; Seth et al., 2015) is a directed FC measure. Granger (1963, 1980) proposed a statistical criterion to infer causality from process X to process Y based on the extra knowledge that can be obtained about the future of Y given the past of X, in a given context Z. This criterion can be formalized as the condition of independence as follows: where the superindex i refers to the whole past in a time window of length T of a process up to and including sample i, V refers to the whole system {X,Y, Z}, and Xⁱ indicates the exclusion of the past of X. That is, X is Granger noncausal to Y given Z if the above equality holds. This condition of independence constitutes a strong form of Granger noncausality (Granger, 1980), which can generally be tested within an information theoretical framework by calculating the Kullback–Leibler divergence (Cover and Thomas, 2006) of the two probabilities in Equation 1. The Kullback–Leibler divergence is a conditional mutual information and, for a bivariate context comprising only X and Y, it is defined as follows: where the conditional entropies are as follows: and The conditional entropies in Equations 3 and 4 quantify the uncertainty of Y_i+1 given the knowledge about the pasts Yⁱ only or given the past of both Xⁱ and Yⁱ, respectively. Accordingly, F_X→Y (Eq. 2) compares the uncertainty in Y_i+1 when using knowledge of only its own past or both Xⁱ and Yⁱ. Such conditional mutual information has been studied in different fields, such as communication theory and econometrics, and it has been formulated under different assumptions (e.g., stationarity vs nonstationarity) and terminology (for details, see Chicharro and Ledberg, 2012). In neuroscience, it is best known as the Transfer entropy (Schreiber, 2000). Although Transfer entropy in its general form (Eq. 2) does not require a model of the interaction and it is powerful in capturing any (linear and nonlinear) conditional dependence between Y_i+1 and Xⁱ (Vicente et al., 2011; Wibral et al., 2014), its calculation is not trivial when high-dimensional probability distributions need to be estimated over short time windows as often encountered in neurophysiological data (for a review of the estimation of information theoretic measures to calculate Transfer entropy, see Hlaváčková-Schindler et al., 2007). To estimate FC measures on short data segments, such as single trials, we implemented an approximation of Transfer entropy that exclusively considers the first term of its expansion (also referred to as the Gram-Charlier expansion as used in Amari et al., 1995), which quantifies the contribution of second-order statistics (i.e., covariance matrix). In other words, our approach is equivalent to implementing a weaker criterion of Granger noncausality, called noncausality in mean (Granger, 1980). This criterion compares the errors of the optimal linear predictors of Y_i+1, minimizing mean squared error, which corresponds to the conditional variance of the future given the past, so that Equation 1 is reduced to the following: where σ²(·|·) denotes the conditional variance. The statistical measure used to test Equation 5 is commonly referred to as the Granger causality. In particular, the classical parametric implementation (Granger, 1963) uses autoregressive models to estimate prediction errors. Given that prediction errors can also be estimated from the covariance matrix (Lütkepohl, 2005; Chicharro, 2011), the conditional variances in Equation 5 determine the first term in the expansion of the entropies of Equation 2. This is because, as previously mentioned, the first term of the expansion quantifies the contribution of second-order statistics to the entropy. Put it into practice, given an interval T defining the portion of data Xⁱ and Yⁱ used to compute covariance measures and a lag determining the number of time points of the pasts Xⁱ and Yⁱ used, submatrices of the covariance matrix of the variables comprising the past and the future, {X_i+1,X_i,X_i-1,…X_i-lag+1, Y_i+1,Y_i,Y_i-1,…Y_i-lag+1}, are used to estimate the entropies (Cover and Thomas, 2006). For example, where Σ denotes a covariance matrix and |·| denotes the determinant. Within this approach, stationarity is assumed for the calculation of covariance matrices, whereas Gaussianity and linearity are not, and only determine how close is the in mean criterion (Eq. 5) to the strong criterion formalized in Equation 1. That is, for Gaussian variables, Equation 6 constitutes an exact expression for the entropy because higher-order moments are determined by the second-order moments, and Transfer entropy can be considered equivalent to Granger causality (Barnett et al., 2009). Whenever such equivalence does not hold (e.g., for non-Gaussian and nonlinear systems), the covariance-based approximation of the Transfer entropy is always equivalent to the parametric estimation of Granger causality. In addition, it has the advantage of not requiring fitting autoregressive models as an intermediate step. Given this equivalence, we will refer to our estimator as covariance-based Granger causality.

The relation between information-theoretic measures of dependence and Granger causality measures is not restricted to the case of Transfer entropy and linear Granger causality. Measures of Granger causality appear as terms in the decomposition of total interdependence between two time series. It has been formulated in information theoretic (Marko, 1973; Rissanen and Wax, 1987) and linear system (Geweke, 1982, 1984; Chicharro, 2011) frameworks, and applied to neural signals in the spectral domain (Ding et al., 2006). In the time domain and for finite time series, this decomposition is expressed as follows: and it represents the total interdependence between X and Y (Geweke, 1982; Chicharro and Ledberg, 2012). F_X,Y quantifies the dynamic increase of the total interdependence between two time-series at a given point in time, in contrast to the static interdependence quantified by linear correlation, and for an infinite lag converges to the mutual information rate. Such total interdependence is the sum of three Granger causality measures: two directed measures plus the “instantaneous” Granger causality term, accounting for unconsidered common influences to the processes. The instantaneous Granger causality is defined as follows: and can be expressed in terms of conditional entropies, The total interdependence is expressed as follows: The four measures in the decomposition of Equation 7 can be estimated with our covariance-based approach, implying the assumptions and approximations discussed above. Covariance-based Granger causality measures (F_X→Y, F_Y→X, F_X·Y and the total interdependence F_X,Y) were computed between pairs of HGAs on a single-trial basis over a time window T ranging from −1 to 0.5 s around movement onset and for the baseline period, from −1.6 to −0.1 s before stimulus onset, similarly to linear correlation analysis, giving a total of 300 samples in each window. The lag used for the calculation of the covariance matrices was 30 time points (i.e., 15% of T). The MATLAB (The MathWorks) code developed for the calculation of the covariance-based Granger causality measures can be downloaded from the website of the corresponding author.

Relation between linear correlation and covariance-based Granger causality measures.

Linear correlation and covariance-based Granger causality measures share common properties. As previously mentioned, linear correlation and total interdependence are symmetric measures quantifying static and dynamic dependencies, respectively. Although these measures are not related by a decomposition analogous to Equation 7, there is a strong relationship between the existence of both types of dependencies. A lack of total interdependence implies a lack of linear correlation; and, if we assume that the future of X and Y causally depends on their own past, respectively, the opposite relation is also true. This occurs because linear correlation ρ is related to the covariance-based approximation of the mutual information I(X_i+1;Y_i+1), which is equal to the logarithm of 1 − ρ², and because conditioning on the past cannot create new dependencies (Chicharro and Panzeri, 2014). It is also clear from Equation 7 that the directed and instantaneous Granger measures are smaller than the total interdependence. Thus, null total interdependence implies the absence of Granger causality measures because they constitute non-negative contributions to the total interdependence. In other words, Granger causality is present if, and only if, both linear correlation and total Granger interdependences are not zero.

Neural mass model (NMM) networks.

Neural signals were simulated using convolution-based neural mass models of local field potentials. Convolution-based NMM consider cortical mesocolumns, rather than a single cell's electrophysiological properties, and they are particularly suited for modeling neurophysiological responses, such as those recorded using the EEG and MEG (David et al., 2006; Moran et al., 2013). The model includes three cell subpopulations: spiny stellate cells in granular layer IV, pyramidal cells, and inhibitory interneurons in extragranular layers (II and III, V and VI). For each node, 13 measured variables include currents and membrane potentials of these three cell subpopulations. A connection between columns A and B corresponds to a link between the pyramidal cell of column A and the spiny stellate cells of column B. Exogenous input can also be included and modeled as an excitation of spiny stellate cells. The generated neural population signals from each node of the network integrate the membrane potentials from three cell subpopulations. We generated NMM signals with two networks depicted in Figure 3a. The first network was composed of 5 nodes with unidirectional links drawn in green. The second network was composed of all 10 nodes with unidirectional links shown in red and green (Fig. 3a). The goal was to simulate two functional networks associated to baseline (the first network) and event-related (the second) activities, respectively. Using such network functional connectivities, we generated neural signals that were subsequently cut into 50 epochs (trials), each of 150 time points. The number of time points was set to simulate the temporal intervals considered for subsequent HGA analyses. Indeed, the temporal window of interest for FC measures was set 1.5 s (from −1 s to 0.5 s with respect to motor response). HGAs were then expected to contain both ongoing and response-locked modulations, leading to potential latency shifts between areas. To specifically examine the influence of nonstationarities in the network across nodes, such as event-related modulations, we “injected” the networks with an exogenous input via node 1. The input was a Gaussian waveform centered at the 75th sample at each trial (the middle of the trial) and variance equal to 10 samples. The amplitude of the Gaussian input was varied from 0 to 0.04 in steps of 0.005. This injected Gaussian made the data nonstationary. To simulate different levels of stochastic background activity (i.e., noise level), we added Gaussian white noise and varied the amplitude from 0.01 to 0.1 in steps of 0.005. We generated datasets for each network type by varying parametrically the amplitude of the exogenous stimulus and noise level.

Statistical inference of FC using linear correlation and Granger causality measures.

For a given dataset with a given input and noise amplitude, we computed three measures of statistical interdependence on a single-trial basis. Given the scope of the numerical simulations, we considered undirected connectivity measures only: (1) linear correlation coefficient, (2) partial correlation coefficient, and (3) total Granger interdependence F_X,Y (see previous sections). These connectivity measures were computed from NMM signals using different window lengths (T ranged from 20 to 110 in steps of 2 samples). For the total Granger interdependence, the lag was varied from 2 to 0.2 × T. The maximal lag value was set to 0.2 × T to avoid small sample size problems in the estimation of covariance matrices, required for the calculation of conditional entropies. Total Granger interdependence was log-transformed to approximate Gaussian distribution.

The goal here was to recover from the simulated neural data the undirected FC (the presence of a link between two nodes, rather than the direction of influence) of the graph depicted in red in Figure 3a. This functional network corresponds to the difference between the second and first network. To do so, we performed a paired two-sample t test between the connectivity measures computed on the first and second network. More formally, we tested the null hypothesis H₀: X_ij(1) = X_ij(2), where X_ij(w) is a given connectivity measure between node i and j of network w. The absolute value of the log-transformed p value associated to the t test was tabulated in a sample connectivity matrix of n-by-n, where n is the number of nodes. The sample connectivity matrix was then compared with the known connectivity matrix corresponding to the red network in Figure 3a using a receiver operating characteristic (ROC) analysis and by computing the area under the ROC curve (AUC). The AUC was calculated for all three connectivity measures (linear correlation, partial correlations, and total Granger interdependence) with varying window lengths T (and lag for the total Granger interdependence) on all datasets generated with varying input amplitude and noise levels. To analyze the differences between free parameters and their associated variances, we performed a three-way ANOVA on the AUC computed on the statistical analyses of linear correlation coefficients. The three factors were stimulus amplitude, noise amplitude, and window length T, respectively. All factors showed statistically significant modulations (Table 1). However, the sum of squares of each source (Table 1) shows that the window length is the factor with the largest variance. This suggests that input amplitude and noise level have negligible effect on AUC compared with the effect of window length. We therefore averaged the AUC values across different input amplitude and noise levels and plotted the mean value with associated SEs as a function of window length (Fig. 3b, blue). The graph shows that the mean AUC rapidly increases as window length increases, and it reaches a plateau at ∼70 samples. In general, the AUC is acceptable (>0.95) for T >40 samples. The red curve displays the mean AUC values for the partial correlation coefficient. Here the curve rises more slowly with window length, and it never reaches stable value >95%. These results suggest that linear correlation analysis is more reliable for estimating FC and less sensitive to window length than partial correlation. Finally, we computed the mean AUC averaged over input amplitude and noise levels for the total Granger interdependence analysis. Figure 3c shows the mean AUC values for different values of window length T and lag. Similarly to the previous results, the AUC values increase with longer window lengths. A stable AUC (>0.95) is achieved for window lengths >40 and lags >4 (i.e., >10% of T). These results suggest that the total Granger interdependence can be used reliably for identifying changes in network FC, and it is comparable to linear correlation analysis. However, in addition to measuring the symmetrical dependences, the total Granger interdependence can be decomposed into the directed and instantaneous terms, which can provide directional information of network FC. Overall, the results of the simulations suggest that our method is a good candidate for characterizing cortical FC and may represent an appropriate tool for the analysis of neural data and HGAs.

Hierarchical approach for the characterization of network's nodes and FC.

We defined a functional cortical network as a set of BAs whose HGA and connectivity measures display significant task-related modulations. To assess statistical significance, we used an LME approach (see Statistical analysis). We used random-intercept and random-slope in LME models, where the “intercept” models the average between baseline and movement-related activity, and the “slope” models the difference between baseline and movement-related activity. Significant changes in movement-related HGA and connectivity measures were identified by testing whether the coefficient of the fixed-effect β₁ was significantly larger than zero. To account for multiple-comparison problems, significance level was set to q < 0.001, corrected for the FDR.

Network characterization was performed at four levels of analysis using a hierarchical approach. First, we identified the group of BAs displaying a significant increase in movement-related HGA with respect to mean baseline HGA (averaged from −0.5 to −0.1 s before stimulus onset). Second, among these BAs, we determined which pairs of BAs displayed a significant increase in linear correlation during the movement-related period with respect to baseline. Third, among those significant pairs, we characterized the pairs of BAs exhibiting a significant increase in total Granger interdependence F_X,Y. Fourth, among those latter pairs of BAs, we searched for significant increases in directed and instantaneous Granger causality. Such a hierarchical approach has several advantages in terms of exhaustiveness and interpretability of the results. It considers a functional cortical network as composed of links between BAs that show significant effects both in static (i.e., linear correlation) and dynamic dependencies (total interdependence and Granger causality). Such an approach may be especially suited for the interpretation of directed Granger causality measures. In general, even though Granger causality can be used to determine which causal structures are compatible with observed dependencies (e.g., Chicharro and Panzeri, 2014), it does not provide a measure of “true” causality (Pearl, 2009). Therefore, Granger causality is often combined with additional measures that can provide complementary information about the properties of the network under study. If spectral representations of Granger causality are used, the most relevant measures are the magnitude squared coherence among nodes and time-frequency power maps. These can reveal the presence of oscillatory activity (i.e., a peak in the average power spectrum) and the spatial extent (i.e., a peak in the average coherence spectrum). This approach has been used, for example, in Brovelli et al. (2004) where the combined application of power, coherence, and Granger causality measures provided a valuable tool for in-depth investigation of the nature of functional coupling of distributed neuronal assemblies. This approach excludes discriminating potential cases in which significant Granger causality is present among nodes without a clear oscillatory activity and/or coherence. Similarly, in the time domain, as in the current study, we interpret directional Granger causality measures only among nodes previously identified as part of task-dependent functional large-scale network of interest, where BAs display an increase in HGA power, linear correlation, and total Granger interdependence. Furthermore, we checked the robustness of our results and ran all-to-all connectivity analysis. Results showed that the core of the identified functional network is robust to the approach used (results not shown).

Detection of functional subnetworks (link communities).

A critical step in the analysis of networks is the detection of communities or groups of nodes potentially corresponding to functional subunits or modules. Communities, however, may overlap (i.e., a given node may belong to more than one group). Large-scale cortical networks may be no exception. For example, a single BA in associative cortical areas, such as BA7 in the superior parietal cortex, may contribute to processing visual inputs in an occipitoparietal network while participating in visuomotor transformations within a frontoparietal network. Link communities (i.e., groups of links rather than nodes) provide an appropriate framework capturing the relationships between overlapping communities while revealing hierarchical organization (Ahn et al., 2010). We adopted such framework for the analysis of adjacency matrices based on linear correlation and Granger causality networks (levels 2, 3, and 4) to identify functional subnetworks. Detection of link communities was based on the calculation of similarity measures between links from undirected and weighted networks. Single-linkage hierarchical clustering was used for the creation of link dendrograms. Link communities were found by cutting the dendrogram at the maximum partition density D, which measures the quality of a link partition. The partition density has a single global maximum along the dendrogram in almost all cases because its value is the average density at the top of the dendrogram (a single giant community with every link and node) and it is very small at the bottom of the dendrogram (most communities consist of a single link). Link community analysis was performed using a Python implementation downloaded from https://github.com/bagrow/linkcomm.