# Whole-brain propagation delays in multiple sclerosis, a combined tractography - magnetoencephalography study

Pierpaolo Sorrentino, Spase Petkoski, Maddalena Sparaco, Emahnuel Troisi

Lopez, Elisabetta Signoriello, Fabio Baselice, Simona Bonavita, Maria Agnese

Pirozzi, Mario Quarantelli, Giuseppe Sorrentino, et al.

► To cite this version:

Pierpaolo Sorrentino, Spase Petkoski, Maddalena Sparaco, Emahnuel Troisi Lopez, Elisabetta Signoriello, et al. Whole-brain propagation delays in multiple sclerosis, a combined tractography - magnetoencephalography study. Journal of Neuroscience, 2022, pp.JN-RM-0938-22. 10.1523/jneurosci.0938-22.2022. [https://doi.org/10.1523/jneurosci.0938-22.2022](https://doi.org/10.1523/jneurosci.0938-22.2022). [hal-03851744](https://HAL.inist-cnrs.fr/hal-03851744)

HAL Id: hal-03851744

https://hal.science/hal-03851744

Submitted on 25 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Whole-Brain Propagation Delays in Multiple Sclerosis, a Combined Tractography-Magnetoencephalography Study

P. Sorrentino,¹,² S. Petkoski,¹ M. Sparaco,³ E. Troisi Lopez,⁴ E. Signoriello,³ F. Baselice,⁵ S. Bonavita,³ M.A. Pirozzi,⁶ M. Quarantelli,⁶ G. Sorrentino,⁴,⁷ and V. Jirsa¹

¹Institut de Neurosciences des Systèmes, Aix-Marseille Université, 13005 Marseille, France, ²Institute of Applied Sciences and Intelligent Systems, National Research Council, 80078 Pozzuoli, Italy, ³Department of Advanced Medical and Surgical Sciences, University of Campania Luigi Vanvitelli, 81100 Caserta, Italy, ⁴Department of Motor Sciences and Wellness, Parthenope University of Naples, 80133 Naples, Italy, ⁵Department of Engineering, Parthenope University of Naples, 80143 Naples, Italy, ⁶Biostructure and Bioimaging Institute, National Research Council, 80145 Naples, Italy, and ⁷Institute for Diagnosis and Cure Hermitage Capodimonte, 80131 Naples, Italy

Two structurally connected brain regions are more likely to interact, with the lengths of the structural bundles, their widths, myelination, and the topology of the structural connectome influencing the timing of the interactions. We introduce an *in vivo* approach for measuring functional delays across the whole brain in humans (of either sex) using magneto/electroencephalography (MEG/EEG) and integrating them with the structural bundles. The resulting topochronic map of the functional delays/velocities shows that larger bundles have faster velocities. We estimated the topochronic map in multiple sclerosis patients, who have damaged myelin sheaths, and controls, demonstrating greater delays in patients across the network and that structurally lesioned tracts were slowed down more than unaffected ones. We provide a novel framework for estimating functional transmission delays *in vivo* at the single-subject and single-tract level.

**Key words:** brain criticality; brain dynamics; brain networks; conduction velocities; magnetoencephalography; multiple sclerosis

## Significance Statement

This article provides a straightforward way to estimate patient-specific delays and conduction velocities in the CNS, at the individual level, in healthy and diseased subjects. To do so, it uses a principled way to merge magnetoencephalography (MEG)/electroencephalography (EEG) and tractography.

## Introduction

Higher-level cognitive functions might be emergent and depend from the regulated interaction of multiple brain areas (Stam, 2014). Hence, network theory has been used to describe such large-scale functional interactions within the brain (Bullmore and Sporns, 2009), as well as the underlying structural

connections (Watts and Strogatz, 1998; Barabasi and Albert, 1999). From the functional standpoint, most studies used functional magnetic resonance imaging (fMRI), thereby focusing on the slow-evolving activity (typically in the timescale of seconds) and using assumptions of stationarity. The presence of time-averaged spatial patterns of correlated and anti-correlated activity was consistently reported (Fox et al., 2005), which also related to the structural topology (Honey et al., 2009; Tewarie et al., 2019). Per contra, brain activities do not exclusively generate periodic patterns, but also aperiodic ones emerge on the large scale (Zalesky et al., 2014). Structural topology imposes constraints on the evolution of brain dynamics, as shown by the fact that the probability of two brain regions activating sequentially is proportional to the coupling intensity along the brain tract linking them (Sorrentino et al., 2021b). However, structural topology alone is necessary, but not sufficient, to generate the observed dynamics, as delays, and noise, are also key elements (Niebur, 1991; Atay et al., 2004; Ghosh et al., 2008; Deco et al., 2009; Cabral et al., 2014; Petkoski and Jirsa, 2019; Tewarie et al., 2019). The time it takes two gray matter regions to sequentially activate depends

Received May 16, 2022; revised Sep. 13, 2022; accepted Sep. 19, 2022.

Author contributions: P.S., S.P., F.B., and V.J. designed research; P.S., M.A.P., and M.Q. performed research; P.S. contributed unpublished reagents/analytic tools; P.S. analyzed data; P.S. and S.P. wrote the first draft of the paper; P.S., S.P., M.S., E.T.L., E.S., F.B., S.B., M.A.P., M.Q., G.S., and V.J. edited the paper; P.S. wrote the paper.

This work was supported by the University of Naples Parthenope "Ricerca locale" grant and by European Union's Horizon 2020 Research and Innovation Program under Grant Agreement No. 945539 (SGA3), Human Brain 570 Project, and VirtualBrainCloud No. 826421.

The authors declare no competing financial interests.

Correspondence should be addressed to P. Sorrentino at pierpaolo.sorrentino@univ-amu.fr or V. Jirsa at viktor.jirsa@univ-amu.fr.

https://doi.org/10.1523/JNEUROSCI.0938-22.2022

Copyright © 2022 Sorrentino et al.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
on multiple factors, such as the distance separating them, the properties of the white matter bundle (structural tract) linking them (such as its diameter and myelination), and the overall topology of the network embedding them (Purves et al., 2001). For this reason, estimating a comprehensive map of the delays poses great challenges, and most modeling studies assume homogeneous velocities of conduction (typically ranging from 2 to 10 m/s; Sanz-Leon et al., 2015; Breakspear, 2017). Hence, the delays are conceptualized as purely dependent on distance (Ghosh et al., 2008; Petkoski et al., 2018; Petkoski and Jirsa, 2019). This is a suboptimal simplification, since we know (e.g., from routine examinations such as visual evoked potentials) that the velocities in highly myelinated, thick tracts can reach up to 150 m/s (Purves et al., 2001). Combining magnetoencephalography (MEG) and tractography one can track the beginning and the spread of a regional perturbation (Sorrentino et al., 2021b). Building on this framework, we hypothesized that the time it takes a signal to go from one region to another might be a candidate proxy for the delay of transmission across any given structural tract. Shorter delays might relate, among other factors, to diameter and myelination. This should result in perturbations spreading with velocities that are not constant but, rather, rise as a function of the length of the white-matter tracts (given that longer tracts are typically also the thicker and more myelinated ones). Conversely, we expected the width of the distribution of the delays to be much narrower than what would be expected from the length of the white matter bundles alone (since the longer tracts are hypothesized to be the fastest). We further reasoned that, if our measurements are related to myelination, then lesions in the myelin sheath should increase the delays. To test this hypothesis, we used source-deconstructed MEG signals to estimate the delays, defined as the average time it takes a perturbation generated in region *i* to spread to region *j*, and then estimated the lengths of the white matter bundles in a cohort of 18 patients affected by multiple sclerosis and 20 healthy subjects using tractography. We used the delay in the subsequent recruitment of brain regions as a proxy for the conduction delay of an impulse between regions and estimated the corresponding velocities by dividing the tract lengths by the corresponding delays. First, we explored the distributions of the delays and velocities in healthy humans. Then, we evaluated the effect of the damage to the myelin sheath on the delays. To further verify the whole workflow, we applied the same procedure to a publicly available dataset based on combined electroencephalogram (EEG) recordings and tractography, and we were able to extract subject-specific delays and, thus, velocities, which had statistics that were similar to those obtained via MEG. The findings were further verified on an independent MEG/Diffusion Tensor Imaging (DTI) dataset, using a different parcellation and a different source-reconstruction algorithm.

## Materials and Methods

### Participants

The participants (of either sex) were recruited at the outpatient clinic of the Institute for Diagnosis and Cure Hermitage Capodimonte. The diagnosis of multiple sclerosis was made according to the revised 2017 McDonald criteria (Thompson et al., 2018a). Exclusion criteria were age <18 years, clinical relapse and/or steroid use in the three months before the study, inability to understand and complete “patient reported outcomes” and cognitive evaluation, or inability to undergo the MRI scan. All patients underwent a neurologic clinical examination, Expanded Disability Status Scale (EDSS) scoring, the Symbol Digit Modalities Test (SDMT) to measure cognitive impairment, the Fatigue Severity Scale (FSS), and the Beck Depression Inventory (BDI). The controls for the MS cohort were selected from among the caregivers and spouses of the

<table><caption>**Table 1.** Features of the multiple sclerosis cohort</caption><thead><tr><th></th><th>Controls</th><th>MS</th><th>p-value</th></tr></thead><tbody><tr><td>Age (years)</td><td>45.8 (±11)</td><td>44.9 (±9.9)</td><td>0.8</td></tr><tr><td>Education (years)</td><td>13.6 (±3.8)</td><td>13.8 (±5)</td><td>0.9</td></tr><tr><td>Gender (m/f)</td><td>6/14</td><td>6/12</td><td>0.3</td></tr><tr><td>Disease duration (months)</td><td>—</td><td>187.7 (±131.8)</td><td>—</td></tr><tr><td>EDSS</td><td>—</td><td>4.5 (±1.9)</td><td>—</td></tr><tr><td>SDMT</td><td>—</td><td>40.3 (±13)</td><td>—</td></tr><tr><td>FSS</td><td>—</td><td>36.1 (±14)</td><td>—</td></tr><tr><td>BDI</td><td>—</td><td>12.8 (±1.3)</td><td>—</td></tr><tr><td>LL</td><td>—</td><td>12,959 (±12,253) mm³</td><td>—</td></tr></tbody></table>

EDSS, Expanded Disability Status Scale; SDMT, Symbol Digit Modalities; FSS, Fatigue Severity Scale; BDI, Beck Depression Inventory; LL, lesion load.

patients. Genetic relatives were not allowed as controls. The subjects for the second, larger healthy cohort were selected as described previously (Sorrentino et al., 2021b). The independent EEG/DTI dataset is described in (Schirner et al., 2018). The demographics and main clinical and radiologic features of the MS cohort are summarized in Table 1. The study was approved by the local Ethics Committee (Prot.n.93C.E./Reg. n.14-17OSS).

### MRI acquisition and processing

Each MRI scan was performed immediately after the MEG recording on the same MRI scanner (1.5 Tesla, Signa, GE Healthcare). Analyzed sequences included echo-planar imaging for DTI reconstruction (TR/TE 12 000/95.5 ms, voxel 0.94 × 0.94 × 2.5 mm³, 32 equally spaced diffusion-sensitizing directions, 5 B0 volumes) and 3D-FLAIR volume for WM lesion segmentation (TR/TE/TI 7000/145/1919 ms, echo train length 170, 212 sagittal partitions, voxel size 0.52 × 0.52 × 0.80 mm³). Preprocessing of the diffusion MRI data were conducted using the software modules provided in the FMRIB software library (FSL; http://fsl.fmrib.ox.ac.uk/fsl). All diffusion MRI datasets were corrected for head movements and eddy current distortions using the “eddy correct” routine (Smith et al., 2004), rotating diffusion sensitizing gradient directions accordingly (Leemans and Jones, 2009), and a brain mask was obtained from the B0 images using the Brain Extraction Tool routine (Smith, 2002). A diffusion-tensor model was fitted at each voxel, and fiber tracts were generated over the whole brain by deterministic tractography using the Fiber Assignment by Continuous Tracking (FACT) algorithm implemented in Diffusion Toolkit (angle threshold 45°, spline-filtered, masking by the FA maps thresholded at 0.2). Two cortical study-specific ROI datasets were obtained by masking the ROIs available in the AAL atlas (Tzourio-Mazoyer et al., 2002) and in an MNI space-defined volumetric version of the Desikan-Killiany-Tourville (DKT) region of interest (ROI) atlas (Klein and Tourville, 2012) using the gray matter (GM) probability map available in the Statistical Parametric Mapping software package (SPM), thresholded at 0.2. This was done to limit the endpoints of the fibers to cortical and adjacent juxtacortical white matter voxels in the subsequent ROI-based analysis of the tractography data. The analysis was replicated twice to ascertain how robust the method was when used with specific brain parcellations. To obtain the corresponding patient-specific ROI sets, each participant’s FA volume was spatially normalized (Friston et al., 1995) to the FA template provided by FSL using SPM12, and the resulting normalization matrices were inverted and applied to the two ROI sets. Additionally, for each subject, the MS lesion map was obtained by segmenting the 3D-FLAIR volume using the lesion prediction algorithm (Schmidt et al., 2012) implemented in the Lesion Segmentation Tool (LST toolbox version 3.0.0; www.statistical-modeling.de/lst.html) for SPM. The 3D-FLAIR volume was then co-registered to the EPI of the patient (Ashburner and Friston, 1997), and the coregistration matrix was applied to the corresponding WM lesion volume, which was resampled by nearest-neighbor interpolation, thus obtaining the patients’ lesion masks, which would be coregistered with the DTI volume. Finally, for each patient the average length of the fibers connecting each pair of ROIs and whether the voxels crossed by those fibers had MS lesions were calculated separately for the AAL and DKT ROI sets, using an in-house routine written in Interactive Data Language (IDL, Harris Geospatial Solutions). To this end, each tract length was computed, for each subject, as the mean of
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top left, The structural connectome. From the tractography, we obtained an estimate of each tract length. Averaging across subjects provided a group-level estimate. Top right, Source-reconstructed MEG time-series. The time points that are below a threshold are represented in blue, while the timepoints above the threshold are in red (the red dots have been magnified). In the example, after region *i* rose above the threshold (and, hence, the neuronal avalanche had started), it took region *j* three timesteps to be recruited and five timesteps for region *z*. The entries of the matrix are expressed in samples. Delays in the manuscript are expressed in seconds. Bottom, Each subject-specific length matrix was divided element-wise by the corresponding estimated delay to obtain a subject-specific velocity matrix.</figcaption></figure>

the physical distances covered by the fibers composing the tract. These in turns were calculated as the sum of the Euclidean distances between the fiber subsequent turning points, including the two extremities, which are the output of the FACT algorithm.

### MEG preprocessing

MEG preprocessing and source reconstruction were performed as previously described (Sorrentino et al., 2018). Preprocessing and source reconstruction operations were carried using the Fieldtrip Toolbox (Oostenveld et al., 2011). Each participant underwent a MEG recording, composed of both eyes-closed resting-state segments of 3'30" each. Four anatomic coils were applied on the head of each participant and their position was recorded along with the position of four head anatomic points, to identify the position of the head during the recording. Eye blinking (if present) and heart activity were recorded through electro-oculogram (EOG) and electrocardiogram (ECG), to identify physiological artifacts (Gross et al., 2013). An expert rater checked for noisy signals and removed them. An anti-alias filter was applied to the MEG signals, acquired at 1024 Hz, before being filtered with a fourth order Butterworth IIR bandpass filter (0.5–48 Hz). We used principal component analysis (Sadasivan and Narayana Dutt, 1996; de Cheveigné and Simon, 2008) and supervised independent component analysis (Barbati et al., 2004) to remove the environmental noise and the physiological artifacts (recorded with EOG and ECG), respectively.

### Source reconstruction

Signal time series were reconstructed using both the AAL and DKT atlases (Tzourio-Mazoyer et al., 2002; Gong et al., 2009), which consist of 116 and 84 ROIs, respectively. The reconstruction took place using the volume conduction model proposed by Nolte (Nolte, 2003). Based on the native MRIs of each subject, the linearly constrained minimum variance (LCMN; Hillebrand et al., 2016) beamformer was applied to reconstruct the signal sources based on the centroids of each ROI

(Hillebrand et al., 2016). ROIs belonging to the cerebellum were excluded because of the low reliability of their source reconstruction (Lardone et al., 2018), for a total of 90 ROIs in the AAL atlas and 84 ROIs in the DKT atlas.

### Neuronal avalanches and branching parameter

To study the dynamics of brain activity, we based our analysis on “neuronal avalanches.” First, the time series for each ROI was downsampled to 512 Hz and discretized by calculating the *z* score; then the positive and negative excursions beyond a threshold were identified (Fig. 1). The main results reported here refer to a threshold equal to three standard deviations (|*z*| = 3), but thresholds of 2.5 and 3.5 were also tested. A neuronal avalanche began when at least one ROI went above the threshold (|*z*| > 3) and ended when all the ROIs were below the threshold (Beggs and Plenz, 2003; Shriki et al., 2013). Before proceeding with the analyses, we binned the data, to ensure that we captured any critical dynamics, if present. To estimate the suitable time bin length, for each subject, for each neuronal avalanche, and for each time bin duration, the branching parameter *σ* was estimated (Harris, 1964; Haldeman and Beggs, 2005). In fact, systems operating at criticality typically display a branching ratio of ~1 (Chialvo, 2010). The branching ratio was calculated as the geometrically averaged (over all the time bins) ratio of the number of events (activations) between the subsequent time bin (descendants) and that in the current time bin (ancestors) and then averaging the avalanche-specific branching ratio over all the avalanches (Bak et al., 1987). More specifically:

$$ \sigma_i = \prod_{j=1}^{N_{\text{bin}}-1} \left( \frac{n_{\text{events}}(j+1)}{n_{\text{events}}(j)} \right)^{\frac{1}{N_{\text{bin}}-1}} \quad (1) $$

$$ \sigma = \prod_{i=1}^{N_{\text{aval}}} (\sigma_i)^{\frac{1}{N_{\text{aval}}}}, \quad (2) $$

where $\sigma_i$ is the branching parameter of the *i*th avalanche in the dataset, $N_{bin}$ is the total number of bins in the *i*th avalanche, $n_{events}(j)$ is the total number of events active in the *j*th bin, and $N_{aval}$ is the total number of avalanches in the dataset. We tested bins from 1 to 5 and picked three for further analysis.

### Estimation of delay matrices and velocity matrices

The delays were estimated for each avalanche. The procedure is schematized in Figure 1. In an avalanche, from the moment region *i* activated, we recorded how long it took region *j* to activate. These are what we considered to be delays. Hence, for each avalanche we obtained a matrix, in which the rows and columns represented brain regions and the entries contained the delays. We then averaged across all the avalanches belonging to one subject, obtaining an average *ijth* delay. The average was performed disregarding zero entries, since each avalanche-specific matrix is very sparse. With this procedure, a subject-specific delay matrix was built. Averaging across subjects (again discarding zero entries) yielded a group-specific matrix. Note that it can happen that no delay estimation is available for a specific edge in a specific subject, since the estimation depends on the fact that an avalanche propagated across the edge. Finally, we obtained a velocity matrix, dividing the tract lengths, obtained from MRI, by the corresponding delays. Velocities were only obtained for the edges where a direct structural connection was present.

### Experimental design and statistical analyses

To build null models to test the functional delays estimates, we randomized the temporal order of the avalanches without changing the spatial structure. To this end we built a surrogate dataset where the avalanche time points were randomly shuffled, but the regions recruited at any
<figure><img src="image_3.png" /><figcaption>**Figure 2.** A, Top, Tract lengths. Rows and columns represent brain regions, and the color code conveys the corresponding length of the tract linking any two regions. Middle, Group average of the delays. Rows and columns represent brain regions, while the color code represents the average time it took region *j* to become active provided region *i* had been active earlier. Values are in seconds and reported on a log-scale to highlight the texture of the delay matrix, which appears to be strongly correlated to the length matrix. Bottom, As before, brain regions are represented as rows and columns, while the matrix entries represent the velocities, expressed in meters/s. B, Top, middle, and bottom, Histograms corresponding to the matrices for the lengths, delays, and velocities, respectively. C, The glass-brain provides an overview on the topography of the functional edges with fast functional velocities. In particular, the blue dots represent brain regions, and the edges signify the velocities. Only structural connections with transmission speed higher than 15 m/s are shown. The width of the edge is proportional to its velocity.</figcaption></figure>

given time point were fixed. By doing this, the time-structure of the
recruitment of regions was disrupted, but the purely spatial component
was retained. After the permutations, the average delays were again com-
puted for each edge. The procedure was repeated 1000 times, yielding a
hundred surrogate delays for each edge. We then used the delays esti-
mated from the random surrogates to compute functional velocities
(dividing the length of the structural edges by the surrogate delay).
Finally, we computed the growth of the velocities as a function of the
length of the structural tracks and compared this with the observed dis-
tribution. To test differences between the distribution of the delays in
the healthy subjects and patients, we used the Kolmogorov–Smirnov
(KS) test. To perform the edge-wise comparison of the delays in the
healthy versus lesioned edges in MS patients, we used permutation test-
ing (Nichols and Holmes, 2002). In short, we tested the null hypothesis
that lesions in the edges would not have an impact on the delays. First,
we calculated the average edge-wise difference between the delays in
each patient and the average delay in the corresponding edge in the con-
trols. Then, we randomly selected a subset of the differences in the delays
with the same size as the number of lesioned edges and computed its av-
erage. We repeated this procedure 1000 times, building a distribution of
the differences in the delays that are to be expected by randomly drawing
a subsample of the edges. Finally, we compared this distribution to the
observed difference in the lesioned edges to obtain the probability of
observing the data under the hypothesis that the edge-specific lesions
would not slow down the functional transmission.

Data and materials availability

Code is available on GitHub. The MEG data and the reconstructed ava-
lanches are available upon request to the corresponding author, conditional

on appropriate ethics approval at the local site. The availability of the
data was not previously included in the ethical approval, and therefore
data cannot be shared directly. In case data are requested, the corre-
sponding author will request an amendment to the local ethical commit-
tee. Conditional to approval, the data will be made available.

Results

Delay estimation

In this study, we noninvasively estimated the functional delays
in transmission across the network of white matter bundles *in vivo*.
To this end, we combined source-reconstructed magneto-
encephalography and tractography. An overview of the pipeline
is shown in Figure 1. The upper row of Figure 2 shows the aver-
age structural matrix for the controls (i.e., the tract lengths)
with the corresponding distributions of the tract lengths on the
right. The middle and bottom rows of Figure 2, respectively,
show the matrices and distributions of the delays and the veloc-
ities. Importantly, the structure of the delay matrix is related to
the length of the structural tracts. However, the width of the
distribution of the delays is much narrower than what would be
expected given constant propagation velocities (see below).
When estimating the velocities, we found that a consistent, fat-
tailed distribution emerged. As expected, because of the hetero-
geneity of the tracts, the velocities appeared to be far from
homogeneous, ranging from ~2 to ~60 m/s. This range appears
to be in accordance with known velocities of major myelin tracts
<figure><img src="image_4.png" /><figcaption>**Figure 3.** A, Group-level relationship between tract lengths and functional delays. Blue dots represent iso-hemispheric edges, while red dots represent cross-hemispheric edges. B, For each of the 20 healthy controls, the red dot represents the intensity of the relationship between the tract lengths and delays, calculated using Spearman's correlation coefficient. For each subject, the avalanches were randomized by shuffling the time sequence while preserving the spatial pattern of the active regions within each timestep. Based on these newly randomized avalanches, the edgewise average delay was calculated and then related to the tract lengths. The procedure was repeated 1000 times per subject. The resulting Spearman's correlations are represented as blue dots. C, For each subject, the delays were averaged according to the percentile of the corresponding length. Hence, each blue dot represents the average delay across all edges whose length belonged to the (patient-specific) nth percentile. The colored lines show the delays that would be expected if the velocities were homogeneous across all edges. In this case, the growth of the delays followed the growth of the tract lengths. D, The red dots represent the observed velocities (averaged according to the percentile of the tract length). The red line represents the average velocity estimated based on the delays derived from the randomly shuffled avalanches. E, Group-level relationship between the velocities. Blue dots represent iso-hemispheric edges, while red dots represent cross-hemispheric edges.</figcaption></figure>

(Waxman and Bennett, 1972; Waxman, 1980; Purves et al., 2001). When plotting the distribution of the fastest edges (see glass brains in Fig. 2), it appears that these were not evenly distributed across the brain and that the cross-hemispheric edges were preferentially selected as the fastest ones, a finding which agrees with the literature (Shen et al., 2015). We found that the delays grew as a function of the length of the tracts, as shown by the Spearman correlation between edge length and delays ($r=0.35, p < 0.0001$; Fig. 3A). The results held using Pearson's correlation ($r=0.31, p < 0.0001$). This relationship held at the individual level (Fig. 3B). However, while the delays were related to the tract lengths, they were not only determined by the lengths. In fact, the tract lengths ranged across an order of magnitude, whereas the corresponding increase in the delays was only moderate. For comparison, in Figure 3C, we show the expected delays given constant velocities and reveal that the observed delays were remarkably steady despite the difference in the lengths of the structural tracts.

## Statistical validation

First, we compared the delays retrieved from the random timeseries to the tract lengths, obtaining 100 Spearman's *rs*, to which we compared the observed correlation, *p* < 0.001. We then used the random surrogates to compute random velocities, i.e., dividing each edge length by the corresponding delay derived from the random surrogates and, using the same validation scheme, we showed that the random velocities appeared to be more

strongly related to the tract lengths than the observed velocities (*p* < 0.001). In other words, when we divided the lengths by the randomized delays, the resulting velocities became a function of the tract length alone. On the other hand, when we divided the tract lengths by the observed delays, the longer tracts appeared to be faster than if the delays were only a function of distance. The mean (red line) and upper and lower bounds (shaded area) of the surrogate velocities derived from surrogate delays (grouped by the percentile of the corresponding tract length) are shown in Figure 3D. The delays derived from the surrogate data led to a higher estimate of the velocities compared with the observed ones. Finally, as shown in Figure 3E, we confirmed, using real data, that the transmission velocities grew as a function of the length of the tracts, such that the longer tracts were also the faster ones. All in all, this part of the analysis showed a finely regulated relationship between the delays and the lengths of the structural tracts, which implies nonhomogeneous functional velocities.

## Delay estimation in multiple sclerosis patients

Next, we sought to test our framework in patients affected by multiple sclerosis, which is a prototypical disease in which myelin in the CNS is selectively attacked by the immune system (Thompson et al., 2018b). We excluded patients with severe functional impairment (Expanded Disability Status Scale, <7) so that we could obtain a picture that is likely to be influenced more by demyelination than by degenerative phenomena (Frischer et al., 2009). We expected to observe greater delays (and lower
A

<figure><img src="image_5.png" /></figure>

B

<figure><img src="image_6.png" /></figure>

C

EDGE-WISE DELAY DIFFERENCES

<figure><img src="image_7.png" /></figure>

D

<figure><img src="image_8.png" /><figcaption>**Figure 4.** A, Each dot represents the average delay for each percentile of tract length. Blue dots represent controls, while red dots represent patients. B, Empirical cumulative distributions of the delays. The blue line corresponds to the average delays of the controls, while the red lines correspond to the average delays of the MS population. C, The violin plot shows the edgewise differences between delays in each patient and the average delay across all the controls for the corresponding edge. D, The vertical red line marks the observed average edgewise difference in the delays calculated only for delays that were lesioned. The distribution to the left shows the average delay difference observed after selecting a random sample based on a thousand randomizations of edgewise delay differences (i.e., ignoring the information about the structural integrity), with the size of each random sample equal to the number of lesioned edges. The results enabled us to reject the null hypothesis that edge-specific lesions would not affect the delays.</figcaption></figure>

velocities) in the patient population compared with the controls. As shown in Figure 4A, the average delay per percentile was consistently higher in the multiple sclerosis patients compared with the controls. The empirical cumulative distribution function confirmed the difference in the delays in patients compared with the controls (KS test, *p* < 0.001; Fig. 4B). Note that, in this analysis, similar to what we did before, we grouped the delays according to percentile across the whole group. Hence, the average value refers to the average delay of all the edges that fall within a given length percentile, with the percentile defined for each subject. We then moved on to an edge-wise comparison. To this end, for each patient, the average delay of each edge was compared with the average delay of each corresponding edge in the controls. Hence, we obtained the difference between the edge-specific delay and the corresponding average delay in the controls for each patient. To base our findings on a more stable estimate, we considered only edges for which a delay estimate was available for each of the 20 controls (i.e., a total of 1202 edges). Figure 4C shows the distribution of the differences between the delays (patients minus controls). The distribution is not centered around zero, as would be expected if the patients did not have longer delays. Instead, it is evident at the visual inspection that the distribution is heavily skewed toward positive values, implying longer delays in the patients compared with controls. Finally, we investigated edge-specific lesions. In fact, while it is reasonable to expect globally greater delays in the patients compared with the controls, the delays corresponding to edges that were lesioned might be more lengthened compared with the delays corresponding to healthy edges. For each patient, we classified edges as healthy or lesioned, based on the presence or absence of structural damage (i.e., the tracts were considered as lesioned if any of the corresponding voxels was lesioned). Under the null hypothesis that the delays corresponding to structurally lesioned edges do not differ from the delays corresponding to nonlesioned edges, we calculated the average delay difference in random samples of nonlesioned edges (in patients), with the size of the sample equal to the number of lesioned edges. We repeated this procedure 1000 times and compared the obtained distribution to the observed difference in the lesioned edges. The results are shown in Figure 4D. The delays in the lesioned edges slowed more (with respect to the corresponding delays in the healthy population) than the healthy edges.

Hence, selecting a subset of edges based on the structural information, we retrieved a difference in the temporal structure of the functional dynamic. This difference would not be expected if the structural damage was unrelated to the delays, as shown by the permutation analysis (*p* < 0.001).

## Replication in independent datasets

The results of the delay estimation were tested using an independent dataset based on co-registered MRI and EEG, and all the main findings were confirmed (Fig. 5). Note that the cleaning and source-reconstruction algorithms were different in this dataset as compared with ours, showing robustness to both the technique used and data processing. Furthermore, we replicated the main results using different parcellations (AAL and DKT), varying the *z* score threshold (2.5 and 3.5 in addition to the 3.0 reported here), varying the binning (see Materials and Methods), and using yet another source-reconstruction algorithm (i.e., the “residual variance” method). The residual variance is a particular dipole-fit approach that involves the minimization of the signal that remains unexplained by a current source model in which dipoles are assumed to have fixed position and fixed or varying orientation. Conversely to LCMV, there is no linear constrain involved in the minimization process (Grech et al., 2008). Furthermore, a larger MEG/tractography dataset, involving 47 young healthy subjects, was also used to further explore the robustness of our findings. These results were replicated using the AAL atlas. Finally, we varied the minimum size of the avalanches used to compute the delays, from taking all avalanches into account, up to selecting only those longer than 10/15 samples (analyses shown in Fig. 5). Finally, we confirmed that the delay matrices showed convergence as the sample size increased, as a further check on the validity of our results (Fig. 5).

# Discussion

## The topochronic map

In this study, we developed a novel method to quickly and non-invasively obtain an estimate of the time delays that occur when locally generated perturbations spread to other brain regions. Our methodology was able to capture the presence of homogeneous delays across brain regions, although the length of the structural tracts grew according to a fat-tail distribution (Sporns and
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Left, Distribution of delays (log-scale). *A*, Delays estimation with binning = 2. *B*, Delays estimation with binning = 4. *C*, Delays estimation using avalanches with size >10. *D*, Delays estimation using avalanches with size >15 bins. *E*, Delays estimation using all available avalanches. *F*, Replication of the delays estimation based on an independent cohort, based on source-reconstructed EEG and tractography. *G*, Delays estimation with a z score >±[2.5]. *H*, Replication of the delays estimation with a z score >±[3.5]. *I*, Delays estimation using the “residual variance” as a source-reconstruction algorithm. *J*, Delays estimation based on an independent cohort of 47 young healthy subjects, based on source-reconstructed MEG and tractography (parcellation according to the AAL atlas is shown). Right, Analysis of convergence. On the x-axis, the image shows the number of patients over which the delays have been averaged. The y-axis shows the mean difference between the average delays (i.e., across all edges) estimated in x subjects and the average delays estimated in x-1 subjects. The differences converge around zero as the delays were averaged across a growing number of participants.</figcaption></figure>

Zwi, 2004; Sporns et al., 2005). If the delays were only dependent on the tract lengths, they should be longer for longer tracts (i.e., it should take longer to cover a greater distance). This effect should be major, given that longer white-matter tracts are roughly one order of magnitude longer than shorter ones (Bullmore and Sporns, 2009). The fact that the delays do not scale with increased tract lengths implies one or more compensatory mechanisms. Such mechanisms might involve varying axonal diameters, myelination, and network effects. In particular, myelination greatly affects the velocity of propagation (Waxman and Bennett, 1972; Waxman, 1980; Horowitz et al., 2015). In fact, given our estimated delays, we found that, in the healthy subjects, the velocities ranged from ~3 to ~60 m/s. This observation is in sharp contrast with the typical simplifying assumptions, made in modeling studies, that velocities are constant, which would be expected to lead to a broad range of delays across the network. The fact that the transmission between the regions connected by the longest edges is faster than those with shorter edges is likely important in terms of the unfolding of the dynamics and should not be surprising considering the experimental results on stimulation, although until this current study solid measurements of personalized links on the whole-brain level had been missing. Hence, noninvasive estimates of the temporal constraints *in vivo* are highly relevant for modeling individual dynamics. Furthermore, many of the fast edges connect regions that are highly central to the brain network, thus, likely causing reverberations that greatly and nonlinearly impact the average transit time (He et al., 2007). Although we can only speculate, we think that it is possible that homogeneous delays, that is, having perturbations reach wide-spread brain regions or arrive at a focal point from wide-spread regions simultaneously, would be favorable for allowing the brain to have simultaneous access to information from multiple locations across the brain, a concept which fits well into the framework of global workspace theory (Dehaene et al., 1998). This supposition finds support in that white matter damage impairs conscious access in multiple sclerosis patients (Reuter et al., 2009; Rosanova et al.,

2012). A reduction in the complexity of the spatiotemporal spreading of such perturbations (referred to as “neuronal avalanches” within the framework of critical dynamics) has been shown to be related to states of reduced consciousness (Tononi and Sporns, 2003) as well as to neurodegeneration (Rucco et al., 2020; Sorrentino et al., 2021a). Importantly, our analyses focused on rare, intermittent, large-scale bursts of activations, which have been consistently observed in human brains (Shriki et al., 2013). The importance of such rapid transients to large-scale brain dynamics is confirmed by recent findings showing that the patterns of functional connectivity are shaped by specific, short moments in time (Zamani Esfahlani et al., 2020). Furthermore, the fact that avalanches preferentially spread along structural tracts indicates that functional delays might be used as a proxy to estimate the velocities across individual tracts (Sorrentino et al., 2021b). As we investigated the topographic distribution of the fast edges, we found that they were spatially distributed nonhomogeneously. Importantly, fast connections seem to be preferentially cross-hemispheric, a finding that is not surprising from a neuroanatomical standpoint, provided that these connections are mediated, for example, via the corpus callosum (Caminiti et al., 2013). While our data were based on broad-band data, the patterns that emerged were not dominated by the occipital α frequency, as it is often the case with M/EEG data. However, we wish to stress that such results might have been biased by the tractography, which may have preferentially estimated tracts in specific anatomic regions (Shen et al., 2015). Hence, further validation is needed to confirm this finding, which should be considered explorative. The replication of our results in an independent MEG dataset as well as in a publicly available multimodal EEG-MRI dataset (Schrirner et al., 2018) indicates the reliability of our findings by increasing the probability that they are not modality specific.

## Implications in neurologic diseases and multiple sclerosis

As explained, myelination is believed to greatly influence conduction velocity and, hence, to modulate delays. To test this, we
measured the delays in multiple sclerosis patients. As expected, the delays were greater in multiple sclerosis patients than in matched healthy controls. Selecting MS patients at a fairly early stage should make the role of demyelinating lesions more prominent than the role of degeneration (Stadelmann et al., 2011). However, this was not quantified and remains a potential source of confounds that needs further investigation. To provide edge-specific information about the time delays, we used subject-specific lesion masks to separate lesioned edges from nonlesioned ones in the patients. We focused the analysis on the delays to avoid potential biases induced by unreliable length estimates because of the demyelinating lesions (Delettre et al., 2019). We observed, for both the lesioned and nonlesioned edges, that functional delays increased in the multiple sclerosis patients, in accordance with the hypothesis that damage to the myelin would provoke longer delays. This also provides support for the claim that our measurements are related to temporal patterns of delays that are imposed on the overall structural connectome. The fact that avalanches propagate more slowly along lesioned edges shows the relevance of the direct pathways in determining the edge-specific delays given the expected relationship between the structural integrity of the tract and the velocity of propagation along it. However, it is important to stress that the role of network effects cannot be easily disentangled (Bullmore and Sporns, 2009). In fact, the transmission in patients was slower even across unaffected tracts. This might be interpreted as an expression of the fact that the delays likely depend from a combination of both direct and indirect paths through which a perturbation can potentially travel between two regions. In this sense, one does not expect two regions that are linked by a healthy edge that is embedded in a diseased network to communicate as quickly as two regions that are also linked by a healthy edge as well as embedded in a healthy network. Other contributing factors may include the erroneous classification of damaged but subthreshold edges as healthy.

## Cross-dataset validation and limitations

The strength of this study is that the tract lengths and the delays were estimated using two different techniques, making it unlikely that the relationship is spurious or tautological. We tested our results by changing both the binning parameter (see Materials and Methods) and the z score threshold to define the spreading of the perturbations and the brain parcellation, again showing the robustness of our findings. Similar results were obtained using EEG, using a different preprocessing pipeline, applying a different algorithm to source-reconstruct, and performing a further parcellation (see Fig. 5). The fact that the relationship between length and velocity was maintained at the subject level is remarkable. One limitation that should be considered, however, is the fact that we used the DKT and the AAL atlases, both of which are coarse grained. However, finer grained parcellations, while optimal for structural MRI, would have been below the resolution for MEG and, hence, might have created spurious results. We propose that, beyond the topology of functional connections imposed by the spatial scaffolding (Sorrentino et al., 2021b), the conjugate property of connectivity is temporal in nature and complementary to the structural topology (Petkoski and Jirsa, 2020). Together, spatial and temporal constraints reveal the topochronic framework from which oscillatory brain activity emerges. Applying tools from statistical mechanics and dynamical system theory to “reverse engineer” the pattern of delays that occur in the living brain, we noninvasively measured the large-scale pattern of functional delays that

occurs in the human brain at rest. The pattern of delays is basically similar across individuals. Hence, the time it takes an (internally generated) local impulse to affect other regions is not only a function of the tract length but is also heavily modulated by the properties of the tract itself and, globally, by the large-scale structure of the network. Including subject-specific delays provides the potential to improve virtual, personalized in-silico brain models (Hashemi et al., 2020). Importantly, the fact that delay estimation can be obtained from EEG, a widely available technique, allows nearly every facility to include subject-specific delays into personalized brain models. Furthermore, large cohorts will be needed to get normative data from large, stratified cohorts, to create priors to be included in personalized models, in case subject-specific estimations are not available. In conclusion, we proposed a simple method for combining multi-modal imaging within the framework of statistical mechanics, to derive subject-specific topochronic maps of large-scale brain dynamics in both healthy and diseased populations.

## References

* Ashburner J, Friston K (1997) Multimodal image coregistration and partitioning—a unified framework. *Neuroimage* 6:209–217.

* Atay FM, Jost J, Wende A (2004) Delays, connection topology, and synchronization of coupled chaotic maps. *Phys Rev Lett* 92:144101.

* Bak P, Tang C, Wiesenfeld K (1987) Self-organized criticality: an explanation of the 1/f noise. *Phys Rev Lett* 59:381–384.

* Barabasi AL, Albert R (1999) Emergence of scaling in random networks. *Science* 286:509–512.

* Barbati G, Porcaro C, Zappasodi F, Rossini PM, Tecchio F (2004) Optimization of an independent component analysis approach for artifact identification and removal in magnetoencephalographic signals. *Clin Neurophysiol* 115:1220–1232.

* Beggs JM, Plenz D (2003) Neuronal avalanches in neocortical circuits. *J Neurosci* 23:11167–11177.

* Breakspear M (2017) Dynamic models of large-scale brain activity. *Nat Neurosci* 20:340–352.

* Bullmore E, Sporns O (2009) Complex brain networks: graph theoretical analysis of structural and functional systems. *Nat Rev Neurosci* 10:186–198.

* Cabral J, Luckhoo H, Woolrich M, Joensson M, Mohseni H, Baker A, Kringelbach ML, Deco G (2014) Exploring mechanisms of spontaneous functional connectivity in MEG: how delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations. *Neuroimage* 90:423–435.

* Caminiti R, Carducci F, Piervincenzi C, Battaglia-Mayer A, Confalone G, Visco-Comandini F, Pantano P, Innocenti GM (2013) Diameter, length, speed, and conduction delay of callosal axons in macaque monkeys and humans: comparing data from histology and magnetic resonance imaging diffusion tractography. *J Neurosci* 33:14501–14511.

* Chialvo DR (2010) Emergent complex neural dynamics. *Nature Phys* 6:744–750.

* de Cheveigné A, Simon JZ (2008) Denoising based on spatial filtering. *J Neurosci Methods* 171:331–339.

* Deco G, Jirsa V, McIntosh AR, Sporns O, Kötter R (2009) Key role of coupling, delay, and noise in resting brain fluctuations. *Proc Natl Acad Sci USA* 106:10302–10307.

* Dehaene S, Kerszberg M, Changeux JP (1998) A neuronal model of a global workspace in effortful cognitive tasks. *Proc Natl Acad Sci USA* 95:14529–14534.

* Delettre C, Messé A, Dell L-A, Foubet O, Heuer K, Larrat B, Meriaux S, Mangin J-F, Reillo I, de Juan Romero C, Borrell V, Toro R, Hilgetag CC (2019) Comparison between diffusion MRI tractography and histological tract-tracing of cortico-cortical structural connectivity in the ferret brain. *Netw Neurosci* 3:1038–1050.

* Fox MD, Snyder AZ, Vincent JL, Corbetta M, Essen DCV, Raichle ME (2005) The human brain is intrinsically organized into dynamic, anticorrelated functional networks. *Proc Natl Acad Sci USA* 102:9673–9678.

* Frischer JM, Bramow S, Dal-Bianco A, Lucchinetti CF, Rauschka H, Schmidbauer M, Laursen H, Sorensen PS, Lassmann H (2009) The
relation between inflammation and neurodegeneration in multiple sclerosis brains. Brain 132:1175–1189.

Friston KJ, Ashburner J, Frith CD, Poline J-B, Heather JD, Frackowiak RJ S (1995) Spatial registration and normalization of images. Hum Brain Mapp 3:165–189.

Ghosh A, Rho Y, McIntosh AR, Kötter R, Jirsa VK (2008) Noise during rest enables the exploration of the brain's dynamic repertoire. PLoS Comput Biol 4:e1000196.

Gong G, He Y, Concha L, Lebel C, Gross DW, Evans AC, Beaulieu C (2009) Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cereb Cortex 19:524–536.

Grech R, Cassar T, Muscat J, Camilleri KP, Fabri SG, Zervakis M, Xanthopoulos P, Sakkalis V, Vanrumste B (2008) Review on solving the inverse problem in EEG source analysis. J Neuroeng Rehabil 5:25.

Gross J, Baillet S, Barnes GR, Henson RN, Hillebrand A, Jensen O, Jerki B, Litvak V, Maess B, Oostenveld R, Parkkonen L, Taylor JR, van Wassenhove V, Wibral M, Schoffelen J-M (2013) Good practice for conducting and reporting MEG research. Neuroimage 65:349–363.

Haldeman C, Beggs JM (2005) Critical branching captures activity in living neural networks and maximizes the number of metastable states. Phys Rev Lett 94:058101.

Harris TE (1964) The theory of branching process. Vol 6. Berlin: Springer.

Hashemi M, Vattikonda AN, Sip V, Guye M, Bartolomei F, Woodman MM, Jirsa VK (2020) The Bayesian virtual epileptic patient: a probabilistic framework designed to infer the spatial map of epileptogenicity in a personalized large-scale brain model of epilepsy spread. Neuroimage 217:116839.

He Y, Chen ZJ, Evans AC (2007) Small-world anatomical networks in the human brain revealed by cortical thickness from MRI. Cereb Cortex 17:2407–2419.

Hillebrand A, Tewarie P, van Dellen E, Yu M, Carbo EWS, Douw L, Gouw AA, van Straaten ECW, Stam CJ (2016) Direction of information flow in large-scale resting-state networks is frequency-dependent. Proc Natl Acad Sci U S A 113:3867–3872.

Honey CJ, Sporns O, Cammoun L, Gigandet X, Thiran JP, Meuli R, Hagmann P (2009) Predicting human resting-state functional connectivity from structural connectivity. Proc Natl Acad Sci U S A 106:2035–2040.

Horowitz A, Barazany D, Tavor I, Bernstein M, Yovel G, Assaf Y (2015) In vivo correlation between axon diameter and conduction velocity in the human brain. Brain Struct Funct 220:1777–1788.

Klein A, Tourville J (2012) Labeled brain images and a consistent human cortical labeling protocol. Front Neurosci 6:171.

Lardone A, Liparoti M, Sorrentino P, Rucco R, Jacini F, Polverino A, Minino R, Pesoli M, Baselice F, Sorriso A, Ferraioli G, Sorrentino G, Mandolesi L (2018) Mindfulness meditation is related to long-lasting changes in hippocampal functional topology during resting state: a magnetoencephalography study. Neural Plast 2018:5340717.

Leemans A, Jones DK (2009) The B-matrix must be rotated when correcting for subject motion in DTI data. Magn Reson Med 61:1336–1349.

Nichols TE, Holmes AP (2002) Nonparametric permutation tests for functional neuroimaging: a primer with examples. Hum Brain Mapp 15:1–25.

Niebur E, Schuster HG, Kammen DM (1991) Collective frequencies and metastability in networks of limit-cycle oscillators with time delay. Phys Rev Lett 67:2753–2756.

Nolte G (2003) The magnetic lead field theorem in the quasi-static approximation and its use for magnetoencephalography forward calculation in realistic volume conductors. Phys Med Biol 48:3637–3652.

Oostenveld R, Fries P, Maris E, Schoffelen J-M, Oostenveld R, Fries P, Maris E, Schoffelen J-M (2011) FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput Intell Neurosci 2011:156869.

Petkoski S, Jirsa VK (2019) Transmission time delays organize the brain network synchronization. Philos Trans A Math Phys Eng Sci 377:20180132.

Petkoski S, Jirsa VK (2022) Normalizing the brain connectome for communication through synchronization. Network Neuroscience, pp 1–56.

Petkoski S, Palva JM, Jirsa VK (2018) Phase-lags in large scale brain synchronization: methodological considerations and in-silico analysis. PLoS Comput Biol 14:e1006160.

Purves D, Augustine GJ, Fitzpatrick D, Katz LC, LaMantia AS, McNamara JO, Williams SM (2001) Increased conduction velocity as a result of myelination. In: Neurosciece, Ed 2. Sunderland: Sinauer Associates.

Reuter F, Del Cul A, Malikova I, Naccache L, Confort-Gouny S, Cohen L, Cherif AA, Cozzone PJ, Pelletier J, Ranjeva JP, Dehaene S, Audoin B (2009) White matter damage impairs access to consciousness in multiple sclerosis. Neuroimage 44:590–599.

Rosanova M, Gosseries O, Casarotto S, Boly M, Casali AG, Bruno MA, Mariotti M, Boveroux P, Tononi G, Laureys S, Massimini M (2012) Recovery of cortical effective connectivity and recovery of consciousness in vegetative patients. Brain 135:1308–1320.

Rucco R, Bernardo P, Lardone A, Baselice F, Pesoli M, Polverino A, Bravaccio C, Granata C, Mandolesi L, Sorrentino G, Sorrentino P (2020) Neuronal avalanches to study the coordination of large-scale brain activity: application to Rett syndrome. Front Psychol 11:550749.

Sadasivan PK, Narayana Dutt D (1996) SVD based technique for noise reduction in electroencephalographic signals. Signal Process 55:179–189.

Sanz-Leon P, Knock SA, Spiegler A, Jirsa VK (2015) Mathematical framework for large-scale brain network modeling in the virtual brain. Neuroimage 111:385–430.

Schriner M, McIntosh AR, Jirsa V, Deco G, Ritter P (2018) Inferring multi-scale neural mechanisms with brain network modelling. Elife 7:e28927.

Schmidt P, Gaser C, Arsic M, Buck D, Förschler A, Berthele A, Hoshi M, Ilg R, Schmid VJ, Zimmer C, Hemmer B, Mühlau M (2012) An automated tool for detection of FLAIR-hyperintense white-matter lesions in multiple sclerosis. Neuroimage 59:3774–3783.

Shen K, Mišić B, Cipollini BN, Bezgin G, Buschkuehl M, Hutchison RM, Jaeggi SM, Kross E, Peltier SJ, Everling S, Jonides J, McIntosh AR, Berman MG (2015) Stable long-range interhemispheric coordination is supported by direct anatomical projections. Proc Natl Acad Sci U S A 112:6473–6478.

Shriki O, Alstott J, Carver F, Holroyd T, Henson RNA, Smith ML, Coppola R, Bullmore E, Plenz D (2013) Neuronal avalanches in the resting MEG of the human brain. J Neurosci 33:7079–7090.

Smith SM (2002) Fast robust automated brain extraction. Hum Brain Mapp 17:143–155.

Smith SM, Jenkinson M, Woolrich MW, Beckmann CF, Behrens TEJ, Johansen-Berg H, Bannister PR, De Luca M, Drobnjak I, Flitney DE, Niazy RK, Saunders J, Vickers J, Zhang Y, De Stefano N, Brady JM, Matthews PM (2004) Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage 23 [Suppl 1]:S208–S219.

Sorrentino P, Rucco R, Jacini F, Trojsi F, Lardone A, Baselice F, Femiano C, Santangelo G, Granata C, Vettoliere A, Monsurrò MR, Tedeschi G, Sorrentino G (2018) Brain functional networks become more connected as amyotrophic lateral sclerosis progresses: a source level magnetoencephalographic study. Neuroimage Clin 20:564–571.

Sorrentino P, Rucco R, Baselice F, De Micco R, Tessitore A, Hillebrand A, Mandolesi L, Breakspear M, Gollo LL, Sorrentino G (2021a) Flexible brain dynamics underpins complex behaviours as observed in Parkinson's disease. Sci Rep 11:4051.

Sorrentino P, Seguin C, Rucco R, Liparoti M, Lopez ET, Bonavita S, Quarantelli M, Sorrentino G, Jirsa V, Zalesky A (2021b) The structural connectome constrains fast brain dynamics. bioRxiv 393017. doi: 10.1101/ 2020.11.25.393017.

Sporns O, Zwi JD (2004) The small world of the cerebral cortex. Neuroinformatics 2:145–162.

Sporns O, Tononi G, Kötter R (2005) The human connectome: a structural description of the human brain. PLoS Comput Biol 1:e42.

Stadelmann C, Wegner C, Brück W (2011) Inflammation, demyelination, and degeneration. Recent insights from MS pathology. Biochim Biophys Acta 1812:275–282.

Stam CJ (2014) Modern network science of neurological disorders. Nat Rev Neurosci 15:683–695.

Tewarie P, Abeysuriya R, Byrne Á, O'Neill GC, Sotiropoulos SN, Brookes MJ, Coombes S (2019) How do spatially distinct frequency specific MEG networks emerge from one underlying structural connectome? The role of the structural eigenmodes. Neuroimage 186:211–220.

Thompson AJ, Banwell BL, Barkhof F, Carroll WM, Coetzee T, Comi G, Correale J, Fazekas F, Filippi M, Freedman MS, Fujihara K, Galetta SL, Hartung HP, Kappos L, Lublin FD, Marrie RA, Miller AE, Miller DH,
Montalban X, Mowry EM, et al. (2018a) Diagnosis of multiple sclerosis: 2017 revisions of the McDonald criteria. Lancet Neurol 17:162–173.

Thompson AJ, Baranzini SE, Geurts J, Hemmer B, Ciccarelli O (2018b) Multiple sclerosis. The Lancet 391:1622–1636.

Tononi G, Sporns O (2003) Measuring information integration. BMC Neurosci 4:31–31.

Tzourio-Mazoyer N, Landeau B, Papathanassiou D, Crivello F, Etard O, Delcroix N, Mazoyer B, Joliot M (2002) Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage 15:273–289.

Watts D, Strogatz S (1998) Collective dynamics of ‘small-world’ networks. Nature 393:440–442.

Waxman SG (1980) Determinants of conduction velocity in myelinated nerve fibers. Muscle Nerve 3:141–150.

Waxman SG, Bennett MVL (1972) Relative conduction velocities of small myelinated and non-myelinated fibres in the central nervous system. Nat New Biol 238:217–219.

Zalesky A, Fornito A, Cocchi L, Gollo LL, Breakspear M (2014) Time-resolved resting-state brain networks. Proc Natl Acad Sci U S A 111:10341–10346.

Zamani Esfahlani F, Jo Y, Faskowitz J, Byrge L, Kennedy DP, Sporns O, Betzel RF (2020) High-amplitude cofluctuations in cortical activity drive functional connectivity. Proc Natl Acad Sci U S A 117:28393–28401.