Dependence of connectivity on geometric distance in brain networks

In any network, the dependence of connectivity on physical distance between nodes is a direct consequence of trade-off mechanisms between costs of establishing and sustaining links, processing rates, propagation speed of signals between nodes. Despite its universality, there are still few studies addressing this issue. Here we apply a recently–developed method to infer links between nodes, and possibly subnetwork structures, to determine connectivity strength as a function of physical distance between nodes. The model system we investigate is brain activity reconstructed on the cortex out of magnetoencephalography recordings sampled on a set of healthy subjects in resting state. We found that the dependence of the time scale of observability of a link on its geometric length follows a power–law characterized by an exponent whose extent is inversely proportional to connectivity. Our method provides a new tool to highlight and investigate networks in neuroscience.


Results: time Scale of observability vs. Distance
The reconstruction of the cortical activity was carried out according to the atlas by Glasser et al. 20 , which provides the position of 360 different nodes. For the sake of computational simplicity, we analyzed 1/5 of the available set, namely 72 randomly selected nodes. For each single node a set of 60 time series is available, corresponding to 20 subjects and 3 recordings per subject. The set of nodes results in 2556 pairs. Due to the slight anatomical differences between the subjects, each pair corresponds to a set of 20 link lengths, so that the total number of d values is approximately 50000 within the range from 5 to 160 mm.
Given a subject, a recording and a pair of nodes, the time scale of observability of the corresponding link was assessed out of the related pair of time series, resulting in approximately 150000 values of W. While ∼ 70000 assessments failed to produce a finite time scale W, the remaining ∼ 83000 ones provided a valid W value within the range from 0.4 s to 48 s.
The first goal of the present work is to verify whether there is a correlation between link length d and time scale of observability W. Figure 1 Fig. 1(b,c), respectively. Figure 1 , which turns out to be significantly nonzero. Consequently, the two variables d and W turn out to be significantly correlated. Besides in the (a) part of the figure, the color map representing f d W ( , ) is shown in Fig. 1(f) where the bin size is reduced by a factor 2 on each direction. The slight asymmetry of the shape hints as well at a correlation between d and W. Figure 1(e) shows the conditional sample distribution The same distribution, upon halving again the bin size on each direction, is shown in Fig. 1(g). The shape of the most likely region suggests that the relationship between W and d is nonlinear. As explained in the Methods section, among different functional forms analyzed, a power-law of the kind ( )  www.nature.com/scientificreports www.nature.com/scientificreports/ of this analysis is shown in Fig. 2. Upon setting the normalization parameter d 0 to 75 mm (see Methods section), the parameters W 0 and γ resulting from a best-fit procedure are = . ± . W (20 9 0 2) s 0 and 0 44 0 01 γ = . ± . . The same analysis explained above was applied to each single subject in order to test whether the previous behaviour is characteristic of a single human brain or, rather, is the spurious effect of a cohort analysis. The results are shown in Fig. 3. The power-law dependence of W on d is indeed present in each subject, although with different values of the parameters W 0 and γ. Most parameters pairs are clustered in a region where W ranges from 15 mm to 30 mm and γ ranges from 0.2 to 0.7. This result is in agreement with the claim by Expert et al. 16 . As far as the exponent γ is concerned, a possible explanation of its variability relies on different levels of connectivity, as discussed in the next section. Interestingly, the average subject behaviour, in terms of average values of the two parameters among the subjects (blue dot), is in a very good agreement with the behaviour extrapolated by a pooled analysis of all subjects (red dot).

