## Abstract

Brain activity during rest displays complex, rapidly evolving patterns in space and time. Structural connections comprising the human connectome are hypothesized to impose constraints on the dynamics of this activity. Here, we use magnetoencephalography (MEG) to quantify the extent to which fast neural dynamics in the human brain are constrained by structural connections inferred from diffusion MRI tractography. We characterize the spatio-temporal unfolding of whole-brain activity at the millisecond scale from source-reconstructed MEG data, estimating the probability that any two brain regions will significantly deviate from baseline activity in consecutive time epochs. We find that the structural connectome profoundly shapes rapid spreading of neuronal avalanches, evidenced by a significant association between these transition probabilities and structural connectivity strengths (r=0.37, p<0.0001). This finding opens new avenues to study the relationship between brain structure and neural dynamics.

## Introduction

The structural scaffolding of the human connectome^{1} constrains the unfolding of large-scale coordinated neural activity towards a restricted *functional repertoire*^{2}. While functional magnetic resonance imaging (fMRI) can elucidate this phenomenon at relatively slow timescales^{3–5}, brain activity shows rich dynamic behaviour across multiple timescales, with faster activity nested within slower ones. Here, we exploit the high temporal resolution of resting-state magnetoencephalography (MEG) data to study the spatial spread of perturbations of local activations (“neuronal avalanches”) in healthy adults, aiming to establish whether the structural connectome constrains the spread of avalanches among regions^{6,7}. We find that avalanche spread is significantly more likely between pairs of grey matter regions that are structurally connected, as inferred from diffusion MRI tractography. This result provides cross-modal empirical evidence suggesting that connectome topology constrains fast-scale transmission of neural information, linking brain structure to brain dynamics.

## Results

Structural connectomes were mapped for 58 healthy adults (26 females, mean age ± SD: 30.72 ± 11.58) using diffusion MRI tractography and regions defined based on the Automated Anatomical Labeling (AAL) and the Desikan-Killiany-Tourville (DKT) atlases. Interregional streamline counts derived from whole-brain deterministic tractography quantified the strength of structural connectivity between pairs of regions (see SI extended methods). Group-level connectomes were computed by averaging connectivity matrices across participants.

MEG signals were pre-processed and source reconstructed for both the AAL and DKT atlases. All analyses were conducted on source-reconstructed signal amplitudes. Each signal amplitude was z-scored and binarized such that, at any time point, a z-score exceeding a given threshold was set to 1 (active); all other timepoints were set to 0 (inactive). An avalanche was defined as starting when any region exceeded this threshold, and finished when no region was active. An avalanche-specific transition matrix (TM) was calculated, where element (*i, j*) represented the probability that region *j* was active at time *t+δ*, given that region *i* was active at time *t*, where *δ*∼3ms. The TMs were averaged per participant, and then per group, and finally symmetrized.

We found striking evidence of an association between avalanche transition probabilities and structural connectivity strengths (Fig. 2), suggesting that regional propagation of fast-scale neural avalanches is partly shaped by the axonal fibers forming the structural connectome (r=0.37, p<0.0001). Specifically, the association was evident for different activation thresholds and both the AAL and DKT connectomes (AAL atlas: for threshold z=2.5, r=0.38; for threshold z=3.0, r=0.37; for threshold z=3.5, r=0.34; DKT atlas: for threshold z=2.5, r=0.32; for threshold z=3.0, r=0.31; for threshold z=3.5, r=0.30; in all cases, p <0.0001), as well as for individual- and group-level connectomes, although associations were stronger for group-level analyses (see Fig. 2, panel A). Associations were also evident within specific frequency bands and for alternative methods to estimate the TM (see SI).

Next, we sought to test whether the associations were weaker for randomized transition matrices computed after randomizing the times of each avalanche while keeping the spatial structure unchanged (see SI extended methods). Randomized transition matrices resulted in markedly weaker associations with structural connectivity, compared to the actual transition matrices (AAL atlas, z-score=3: mean r = 0.25, observed r= 0.38, p<0.001). This suggests that the empirical organization of the connectome significantly shapes the temporally resolved propagation of neural activity. We replicated these findings for a group-level connectome derived using structural data from 200 healthy adults in the Human Connectome Project (r=0.11, p<0.001, z-score=3; Methods). Our results were thus robust to multiple connectome mapping pipelines and parcellation atlases, significant for both group-averaged and individual connectomes, and could not be explained by chance transitions.