Discussion
As asserted in the Introduction, due to the presence of physical constraints, connectivity has to depend on distance 21 . This property appears to be universal, i.e. independent of the system under investigation. As an example, parallel to the work by Bullmore and Sporns 1 concerning trade-off issues in brain connectivity, the work by Gastner and Newman 22 addresses distributions of geometric properties in terms of costs and benefits within the framework of geographical networks. However, despite its universality, studies addressing the dependence of connectivity on physical distance are still few. Among these ones, a recent work by Hens et al. 23 discusses a general model for signal propagation in networks to classify them in families depending on "the interplay between network paths, degree distribution and interaction dynamics".  In neuroscience, to separate local-scale and large-scale regimes, Bellec et al. 24 described distance-dependent correlation between fMRI time series by relying on variograms, a tool from spatial statistics to quantify correlation as a function of distance 25 . Variograms allowed to empirically extract information on the spatial extent of correlations in fMRI connectivity [26][27][28] and to account for the complicated characteristics of fMRI data 29 . In addition, variograms were applied to remove spurious correlations due to voxel proximity in fMRI studies within the auditory cortex 30 and to monitor the spatial distribution of cellular activity in the brain 31 . In general, correlation is shown to quickly decrease down to a critical distance and then to saturate. A crucial issue in the investigation of human brain connectivity is to establish whether, and to what extent, structural connectivity, assessed by diffusion tensor imaging and tractography 32 , determines functional connectivity [33][34][35][36] . This issue is ultimately linked to how the neuron wiring is related to brain cognitive functions 4,37 and how it is possible to reconstruct physical links out of temporal correlations detected through electrophysiological measurements.
Our investigation tackles the problem of determining the dependence of the connectivity strength on the geometric distance in a link between two nodes. Connectivity strength is expressed in terms of time scale of observability assessed by exploiting an analytical tool recently developed that relies on the analysis of time series each stemming from a single node. In the present case, time series are cortical activities reconstructed out of MEG recordings. . Correlation diagram (left) and p-value diagram (right) for the R-TF and the R-s32 brain regions (see Table 1) computed on one recording of the second subject.

Nr.
Atlas area (hemisphere) Nr. www.nature.com/scientificreports www.nature.com/scientificreports/ What we observe, by using an approach based on the analysis of distributions similar to that discussed by Bialek et al. 38 , is that the dependence appears to be a power-law of the kind W ∼ d γ , where W is the time scale of observability and d is the geometric distance. The exponent γ takes on values ranging from 0.2 to 0.7. Lower values of γ corresponds to higher levels of connectivity, as explained in the following.
The quantity W measures the time scale at which the cross-correlation between time series generated by two nodes becomes significantly visible. The source of cross-correlation are typically peak-like 17 events that occur in both nodes at the same time. This process is countered by noise, which tends to wash out cross-correlation. If nodes are directly connected by physical links, and if the propagation speed of signals between nodes is much faster than the time scale of observability of cross-correlation -as it is the case of neural signals, which propagate in ms, whereas W is at least of order 1 s -the time scale of observability of links is not expected to depend on their length. In this case, γ is expected to be ∼ 0. On the other hand, if no direct link between two nodes exists, a working link has to rely on intemediate nodes that act as relay hubs. The relay process possibly introduces a noise component, which leads to a progressive signal degradation as distance increases. Consequently, the larger the distance d, the less frequent a peak-like co-activating process occurs, and thus the longer is the time window W required to observe peak-like events that occur at the same time.
On the basis of the large variety of neural connections occurring in the human brain, it is possible that different mechanisms like the two ones mentioned above simultaneously contribute to the observed power-law behaviour. In addition, different behaviours can be expected if particular sets of nodes, for example making up a subnetwork (like the default mode network 39,40 ), are considered, as well as if correlation between the activity of nodes is not due to a direct information link but it is rather the manifestation of simultaneous responses to a common stimulation. The approach presented in this paper can be used to identify sets of nodes that form a subnetwork, for example by characterizing them on the basis of a specific behaviour in the W ( , ) 0 γ parameter space. The exact identification and quantification of these mechanisms is beyond the scope of the present work. One possible way of tackling this issue is to study the complexity of neurophysiological signals 41 and its influence on the time scale of observability when pairs of signals are analyzed. This approach requires analytical techniques and observables typical of nonlinear time series analysis like embedding 42,43 , correlation dimension 44 , maximum Lyapunov exponent 45,46 and permutation entropy 47 .
In conclusion, we found that the the link strength, in terms of time scale of observability, significantly depends on the geometric link length. The method discussed in this work can be used to highlight the presence of an underlying subnetwork structure between subsets of nodes.

Methods observability of a link inferred out of zero-delay cross-correlation analysis of the constituent nodes.
In this work the assessment of connectivity between brain regions is carried out by applying a recently-introduced zero-delay cross-correlation method 17 . The aim of the algorithm 48 is to assess the existence of links between nodes of a possible subnetwork structure out of time series recorded at each node and to provide an estimate of the time scale on which an existing link is observable.
The input of the analysis is a pair of time series, each associated to one of the two nodes. The first step to provide an evidence of a link between the two nodes consists of evaluating the zero-delay cross-correlation between the two time series. Cross-correlation is computed as the sample Pearson correlation coefficient over moving time windows of different widths. Therefore, correlation coefficients turn out to depend on both the window position and width, as displayed by means of two-dimensional correlation diagrams shown in Fig. 4 (left). In the present work, the window width was set to span a time interval from 400 ms to 48 s.
To assess the significance of the correlation coefficients, i.e. to associate a p-value to each correlation coefficient, a surrogate-based approach is followed 49 Figure 4 (right) shows the p-value diagram corresponding to the correlation diagram displayed in Fig. 4 (left).
A p-value diagram is then further processed to assess the existence of a link between the two nodes. This step requires the evaluation of the efficiency η corresponding to a window width w, i.e. of the function η = η(w): given the p-value diagram for the pair of nodes, and given a value w of the window width, the efficiency η(w) is defined as the fraction of the running windows of width w that exhibit a p-value smaller than a given significance threshold. In this work, this significance threshold is set to 5%. Efficiency is typically a growing function of w 17 . A link between two nodes is deemed to exist at a time scale W if the efficiency at the window width W overcomes a second threshold that is here set to 0.5. If the efficiency fails to overcome the threshold, no link is attributed to the pair of nodes.
The window width W defines the minimum time scale at which a link starts to be observable: hereafter, the link is supposed to exist for any time scale larger than W, at least up to the maximum window width of 48 s. It should be noted that observability eventually fades out -so that the corresponding link disappears -once the observation window becomes so wide that the noisy contributions to the time series become dominant again 17 .
The time scale of observability of a link is a measure of the minimum observation window such that two nodes are deemed to be correlated. For example, in the case of two nodes showing an identical activity over time, the minimum time scale of observability is zero. More realistically, the activity of two nodes turns out to be co-activated only for short periods of time, for example if a subnetwork structure is established for a given purpose and (2019) 9:13412 | https://doi.org/10.1038/s41598-019-50106-2 www.nature.com/scientificreports www.nature.com/scientificreports/ then reallocated after the purpose is accomplished. In this case, correlation is observed only if a sufficiently large window is used, and the time scale of observability turns out to correspond to the repetition time of this co-activation. The value W thus provides an estimate of the time scale of the process underlying the activation of links. These time scales are not necessarily related to the information transfer speed across the subnetwork: in the case of brain networks, information transfer occurs at the millisecond scale while the activation of links spans time intervals of second or tens of seconds 17,40,50 .
Given a number N of candidate nodes of a possible subnetwork structure, there are − N N ( 1)/2 possible links, each corresponding to a node pair. The analysis described above has then to be carried on each of these pairs, and the results further processed in order to assess the possible presence of an underlying subnetwork structure 17 . Dataset and preprocessing. The dataset used in this work consists of MEG resting state recordings of 20 healthy subjects (age between 22 and 35, 16 males, 4 females) blindly extracted from the public database of the Human Connectome Project (HCP) 18,19 . The HCP provides the required ethical approval and consent needed for study and dissemination. Procedures for subject recruitment, including informed consent forms and consent to share de-identified data, were approved by the Institutional Review Board of the Washington University in St. Louis. All experimental procedures were performed under the guidelines of the HCP, which adhered to the relevant IRB processes related to that project.
In brief, for each subject, three MEG resting state sessions of about 5 minutes each are available. Data were recorded with participants lying in supine position in a whole-head 248 magnetometers MAGNES 3600 scanner (4D Neuroimaging, San Diego, CA). Participants were instructed to rest with open eyes and maintain fixation on a projected red crosshair on a dark background. MEG sensor data, sampled at 2035 Hz, were cleaned by excluding bad channels and other artifacts and removing ocular/cardiac/myogenic activity by means of independent component analysis. The public HCP database provides single-shell volume conduction models 51 computed out of a brain-enclosing surface mesh with 5000 points, as well as surface reconstructions of the mid-thickness cortical mantle, both segmented from individual anatomical T1-weighted MRI scan (Siemens Trio 3 T -Siemens Healthcare GmbH, Erlangen, Germany). All meshes coordinates are standardized to the MNI space and co-registered to the sensor array. Further details can be found on the HPC website (MEG connectome pipeline version 3.0) 19 . Cortical activity was reconstructed by means of a minimum norm algorithm 52 with unconstrained dipole orientations. We used an 8004 points cortical mesh as a source model, resulting in a grid resolution of approximately 5 mm. Noise covariance was estimated from the available empty-room recordings and no regularization was applied. Preprocessing and source reconstruction were carried out by means of FieldTrip routines 53 . After the reconstruction process, the time series were resampled from 2035 Hz down to 250 Hz. In order to get equallylong time series of 300 s duration (75000 samples), the first 4 seconds of each time series as well as a final segment of variable length were discarded.  Table 1 on a default anatomy. Colors are consistent with those used in the atlas by Glasser et al. 20 and are related to the functional group to which each area belongs. Table 1 lists the 72 brain areas randomly selected out of the 360 areas defined in the atlas by Glasser et al. 20 , which was built by combining structural, diffusion, functional and resting state MRI data from 210 healthy young individuals. The random selection was carried out by the following procedure: 1. number the areas between 1 and 360; 2. toss a number between 1 and 360 by means of a uniform random number generator and thus select the first area; 3. toss another number k between 1 and 360; 4. check whether k is equal to anyone of the previously tossed numbers; if yes, repeat operation 3, otherwise jump to the next step; 5. check whether the new area lies within 1 cm of anyone of the previously selected areas; if yes, repeat operation 3, otherwise accept the new area and jump to the next step; 6. check whether the total number of areas is less than 72; if yes, repeat operation 3, otherwise stop.
The MNI coordinates of the centroid of each area provide the locations for the 72 sources that identify the respective nodes. Figure 5 shows the anatomical position of the selected regions. For each of the 72 locations, the analyzed time series corresponds to the norm of the current dipole vector reconstructed at that location. For each pair of nodes, the geometric distance between the two nodes is computed out of their MNI coordinates. Figure 6(a) shows the distribution of link lengths evaluated by considering the 20 subjects and the 72 selected nodes for each subject. Two histograms are shown: the first one (red line) refers to the whole set of link lengths while the second one (blue line), which is more shifted to lower values than the previous one, corresponds to those links for which an assessment of the time scale of observability W provided a valid result.