## Discussion

Our results provide new insight into the propagation of fast-evolving brain activity in the human connectome. We show that the spatial unfolding of neural dynamics at the millisecond scale is shaped by the network of large-scale axonal projections comprising the connectome, thereby constraining exploration of the brain’s putative functional repertoire. While previous studies provide evidence of coupling between structural connectivity and functional MRI activity ^{3,8,9}, the neural signals measured with MEG in the present study are orders of magnitude faster, enabling investigation of intrinsic neural dynamics nested in slow activity^{10}. Our findings suggest that long-term structure-function coupling previously uncovered with functional MRI occurs against a backdrop of faster fluctuations, which are also constrained by the connectome and may enable individuals to rapidly respond to changing environments and new cognitive demands ^{11}. Finally, our results explain how the large-scale activity unfolding in time might lead to the previous observation that average resting-state functional connectivity has topological features that mirror those of the structural connectome^{12}. The proposed framework links the large-scale spreading of aperiodic, locally generated perturbations to the structural connectome, and might be further exploited to investigate polysynaptic models of network communication, which aim to describe patterns of signalling between anatomically unconnected regions^{13,14}. Therefore, our work provides a foundational step towards elucidating the mechanisms governing communication in the human connectome. In turn, this can be exploited to predict the effects of structural lesions on behaviour and/or clinical phenotypes, under the above-mentioned hypothesis that structure influences behavioural outcomes by constraining global dynamics. In conclusion, using MEG to study neuronal avalanches, we provide a new framework to link fast neural dynamics to the structural connectome.

## Methods

### MEG pre-processing

MEG pre-processing and source reconstruction were performed as in^{15}. In short, the MEG registration was divided in two eyes-closed segments of 3:30 minutes each. To identify the position of the head, four anatomical points and four position coils were digitized. Electrocardiogram (ECG) and electro-oculogram (EOG) signals were also recorded. The MEG signals, after an anti-aliasing filter, were acquired at 1024 Hz, then a fourth order Butterworth IIR band-pass filter in the 0.5-48 Hz band was applied. To remove environmental noise, measured by reference magnetometers, we used Principal Component Analysis. We adopted Independent Component Analysis to clean the data from physiological artifacts, such as eye blinking (if present) and heart activity (generally one component). Noisy channels were identified and removed manually by an expert rater. 47 subjects were selected for further analysis. The time series of neuronal activity were reconstructed based on the Automated Anatomical Labeling (AAL) and the Desikan-Killiany-Tourreville (DKT) atlases. To do this, we used the Linearly Constrained Minimum Variance (LCMV) beamformer algorithm based on the native MRIs. Finally, we excluded the ROIs corresponding to the cerebellum because of their low reliability in MEG. However, when these regions were included, the results were replicated. All the preprocessing and the source reconstruction were performed using the Fieldtrip toolbox.

### Transition matrices