Assessment of a functional relationship between W and d.
We also analyzed, for each single subject, the matching of the histogram of link lengths d corresponding to valid values of W with the histogram of the link length d, also corresponding to valid values of W, assessed on the set of all other subjects. The analysis was carried out by using the Kolmogorov-Smirnov test. In all 20 cases, the pvalue turned out to be close to unity. Link length thus follows the same distribution independently of the subject.   While the shape of the histogram shows a maximum in the center of the range as in Fig. 6(a), the presence of a frequency offset of approximately 0.013 in the histogram of Fig. 7(a) forbids the formulation of any linear relationship between W and d. On the other hand, a linear mapping of the abscissa axes appears to be possible in order to (approximately) overlap the histograms of the logarithm of the two variables d and W, as it results from the plots of Figs 6(b) and 7(b), and despite W being truncated at 48 s, or equivalently = . W log [ (s)] 3 87, because of experimental reasons. It has also to be noted that no linear mapping can lead to an overlap between d and W log( ) and, viceversa, between d log( ) and W, thus ruling out the possibility of exponential or logarithmic functional relationships between d and W.
Consequently, the dependence of the time scale of observability W on the link length d can be described by means of a power-law curve defined by ( ) where the two parameters W 0 and d 0 have the dimension of time and distance, respectively, and the exponent γ is dimensionless. To describe the power-law curve, either W 0 or d 0 can be arbitrarily set. The choice was to set d 0 to 75 mm, which approximately corresponds to the average link length (see Fig. 6).