Each source reconstructed signal was binned (such as to obtain a branching ratio ∼1, see SI) and then z-scored and binarized, such that, at any time bin, a z-score exceeding ± 3 was set to 1 (active); all other time bins were set to 0 (inactive). See SI for further details. Alternative z-score thresholds (i.e.2.5 and 3.5) were tested. An avalanche was defined as starting when any region is above threshold, and finishing when no region is active, as in^{15}). Avalanches shorter than 10 time bins (∼30 msec) were excluded. However, the analyses were repeated including only avalanches longer than 30 time bins (∼90 msec), to focus on rarer events that are highly unlikely to be noise, and including all avalanches, and the results were unchanged. The results reported refer to a branching ratio of 1. However, they were unchanged for binning ranging from 1 to 5. An avalanche-specific transition matrix (TM) was calculated, where element (*i, j*) represented the probability that region *j* was active at time *t+δ*, given that region *i* was active at time *t*, where *δ*∼3ms. The TMs were averaged per participant, and then per group, and finally symmetrized. The introduction of a time-lag makes it unlikely that our results can be explained trivially by volume conduction (i.e. the fact that multiple sources are detected simultaneously by multiple sensors, generating spurious zero-lags correlations in the recorded signals). For instance, for a binning of 3, as the avalanches proceed in time, the successive regions that are recruited do so after roughly 3 msecs (and 5 msecs for the binning of 5). Activations occurring simultaneously do not contribute to the estimate of the transition matrix. To further rule out this possibility, we have computed the transition matrices in a different way, to include even longer delays (see SI for details). In short, after the initial time-bin of an avalanche, we keep track of what other regions were recruited after the first perturbation. Importantly, we did not scroll through the avalanche in time, as previously described, so as to include time delays as long as the avalanche itself. The results were confirmed with this procedure as well (see SI for details). Finally, we explored transition matrices estimated using frequency-specific signals (see SI). In short, we filtered the source-reconstructed signal in the classical frequency bands (delta, 0.5 – 4 Hz; theta 4 – 8 Hz; alpha 8 – 13 Hz; beta 13 – 30 Hz; gamma 30 – 48 Hz), before computing neuronal avalanches and the transition matrix. The results show that the structural connectome mediates the propagation of perturbations in all frequency bands (delta, r=0.39; theta, r=0.29; alpha, r=0.32; beta, r=0.32; gamma, r=0.32; p<0.0001 in all bands). Further studies are warranted to investigate frequency-specific behaviour of large-scale avalanche spreading.

### Diffusion MRI pre-processing and structural connectome mapping

Diffusion MRI data were acquired for the same individuals using a 1.5 Tesla machine (Signa, GE Healthcare). Preprocessing was performed using the software modules provided in the FMRIB Software Library (FSL, http://fsl.fmrib.ox.ac.uk/fsl). Data were corrected for head movements and eddy current distortions using the “eddy_correct” routine, rotating diffusion sensitizing gradient directions accordingly, and a brain mask was obtained from the B0 images using the Brain Extraction Tool routine. A diffusion-tensor model was fitted at each voxel, and fiber tracks were generated over the whole brain using deterministic tractography, as implemented in Diffusion Toolkit (FACT propagation algorithm, angle threshold 45°, spline-filtered, masking by the FA maps thresholded at 0.2). For tractographic analysis, the ROIs of the AAL atlas and of a MNI space-defined volumetric version of the Desikan-Killiany-Tourville (DKT) ROI atlas were used, both masked by the GM tissue probability map available in SPM (thresholded at 0.2). To this end, for each subject, FA volumes were normalized to the MNI space using the FA template provided by FSL, using the spatial normalization routine available in SPM12, and the resulting normalization matrices were inverted and applied to the ROIs, to apply them onto each individual. The quality of the normalization was assessed visually. For each individual, the number of streamlines interconnecting each pair of regions was enumerated using custom software written in Interactive Data Language (IDL, Harris Geospatial Solutions, Inc., Broomfield, CO, USA). Results were replicated using both the AAL and the DKT atlases. In supplementary analyses, connectomes were also mapped using diffusion MRI data for 200 participants from the Human Connectome Project using an alternative workflow. The resulting individual connectomes were then averaged to yield a group-consensus connectome. Further details are available in SI. See Fig.1 for an overview on the methods.

### Statistical analysis

The Spearman rank correlation coefficient was used to assess the association between transition probabilities and structural connectivity. A correlation coefficient was computed separately for each individual across all pairs of regions. Transition matrices were symmetrized before this computation.

Randomized transition matrices were generated to ensure that associations between transition probabilities and structural connectivity could not be attributed to chance. Avalanches were randomized across time, without changing the order of avalanches at each time step. We generated a total of 1000 randomized transition matrices and the Spearman rank correlation coefficient was computed between each randomized matrix and structural connectivity. This yielded a distribution of correlation coefficients under randomization. The proportion of correlation coefficients that were greater than, or equal to, the observed correlation coefficient provided a p-value for the null hypothesis that structure-function coupling was attributable to random transition events.

## Footnotes

↵* Co-senior authors

All sections of the paper have been modified and implemented to facilitate a better understanding of the study.