The structural connectome constrains fast brain dynamics
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 spatiotemporal unfolding of wholebrain activity at the millisecond scale from sourcereconstructed 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 relates to, and likely affects, the 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 (Sporns et al., 2005) constrains the unfolding of largescale coordinated neural activity towards a restricted functional repertoire (Deco et al., 2011). While functional magnetic resonance imaging (fMRI) can elucidate this phenomenon at relatively slow timescales (Honey et al., 2007; Goñi et al., 2014; Zalesky et al., 2014), brain activity shows rich dynamic behaviour across multiple time scales, with faster activity nested within slower scales. Here, in healthy young adults, we exploit the high temporal resolution of restingstate magnetoencephalography (MEG) data to study the spatial spread of perturbations of local activations representative of neuronal avalanches. We aim to establish whether the structural connectome constrains the spread of avalanches among regions (Beggs and Plenz, 2004; Shriki et al., 2013). 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 crossmodal empirical evidence suggesting that connectome topology constrains fastscale 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 wholebrain deterministic tractography quantified the strength of structural connectivity between pairs of regions. Streamline counts were normalized by regional volume. Grouplevel connectomes were computed by averaging connectivity matrices across participants.
MEG signals were preprocessed and source reconstructed for both the AAL and DKT atlases. All analyses were conducted on sourcereconstructed signal amplitudes. Each signal amplitude was zscored and binarized such that, at any time point, a zscore exceeding a given threshold was set to 1 (active); all other time points 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 avalanchespecific 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 δ~3 ms. The TMs were averaged per participant, and then per group, and finally symmetrized. Figure 1 provides an overview of the pipeline.
We found striking evidence of an association between avalanche transition probabilities and structural connectivity strengths (Figure 2), suggesting that regional propagation of fastscale neural avalanches is partly shaped by the axonal fibres forming the structural connectome (r = 0.40, 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.41; for threshold z = 3.0, r = 0.40; for threshold z = 3.5, r = 0.39; DKT atlas: for threshold z = 2.5, r = 0.38; for threshold z = 3.0, r = 0.37; for threshold z = 3.5, r = 0.35; in all cases, p<0.0001), as well as for individual and grouplevel connectomes, although associations were stronger for grouplevel analyses (see Figure 2A).
We also investigated this phenomenon within specific frequency bands. Associations were evident in all the classical frequency bands: delta (0.5–4 Hz; r = 0.39), theta (4–8 Hz; r = 0.29), alpha (8–13 Hz; r = 0.32), beta (13–30 Hz; r = 0.32), and gamma (30–48 Hz; r = 0.32), with p<0.0001 for all bands (see Supplementary file 1). Supplementary analyses suggested that these results could not be attributable to volume conduction confounds (see section Field spread analysis).
Next, we sought to test whether the associations were weaker for randomized TMs computed after randomizing the times of each avalanche while keeping the spatial structure unchanged. Randomized TMs resulted in markedly weaker associations with structural connectivity compared to the actual TMs (AAL atlas, zscore = 3: mean r = 0.26, observed r = 0.40, p<0.001). Note that the mean correlation coefficient was greater than zeros for the randomized data because the randomization process preserved basic spatial attributes in the data. We also found that the findings remained significant after excluding subcortical regions (with lower signaltonoise ratios). Finally, we replicated these findings for a grouplevel connectome derived using diffusion MRI acquired from 200 healthy adults in the Human Connectome Project (r = 0.11, p<0.001, zscore = 3; see Materials and methods). Our results were thus robust to multiple connectome mapping pipelines and parcellation atlases, significant for both groupaveraged and individual connectomes, and could not be explained by chance transitions and/or volume conduction effects. Collectively, these results suggest that connectome organization significantly shapes the propagation of neural activity.
Discussion
Our results provide new insight into the propagation of fastevolving brain activity in the human connectome. We show that the spatial unfolding of neural dynamics at the millisecond scale relates to the network of largescale axonal projections comprising the connectome, likely constraining the exploration of the brain’s putative functional repertoire. The short time scale of several milliseconds biases the constraint to direct connections, which is the focus of this paper. Longer delays may impose constraints upon largerscale motifs of the network and further characterize the subspaces, in which brain dynamics unfold.
Previous functional MRI studies provide evidence of coupling between structural connectivity and slow activations (Honey et al., 2007; Honey et al., 2010; Honey et al., 2009). However, intrinsic neural dynamics evolve quickly and are nested within slow activity (Saggio et al., 2017). Our findings suggest that longterm structurefunction coupling 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 (McIntosh and Jirsa, 2019).
Consistent with our findings, two recent M/EEG studies showed that functional connectivity, as estimated using amplitudeenvelope coupling (AEC), relates to structural connectivity (Glomb et al., 2020; Tewarie et al., 2019). However, in contrast to AEC, we conducted timeresolved analyses, characterizing avalanche dynamics at high temporal resolution. Further work is needed to determine the extent to which structurefunction coupling is dynamic. To this regard, our results suggest that coupling is strongest during avalanche events, consistent with established theories (Dehaene et al., 1998). Finally, our results might explain how the largescale activity unfolding in time might lead to the previous observation that average restingstate functional connectivity displays topological features that mirror those of the structural connectome (Bullmore and Sporns, 2009). Our proposed framework links the largescale 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 (Seguin et al., 2018; Seguin et al., 2019). In fact, our results show that transitions of activations are observed across regions that do not appear to be directly linked in the structural connectome. This provides evidence for polysynaptic communication.
Neuronal avalanches have been previously observed in MEG data (Shriki et al., 2013), and their statistical properties, such as a size distribution that obeys a powerlaw with an exponent of −3/2, reported. These features are compatible with those that would be predicted starting from a process operating at criticality with a branching ratio equal to 1. While beyond the scope of this paper, our framework might contribute to elucidating the role of the structural scaffolding (and its topological properties) to the emergence of the observed largescale, scalefree critical dynamics. In turn, this might be exploited to predict the effects of structural lesions on behaviour and/or clinical phenotypes.
While our findings were replicated across multiple frequency bands, structural connectivity can potentially impose frequencydependent constraints on avalanche spread. Future work should investigate frequencyspecific data to understand what leads to the emergence of avalanches and, most importantly, to the specific spatiotemporal patterns of recruited regions that defines individual (or at least groups of) avalanches in each specific frequency band.
For the present application, we reconstructed the structural connectome using a deterministic tractography algorithm. While probabilistic algorithms can provide advantages in some applications, they are prone to reconstruction of spurious connections (false positives), compared to deterministic methods, reducing connectome specificity (Sarwar et al., 2019; Zalesky et al., 2016). We used deterministic tractography because previous functional MRI studies report that structurefunctional coupling is greater for connectivity matrices inferred from deterministic tractography compared to probabilistic methods (Abeyasinghe et al., 2021). Nonetheless, additional studies are needed to clarify if and to what extent the present results are influenced by the structural connectome reconstruction method. While we replicated our findings using alternative datasets (i.e., Human Connectome Project [HCP]) and parcellations, further replication using alternative connectome mapping pipelines is warranted.
In conclusion, using MEG to study fast neuronal dynamics and diffusion MRI tractography to map connectomes, we found that the connectome significantly constrains the spatial spread of neuronal avalanches to axonal connections. Our results suggest that largescale structurefunction coupling is dynamic and peaks during avalanche events.
Materials and methods
Participants
We recruited 58 young adults (males 32/females 26, mean age ± SD was 30.72 ± 11.58) from the general community. All participants were righthanded and native Italian speakers. The inclusion criteria were (1) no major internal, neurological, or psychiatric illnesses; and (2) no use of drugs or medication that could interfere with MEG/MRI signals. The study complied with the Declaration of Helsinki and was approved by the local Ethics Committee. All participants gave written informed consent.
MRI acquisition
Request a detailed protocol3D T1weighted brain volumes were acquired at 1.5 Tesla (Signa, GE Healthcare) using a 3D MagnetizationPrepared GradientEcho BRAVO sequence (TR/TE/TI 8.2/3.1/450 ms, voxel 1 × 1 × 1 mm^{3}, 50% partition overlap, 324 sagittal slices covering the whole brain), and diffusion MRI data for individual connectome reconstruction were obtained using the following parameters: EchoPlanar Imaging, TR/TE 12,000/95.5 ms, voxel 0.94 × 0.94 × 2.5 mm^{3}, 32 diffusionsensitizing directions (5 B0 volumes). The MRI scan was performed after the MEG recording. Preprocessing of the diffusion MRI data was carried out 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 currents 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 diffusiontensor model was fitted at each voxel, and streamlines were generated over the whole brain by deterministic tractography using Diffusion Toolkit (FACT propagation algorithm, angle threshold 45°, splinefiltered, masking by the Fractional Anisotropy (FA) maps thresholded at 0.2). For tractographic analysis, the regions of interest (ROIs) of the AAL atlas and of an MNI spacedefined volumetric version of the DKT ROI atlas were used, both masked by the gray matter (GM) tissue probability map available in Statistical Parametric Mapping Software  SPM (thresholded at 0.2). To this end, for each participant, 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 subject. The quality of the normalization was assessed visually. From each subject's wholebrain tractography and corresponding GM ROI set, the number of streamlines connecting each couple of GM ROIs and the corresponding mean tract length was calculated using an inhouse software written in Interactive Data Language (IDL, Harris Geospatial Solutions, Inc, Broomfield, CO).
Connectomes in the replication dataset were constructed using an alternative mapping pipeline and diffusion MRI data from the HCP. Deterministic tractography was performed using MRtrix3 (Tournier et al., 2019) under the following parameters: FACT algorithm, 5 million streamlines, 0.5 mm propagation step size, 400 mm maximum propagation length, and 0.1 FA threshold for the termination of streamlines (Seguin et al., 2019). The number of streamlines connecting any couple of regions was normalized by the combined volume of the two regions. Structural matrices were constructed for 200 HCP participants using the AAL atlas and averaged to derive a grouplevel connectome.
MEG preprocessing
Request a detailed protocolMEG preprocessing and source reconstruction were performed as in Sorrentino et al., 2021. The MEG system was equipped with 163 magnetometers and was developed by the National Research Council of Italy at the Institute of Applied Sciences and Intelligent Systems (ISASI). All technical details regarding the MEG device are reported in Rombetto et al., 2014. In short, the MEG registration was divided into two eyesclosed segments of 3:30 min each. To identify the position of the head, four anatomical points and four position coils were digitized. Electrocardiogram (ECG) and electrooculogram (EOG) signals were also recorded. The MEG signals, after an antialiasing filter, were acquired at 1024 Hz, then a fourthorder Butterworth IIR bandpass 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 supervised independent component analysis to clean the data from physiological artefacts, such as eye blinking (if present) and heart activity (generally one component). Noisy channels were identified and removed manually by an expert rater (136 ± 4 sensors were kept). 47 subjects were selected for further analysis.
Source reconstruction
Request a detailed protocolThe time series of neuronal activity were reconstructed in 116 ROIs based on the AAL atlas (TzourioMazoyer et al., 2002; Hillebrand et al., 2016); and in 84 regions of interest based on the DKT atlas. To do this, we used the volume conduction model proposed by Nolte, 2003 applying the linearly constrained minimum variance (LCMV) beamformer algorithm (Van Veen et al., 1997) based on the native structural MRIs. Sources were reconstructed for the centroid of each ROI. Finally, we considered a total of 90 ROIs for the AAL atlas since we have excluded 26 ROIs corresponding to the cerebellum because of their low reliability in MEG (Lardone et al., 2018). All the preprocessing steps and the source reconstruction were made using the Fieldtrip toolbox (Oostenveld et al., 2011).
Neuronal avalanches and branching parameter
Request a detailed protocolTo study the dynamics of brain activity, we estimated ‘neuronal avalanches’. Firstly, the time series of each ROI was discretized calculating the zscore, then positive and negative excursions beyond a threshold were identified. The value of the threshold was set to three standard deviations (z = 3), but we tested the robustness of the results changing this threshold from 2.5 to 3.5. A neuronal avalanche begins when, in a sequence of contiguous time bins, at least one ROI is active (z >3) and ends when all ROIs are inactive (Beggs and Plenz, 2003; Shriki et al., 2013). The total number of active ROIs in an avalanche corresponds to its size.
These analyses require the time series to be binned. This is done to ensure that one is capturing critical dynamics, if present. To estimate the suitable time bin length, for each subject, each neuronal avalanches, and each time bin duration, the branching parameter σ was estimated (Haldeman and Beggs, 2005; Harris, 1964). In fact, systems operating at criticality typically display a branching ratio ~1. The branching ratio is 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 it over all the avalanches (Bak et al., 1987). More specifically:
where ${\sigma}_{i}$ is the branching parameter of the ith avalanche in the dataset, ${N}_{bin}$ is the total amount of bins in the ith avalanche, ${n}_{events}\left(j\right)$ is the total number of events active in the jth bin, and ${N}_{aval}$ is the total number of avalanche in the dataset. We tested bins from 1 to 5, and picked 3 for further analyses, given that the branching ratio was 1 for bin = 3. However, results are unchanged for other bin durations, and the branching ratio remains equal to 1 or differences were minimal (range: 0.999–1.010 – data not shown). Bins of longer duration would violate the Nyquist criterion and were thus not considered. The results shown are derived when taking into account avalanches longer than 10 time bins. However, we repeated the analysis taking into account avalanches longer than 30 time bins, as well as taking all avalanches into account, and the results were unchanged.
Transition matrices
Request a detailed protocolThe amplitude of each binned, zscored sourcereconstructed signal was binarized, such that, at any time bin, a zscore exceeding ±3 was set to 1 (active); all other time bins were set to 0 (inactive). Alternative zscore 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 Sorrentino et al., 2021. Avalanches shorter than 10 time bins (~30 ms) were excluded. However, the analyses were repeated including only avalanches longer than 30 time bins (~90 ms) to focus on rarer events (sizes of the neuronal avalanches have a fattailed distribution) that are highly unlikely to be noise, and including all avalanches, and the results were unchanged. An avalanchespecific 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 δ~3 ms. 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 lag 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 ms (and 5 ms for the binning of 5). Hence, activations occurring simultaneously do not contribute to the estimate of the TM. See below for further analyses addressing the volume conduction issue. Finally, we explored TMs estimated using frequencyspecific signals. To this end, we filtered the sourcereconstructed 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 TM, by applying a fourthorder Butterworth passband filter to the sourcereconstructed data, before proceeding to the further analysis as previously described. The results remained significant in all the explored frequency bands. This analysis was carried out for the DKT atlas, binning = 3, zscore threshold = ± 3.
Field spread analysis
Request a detailed protocolVolume conduction alone is an unlikely explanation of our results, given that simultaneous activations do not contribute to the TM, due to the time lags introduced. To confirm that volume conduction effects were negligible, the TMs were recomputed using longer delays. In short, we identified the regions that were recruited in an avalanche after the first perturbation (i.e., the initial time bin of an avalanche). Since we did not scroll through the avalanche in time, as previously described, we considered time delays as long as the avalanche itself, while minimizing the influence of short delays. This means that the avalanchespecific TM is now binary, and the ijth element is equal to 1 if region i started the avalanche (i.e., it was active at the first time bin) and region j was recruited in the avalanche at any subsequent time point, and 0 otherwise. This alternative procedure for the estimation of the TMs was carried out for the AAL atlas, in the case of binning = 3, zscore threshold = ± 3. In this case, a significant association remained between transition probabilities and structural connectivity (r = 0.36; p<0.0001). Figure 2—figure supplement 1 provides further details.
To further rule out the possibility that field spread might introduce spurious correlations that might drive the relationship between the TM and the structural connectivity matrix, we conducted further analyses involving surrogate data. We generated n white Gaussian processes, with n = 66, that is, the number of cortical regions, and we smoothed them using a zerophase polynomial filter. Then, we added 100 perturbations, where each perturbation was assigned to a randomly chosen regions and random time point, subject to the following constraints. Perturbations were separated by at least 200 samples (no overlap was allowed, i.e., the perturbations could only occur in one region at a time), their length was randomly selected among 5, 10, or 100 samples, their amplitude between 50 and 400. This procedure was carried out 47 times to obtain an independent surrogate dataset for each one of the 47 participants, which will be referred to as the ‘uncoupled’ dataset. The uncoupled dataset was then transformed using the subjectspecific leadfield matrix, yielding new surrogate sensorlevel time series, where each sensor is a weighted sum of all the sources, according to the same leadfield matrix that was used to reconstruct the real data. Noise, correlated as 1/distance among sensors, was then added to the sensorlevel time series, with an SNR = 4. Then, new sourcereconstructed time series were computed for each subject. Based on these new time series, we performed the same procedure to compute the TM as described above. Specifically, we zscored the time series, thresholded them (threshold z = ±3), retrieved the avalanchespecific TMs, averaged these within each subject and then across the group, and finally symmetrized the matrix. We then investigated the extent of the correlation between the new TM and the structural connectivity matrix. We repeated the entire procedure reported above 100 times and show that is unlikely that linear mixing alone can explain the significant association between transition probabilities and structural connectivity (p<0.001).
Statistical analysis
Request a detailed protocolThe 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. TMs were symmetrized before this computation. Randomized TMs 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 active regions at each time step. We generated a total of 1000 randomized TMs 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 pvalue for the null hypothesis that structurefunction coupling was attributable to random transition events.
Data availability
The MEG data and the reconstructed avalanches are available upon request to the corresponding author (Pierpaolo Sorrentino), 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 corresponding author will request an amendment to the local ethical committee. Conditional to approval, the data will be made available. The Matlab code is available at https://github.com/pierpaolosorrentino/TransitionMatrices (copy archived at https://archive.softwareheritage.org/swh:1:rev:3846bb80457ddd12cf8c29f36d714b5b8f64eefc).
References

Selforganized criticality: An explanation of the 1/f noisePhysical Review Letters 59:381–384.https://doi.org/10.1103/PhysRevLett.59.381

Neuronal avalanches in neocortical circuitsThe Journal of Neuroscience 23:11167–11177.https://doi.org/10.1523/JNEUROSCI.233511167.2003

Complex brain networks: graph theoretical analysis of structural and functional systemsNature Reviews Neuroscience 10:186–198.https://doi.org/10.1038/nrn2575

Emerging concepts for the dynamical organization of restingstate activity in the brainNature Reviews Neuroscience 12:43–56.https://doi.org/10.1038/nrn2961

The hidden repertoire of brain dynamics and dysfunctionNetwork Neuroscience 3:994–1008.https://doi.org/10.1162/netn_a_00107

FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological dataComputational Intelligence and Neuroscience 2011:1–9.https://doi.org/10.1155/2011/156869

FastSlow bursters in the unfolding of a high codimension singularity and the Ultraslow transitions of classesThe Journal of Mathematical Neuroscience 7:1–47.https://doi.org/10.1186/s1340801700508

Mapping connectomes with diffusion MRI: deterministic or probabilistic tractography?Magnetic Resonance in Medicine 81:1368–1384.https://doi.org/10.1002/mrm.27471

Neuronal avalanches in the resting MEG of the human brainJournal of Neuroscience 33:7079–7090.https://doi.org/10.1523/JNEUROSCI.428612.2013

The human connectome: A structural description of the human brainPLOS Computational Biology 1:e42.https://doi.org/10.1371/journal.pcbi.0010042

Localization of brain electrical activity via linearly constrained minimum variance spatial filteringIEEE Transactions on Biomedical Engineering 44:867–880.https://doi.org/10.1109/10.623056
Decision letter

Timothy E BehrensSenior Editor; University of Oxford, United Kingdom

Diego VidaurreReviewing Editor; University of Oxford, United Kingdom

Andrew J QuinnReviewer; Oxford Centre for Human Brain Activity, Wellcome Centre for Integrative Neuroimaging, United Kingdom

George O'NeillReviewer; University College London, United Kingdom
Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Acceptance summary:
This paper addresses the relationship between the electrophysiological and the anatomical connectomes, utilising a method to describe avalanches of activity. In these rapid events brain activity cascades across cortex. The current paper shows that these avalanches follow the routes of anatomical connections. The result also implies more spatial precision that most would assume possible, which makes the manuscript particularly interesting to M/EEG researchers. The reviewers therefore agree that the paper has broad interest.
Decision letter after peer review:
Thank you for submitting your article "The structural connectome constrains fast brain dynamics" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by Diego Vidaurre as a Reviewing Editor and Timothy Behrens as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Andrew J Quinn (Reviewer #1); George O'Neill (Reviewer #2).
Essential Revisions:
1. Resolve if the results are fundamentally driven by leakage
2. Determine the influence of the region size and SNR on the results.
3. Clarify the description of the method
4. Frame the work into the current literature
Reviewer #1 (Recommendations for the authors):
1) The following matlab code illustrates my concern with volume conduction (written in version R2019a). This generates a simple aperiodic signal and linearly weights it into 20 regions before computing an approximation of the avalanches analysis.
%%
rng(42);
a = poly(0.98);
x = filtfilt(1,a,randn(5000,1));
x = x(500:4500); % trim filter edges
% Add two larger pertubations
x(1000:2000) = x(1000:2000) + 400*sin(2*pi*linspace(0,0.5,1001)');
x(2500:3000) = x(2500:3000) + 1000*sin(2*pi*linspace(0,0.5,501)');
% Create + apply linear weights vector
weights = [linspace(0,1,10), linspace(1,0,10)].^3;
y = x * weights;
% Add noise, ztransform and smooth
y = y + randn(4001,1)*60;
y = zscore(y);
y = movmean(y,41,1);
% Illustrate
threshold = 1; % bit low but for illustration
figure;
subplot(211);hold on
plot(y);
plot([0 4001],[threshold threshold],'k:')
plot([0 4001],[threshold threshold],'k:')
axis('tight')
title('Linearly mixed signals')
subplot(212)
imagesc(abs(y)'>threshold)
colormap('bone')
title('Avalanches');
%%
Note that in panel 2 on the bottom the onset of the 'active' periods is temporally lagged across regions – sometimes by up to 10s or 100s of samples. The combination of amplitude weighting, noise and ztransforming can therefore induce apparent timelagged interactions. I don't know how much of the main effect in the paper is driven by this effect but it absolutely must be explored and accounted for. It is not sufficient to say that "The introduction of a timelag makes it unlikely that our results can be explained trivially by volume conduction"
2) A lot of discussion about the results uses quite forceful causal language.
line 25 "We find that the structural connectome profoundly shapes rapid spreading of neuronal avalanches"
line 100 "We show that the spatial unfolding of neural dynamics at the millisecond scale is shaped by the network of largescale axonal projections comprising the connectome, thereby constraining exploration of the brain's putative functional repertoire."
Whilst I agree that a structure>function relationship is perhaps the more likely interpretation, this is a correlational study which does not assess whether structure shapes function or vice versa – simply that there is an association. The writing would be improved by acknowledging this and adopting a more cautious interpretation.
3) Some parts of the methods writing is unclear, perhaps as the whole manuscript is so short. For example, I do not understand the second method for estimating transition probabilities outlined here:
"After the initial timebin of an avalanche, we kept 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."
4) I am unsure how to interpret the finding that the results are almost completely unchanged by bandpass filtering the data – α, β and γ even have identical Rvalues. On one hand this could indicate that the results rest completely on wide broadband effects but given the large and well established functional and topographical differences between these oscillations it seems likely that we would expect at least some difference. Particularly as similar past papers do show different structurefunction relationships between different MEG frequency bands (https://www.sciencedirect.com/science/article/abs/pii/S1053811918320603). It would be useful for the authors to discuss this point and its potential meaning in more detail.
5) It would be good to see full use of the large datasets presented here. For instance, do the results have testretest reliability across the two scans analysed per participant? Why only reproduce the finding using the HCP DWI data but not the HCP MEG data?
6) The data appears to be available from the authors 'upon request and conditional on ethics approval' but the analysis code does not appear to be available online.
Reviewer #2 (Recommendations for the authors):
I think you are definitely onto something here, but I am not quite sure if this is necessarily the full picture – yet. It's not too difficult to remedy I imagine, but I think on a technical level I'd like to query two main points on this.
I realise that this is a short communication, but it does appear to assume that no work in the M/EEG field on relating structure and function exists, any previous work been done using fMRI. I think it might be important to just give a bit more context where in the literature this might sit? Just looking around for a couple of quick examples of relating white matter tractography data to functional connectivity, there's a study relating EEG connectivity to structural connectivity (Glomb et al., 2020, Network Neurosci. 4 (3): 761787) And plenty of modelling papers which trying to look at how structure would generate observable M/EEG connectomes (eg. Tewarie et al., 2019 NeuroImage 186). I think a slightly wider look at the literature on this would be really helpful. Crucially, it's not that I don't think that this isn't novel, but rather it seems to currently seems to not acknowledge the exitance of previous work in the M/EEG field?
The use of neuronal avalanches method, is it as invariant to leakage as you initially think it is? Leakage is a product of signal variance (c.f. O'Neill et al., 2015 Phys. Med. Biol R217). So if a signal is stronger in amplitude, the leakage increases and can propagate further. So in principle, if a seed area reaches over a threshold zscore, it could continue to raise in amplitude such that a test area might also then reach a threshold? Could other (i guess more established) connectivity metrics be able to resolve this also? I was wondering even if you assume stationary variance in your signal, you might be able to analytically estimate 'leakage' between two areas (see the O'neill paper above for an analytic expression). Could you then (potentially) have a leakage connectome and see how well that predicts your avalanches perhaps? If its considerably less than the structural, I'd be a lot more confident in what I am seeing.
Reviewer #3 (Recommendations for the authors):
Is the inhouse software written in Interactive Data Language (used for calculating the number of streamlines connecting each pair of ROIs and the corresponding mean tract length) publicly available? If so, please provide a web link to it.
A bit of extra information about the MEG scanner would also be useful, in addition to reporting the actual number of channels that were manually removed.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "The structural connectome constrains fast brain dynamics" for further consideration by eLife. Your revised article has been evaluated by Timothy Behrens (Senior Editor) and Diego Vidaurre (Reviewing Editor).
The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:
As you can see, all Reviewers are satisfied with the responses, but Reviewer #1 requests some further clarifications on the simulations, and a possible improvement. Since this seems very relevant to the contents of the manuscript, we would like to ask you to go through one more revision.
Reviewer #1 (Recommendations for the authors):
Many thanks to the authors for a clear and thorough response to my queries. Overall, the revised manuscript is greatly improved and the changes largely address my concerns.
I must ask for one final clarification, and possible addition, to the simulation scheme looking at leakage. The revision doesn't state whether any noise is added at the sensor level between projection through the leadfields and subsequent beamforming – is any noise added here? If so, please in include any details (and my apologies if I missed this detail somewhere)
If not, then this additional step should be added to the pipeline. As it stands, the simulation without sensor noise creates unusually favourable conditions for the beamformer. Environmental and sensor noise have their own dynamics and can be correlated in time and across sensors, these correlations mask neuronal signals and can be tricky for the beamformer to completely remove.
A gold standard would be to add sensor noise from an empty room recording, or if this is impractical then some modestly scaled noise whose sensor x sensor covariance matrix an example empty room recording.
The correlations of both the 'coupled' and 'uncoupled' simulations are substantially lower than the observed value in the 'correlation of lineally mixed surrogates' figure. I suspect this is partially due to favourable beamforming from missing sensor noise, as it stands the beamformer is able to remove nearly all spurious events. Including the sensor noise would likely induce some more spurious connections in both the coupled and uncoupled cases, it would be a very strong validation of the method to demonstrate that these are still distinguishable from the observed effects.
Reviewer #2 (Recommendations for the authors):
The authors have met my initial concerns with the submission. Assuming other reviewers feel the same I am happy with this being published.
Reviewer #3 (Recommendations for the authors):
Thank you very much for your revisions. My comments have been addressed satisfactorily.
https://doi.org/10.7554/eLife.67400.sa1Author response
Essential Revisions:
1. Resolve if the results are fundamentally driven by leakage
2. Determine the influence of the region size and SNR on the results.
3. Clarify the description of the method
4. Frame the work into the current literature
To explore the effects of leakage, we have used a model very similar to the one proposed by reviewer one, and used the leadfield matrices that were used in our study to linearly mix the insilico data. We explore how this independent but linearlymixed time series relate to structural connectivity, and we provide evidence that linear mixing alone is unlikely to explain our results. These analyses are shown in the response to reviewer one.
To address the effect of region size, we have normalized the streamline counts by the volume of the regions and recomputed the correlation between the MEG derived transition matrices and the structural connectivity matrices. We find that for both the DKT and the AAL parcellations, the correlations improve after this normalization step. The improvement is more marked for the DKT atlas, which now reaches levels of correlation similar to the ones obtained by the AAL atlas. The analyses reported in the revised manuscript refer to these updated findings, and previous results for the correlations with the nonnormalized streamline counts are no longer reported. These analyses are reported in the response to reviewer three.
To address the issue of the potential effect of SNR, we have repeated the analysis including only cortical regions for both the DKT and the AAL. This is done under the hypothesis that, by eliminating deep sources, one takes into account regions that have a comparable SNRs. Further analyses addressing this point are reported in the response to reviewer three.
Finally, we have now framed our manuscript within the current literature, trying to highlight similarities and differences in the interpretation of the data as compared to other similar studies. These changes are mainly reported in the response to reviewer two. In what follows, one can find a detailed, pointtopoint response to the reviewers.
Reviewer #1 (Recommendations for the authors):
1) The following matlab code illustrates my concern with volume conduction (written in version R2019a). This generates a simple aperiodic signal and linearly weights it into 20 regions before computing an approximation of the avalanches analysis.
%%
rng(42);
a = poly(0.98);
x = filtfilt(1,a,randn(5000,1));
x = x(500:4500); % trim filter edges
% Add two larger pertubations
x(1000:2000) = x(1000:2000) + 400*sin(2*pi*linspace(0,0.5,1001)');
x(2500:3000) = x(2500:3000) + 1000*sin(2*pi*linspace(0,0.5,501)');
% Create + apply linear weights vector
weights = [linspace(0,1,10), linspace(1,0,10)].^3;
y = x * weights;
% Add noise, ztransform and smooth
y = y + randn(4001,1)*60;
y = zscore(y);
y = movmean(y,41,1);
% Illustrate
threshold = 1; % bit low but for illustration
figure;
subplot(211);hold on
plot(y);
plot([0 4001],[threshold threshold],'k:')
plot([0 4001],[threshold threshold],'k:')
axis('tight')
title('Linearly mixed signals')
subplot(212)
imagesc(abs(y)'>threshold)
colormap('bone')
title('Avalanches');
%%
Note that in panel 2 on the bottom the onset of the 'active' periods is temporally lagged across regions – sometimes by up to 10s or 100s of samples. The combination of amplitude weighting, noise and ztransforming can therefore induce apparent timelagged interactions. I don't know how much of the main effect in the paper is driven by this effect but it absolutely must be explored and accounted for. It is not sufficient to say that "The introduction of a timelag makes it unlikely that our results can be explained trivially by volume conduction"
We agree on the importance of this issue. We explored this point to check if spatial leakage could explain the correlation we found between the transition and structural connectivity matrices.
Firstly, on the same lines suggested by referee one, we simulated n simple aperiodic signals, with n = number of regions (sources) = 66, since we focused on cortical regions. Consistent with the referee’s suggestions, we generated n white Gaussian processes and introduced a temporal correlation in each, using a zerophase, allpole polynomial filter. Then, for 100 distinct times, we added a large perturbation to randomly selected regions, with a duration selected randomly among 5, 10 or 100 samples, and with an amplitude randomly selected between 50 and 400, as to explore the effects of perturbations of different intensities. We ensured that at least 200 samples separate any two subsequent perturbations in time. In other words, two perturbations cannot overlap in this surrogate data. However, the same region might have more than one perturbation over time. As such, we refer this surrogate data as “uncoupled”. We repeated this procedure for each subject, such that each of the 47 participants had a set of “uncoupled”, time series, each containing 100, nonoverlapping random perturbations generated at random locations. The “uncoupled” dataset represents local perturbations that occur at random without any influence from brain structure.
Based on the “uncoupled” surrogate data, we simulated, for each subject, timeseries coupled by the subjectspecific structural connections inferred from tractography. To do this, we modified the uncoupled timeseries as follows: the activity of each region i at time t was given by the summation of the activity of every other region j, at time t – $\delta $, weighted by the normalized streamline count between regions i and , as in :
where i_{coupled}(t) is the value of the new coupled time series at time t, i_{uncoupled}(t) is the corresponding value of the uncoupled i^{th} series at time t, N is the number of regions, j(tδ) is the value of the j^{th} uncoupled time series at time t – δ, with δ = 15 samples, W_{ij} is the streamline count of the ij^{th} structural edge.
We named this set of surrogate data “coupled”. While this only phenomenologically simulates a firstorder dependence with homogeneous delays, it suffices to make a point about field spread, and it is not intended (of course) to simulate any mechanism.
We applied the same statistical procedures described in our original manuscript to the surrogate data (uncoupled and coupled). In particular, we zscored the surrogate data, applied a threshold (z>±3), computed neuronal avalanches and the corresponding transition matrix (including all avalanche sizes). Finally, as per our original analysis, we averaged transition probabilities across individuals and then symmetrized the obtained transition matrix. In Author response image 1, one can see the resulting transition matrices for the uncoupled and coupled surrogate data, as well as the structural connectivity matrix to the right.
We correlated, using Spearman’s correlation, each of the transition matrices (coupled and uncoupled) with the structural connectivity matrix. As expected, in the case of the uncoupled regions, the corresponding TM does not relate to structural connectivity. However, the TM corresponding to the coupled time series, showed a direct correlation with the dti (r=0.33, p<0.0001).Based on these surrogate models, we went on to study the effect of leakage. To this end, we considered the subjectspecific leadfield matrices, that had been used to beamform the original data. Hence, each one of the 47 subjectspecific leadfields (i.e. one per each of the 47 participants) was used to project the surrogate data to the sensor space, obtaining 47 sets of timeseries in the sensor space, derived from the uncoupled data, and 47 sensor space time series based on the coupled data. Note that the leadfield estimates the signal at the location of each sensor as a linear combination of all the sources. Finally, we estimated the beamformer weights from this newly generated sensorspace timeseries, and reconstructed new sourcelevel timeseries. This approach enables the application of the mixing matrix (leadfield) that was actually used per each subject to compute the data obtained in the main results. As such, so we believe this allows us to be as adherent as possible to the procedure that was actually used to estimate the main results presented in the manuscript.
Using these newly reconstructed sourcelevel timeseries we applied the same statistical procedures to map transition matrices for the surrogate datasets, , which are shown, also in log scale, in Author response image 2.
Similarly as before, we correlated the transition matrix with the structural connectivity matrix using Spearman’s correlation. As expected, one can see that a correlation exists between the coupled surrogate and the structural matrix (r=0.12, p=0.001), but not for the uncoupled surrogate (r=0.02, p=0.59). The absence of an association between transition and structural connectivity in the uncoupled surrogate data suggests that the potential effects of leakage are an unlikely explanation for the associations reported in the real data. Importantly, the relationship that we observed in real data is stronger than the one we obtained from the simulated coupled data. Multiple reasons could explain this fact, including that our coupled model is limited to firstorder dependencies in structure, while multiple lines of evidence show that white matter bundles also induce dependencies of higher order (Seguin et al., 2020).The entire procedure described above was carried out one hundred times, whereby a new random sequence (generated as described above) per each participant was generated, multiplied by the subjectspecific leadfield matrix, sourcereconstructed, then the subjectspecific transition matrices were computed, and averaged per group. The distribution of the retrieved 100 correlations is reported in Author response image 3 for the case of coupled surrogates:
It is easy to see that the randomly generated sequences of uncoupled perturbations, which are then mixed linearly exactly as they real data were mixed, do not reproduce high correlations with structural connectivity, hence making it unlikely that linear mixing alone is a plausible explanation of the relationship that we found.This analysis has been partly reported as a new section in the methods, as follows:
“Field spread analysis
Volume conduction alone is an unlikely explanation of our results, given that simultaneous activations do not contribute to the transition matrix, due to the time lags introduced. […] We repeated the entire procedure reported above one hundred times, and show that is unlikely that linear mixing alone can explain the significant association between transition probabilities and structural connectivity (p<0.001).”
Finally, a small note on the code that reviewer 1 sent us. Firstly, we wish to acknowledge the high quality of the review, and the significant efforts of the reviewer. We are sincerely grateful. We suggest that a small modification to the code would better align this toy example with the procedure used in the manuscript. Namely, we would change y = x * weights to y =(x + randn(4001,1)*60) * weights; and eliminate the following line: y = y + randn(4001,1)*60. This would better conform to our case since noise also should be scaled by the weights. When running the code with such modifications, the output figure becomes – see Author response image 4:
In other words, the zscore has now normalized all signals equally, as in the procedure we used. The constant turquoise line at y=0 is due to one weight that was set to zero. We understand that these are toy models, and we are only hoping to illustrate how subtle changes to the model can lead to significant changes. However, the shape of the avalanches and some of the statistical properties that appear phenomenologically could not be explained easily by such a toy model (in our opinion, shape collapses are particularly convincing to this regard – please refer to the answer to referee three as well for more considerations and analyses on this point).
2) A lot of discussion about the results uses quite forceful causal language.
line 25 "We find that the structural connectome profoundly shapes rapid spreading of neuronal avalanches"
line 100 "We show that the spatial unfolding of neural dynamics at the millisecond scale is shaped by the network of largescale axonal projections comprising the connectome, thereby constraining exploration of the brain's putative functional repertoire."
Whilst I agree that a structure>function relationship is perhaps the more likely interpretation, this is a correlational study which does not assess whether structure shapes function or vice versa – simply that there is an association. The writing would be improved by acknowledging this and adopting a more cautious interpretation.
The sentence on line 25 was changed to: “We find that the structural connectome relates to, and likely affects, the rapid spreading of neuronal avalanches”.
The sentence on line 100 was changed to: “We show that the spatial unfolding of neural dynamics at the millisecond scale relates to the network of largescale axonal projections comprising the connectome, likely constraining the exploration of the brain’s putative functional repertoire.”
The language in the sentences has been adjusted as suggested. We have aimed to avoid jargon and keep the language as simple as possible to ensure that the paper is accessible to clinicians and neuroscientists. We believe that this framework might allow clinical researchers to test the effects of localized lesions on the largescale dynamics and, perhaps, on symptoms.
3) Some parts of the methods writing is unclear, perhaps as the whole manuscript is so short. For example, I do not understand the second method for estimating transition probabilities outlined here:
"After the initial timebin of an avalanche, we kept 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."
We apologize, and agree that the explanation was not clear. We have now added a more through explanation that reads:
“In short, we identified regions that were recruited in an avalanche after the first perturbation (i.e. the initial timebin of an avalanche). Since we did not scroll through the avalanche in time, as previously described, we considered time delays as long as the avalanche itself, while minimizing the influence of short delays. This means that the avalanchespecific transition matrix is now binary, and the ij^{th} element is equal to 1 if region i started the avalanche (i.e. it was active at the first timebin) and region j was recruited in the avalanche at any subsequent timepoint, and 0 otherwise.”
4) I am unsure how to interpret the finding that the results are almost completely unchanged by bandpass filtering the data – α, β and γ even have identical Rvalues. On one hand this could indicate that the results rest completely on wide broadband effects but given the large and well established functional and topographical differences between these oscillations it seems likely that we would expect at least some difference. Particularly as similar past papers do show different structurefunction relationships between different MEG frequency bands (https://www.sciencedirect.com/science/article/abs/pii/S1053811918320603). It would be useful for the authors to discuss this point and its potential meaning in more detail.
We agree with this point. Our results show that a correlation between perturbations and structural scaffolding exists in all frequency bands. This is only meant, again, as a demonstration that structure affects the large scale organization of the brain at many levels, and it is not only limited to specific frequencies, as shown earlier. However, the avalanches are distinct across frequency bands. For example, avalanches do not occur (necessarily) simultaneously in the different frequency bands.
We have added and briefly discussed the reference that the referee indicated. That is indeed very relevant, thank you. However, as the reviewer noted, our results are solely based on the aperiodic part of activity. In a sense, one could say that the part of the signal we used is the one that Tewarie et al. did not consider (by focusing on the AEC and, more generally, making assumptions of stationarity), and vice versa. How aperiodic, scalefree activity reconciliates with coordination that relates to steady oscillations and synchrony remains to be elucidated. We hope that this evidence can help to provide observables related to this, against which theoretical predictions could be tested. We have integrated this within a broader review of the current literature included in the discussion, as follows:
“While our findings were replicated across multiple frequency bands, structural connectivity can potentially impose frequencydependent constraints on avalanche spread. Future work should investigate frequencyspecific data to understand what leads to the emergence of avalanches and, most importantly, to the specific spatiotemporal patterns of recruited regions that defines individual (or at least groups of) avalanches in each specific frequencyband.”.
Please also consider the response to question one of reviewer two, where further addition to the discussion also addressed the point raised here (and cite the literature suggested by the reviewers).
5) It would be good to see full use of the large datasets presented here. For instance, do the results have testretest reliability across the two scans analysed per participant? Why only reproduce the finding using the HCP DWI data but not the HCP MEG data?
For the DKT analysis, we have checked the correlation, both at the individual and at the group level, for the two scans separately. The results are reported in Author response image 5:
Interestingly, when both halves are considered separately, the correlation between transition matrices and structural connectivity is comparable, providing evidence of robustness of the data. However, in both analyses, it is somewhat lower than the one we observed when taking both halves into account. Hence, we suggest this is evidence that, if the timeseries are too short, some of the structure of the data might become difficult to measure and, accordingly, more data allows further convergence due to improved SNR. With regard to the HCP MEG data, we chose to replicate the results using the HCP tractography estimates so as to rule out that tractography performed on a 1.5 tesla machine might have introduced biases invalidating the results. Hence, this analysis is not entirely a replication of the results, but was done to rule out a specific concern. On the other hand, a subjectwise replication of the tractography and MEG correlation is a replication, which is desirable, but not strictly hindering the validity of the results of the current study. We have added a short sentence mentioning this, that reads: “However, further studies in independent datasets will have to replicate the current results”.
6) The data appears to be available from the authors 'upon request and conditional on ethics approval' but the analysis code does not appear to be available online.
We apologize for this. Actually, the code for the main analysis is available on GitHub, at https://github.com/pierpaolosorrentino/TransitionMatrices. This was stated in the submission process, we thought the reviewer would have access to this information. We apologize again for the inconvenience.
Reviewer #2 (Recommendations for the authors):
I think you are definitely onto something here, but I am not quite sure if this is necessarily the full picture – yet. It's not too difficult to remedy I imagine, but I think on a technical level I'd like to query two main points on this.
I realise that this is a short communication, but it does appear to assume that no work in the M/EEG field on relating structure and function exists, any previous work been done using fMRI. I think it might be important to just give a bit more context where in the literature this might sit? Just looking around for a couple of quick examples of relating white matter tractography data to functional connectivity, there's a study relating EEG connectivity to structural connectivity (Glomb et al., 2020, Network Neurosci. 4 (3): 761787) And plenty of modelling papers which trying to look at how structure would generate observable M/EEG connectomes (eg. Tewarie et al., 2019 NeuroImage 186). I think a slightly wider look at the literature on this would be really helpful. Crucially, it's not that I don't think that this isn't novel, but rather it seems to currently seems to not acknowledge the exitance of previous work in the M/EEG field?
We apologize if we gave the impression of not acknowledging previous work. This was not our intention. We chose to aim at a brief report and not a full paper since we would like to convey our novel findings as succinctly as possible. We have now cited some of the work indicated by the reviewer, and, in the discussion, we now frame our work in more detail within the current literature, as follows:
“Consistent with our findings, two recent M/EEG studies showed that functional connectivity, as estimated using amplitudeenvelope coupling (AEC), relates to structural connectivity (Glomb et al., 2020; Tewarie et al., 2019). However, in contrast to AEC, we conducted timeresolved analyses, characterizing avalanche dynamics at high temporal resolution. Further work is needed to determine the extent to which structurefunction coupling is dynamic. To this regard, our results suggest that coupling is strongest during avalanche events, consistent with established theories (Dehaene et al., 1998)”.
The use of neuronal avalanches method, is it as invariant to leakage as you initially think it is? Leakage is a product of signal variance (c.f. O'Neill et al., 2015 Phys. Med. Biol R217). So if a signal is stronger in amplitude, the leakage increases and can propagate further. So in principle, if a seed area reaches over a threshold zscore, it could continue to raise in amplitude such that a test area might also then reach a threshold? Could other (i guess more established) connectivity metrics be able to resolve this also? I was wondering even if you assume stationary variance in your signal, you might be able to analytically estimate 'leakage' between two areas (see the O'neill paper above for an analytic expression). Could you then (potentially) have a leakage connectome and see how well that predicts your avalanches perhaps? If its considerably less than the structural, I'd be a lot more confident in what I am seeing.
We thank the reviewer for this relevant comment. Please refer to our response to reviewer one (comment 1) for additional analysis exploring the effects of leakage. We chose not to pursue the idea of using more traditional metrics (plv, pli, plm etc.) since they are based on the framework of synchronization and on the assumption of stationarity. In contrast, our analysis focuses solely on the aperiodic, fat tailed distributed part of the activity, which is what is not to be expected in a purely oscillatory system. A comment on this was added to the manuscript, as part of the broader review of the literature in the discussion.
Reviewer #3 (Recommendations for the authors):
Is the inhouse software written in Interactive Data Language (used for calculating the number of streamlines connecting each pair of ROIs and the corresponding mean tract length) publicly available? If so, please provide a web link to it.
The software is not publicly available. However, it can be provided for research use upon request to the author [MQ]
A bit of extra information about the MEG scanner would also be useful, in addition to reporting the actual number of channels that were manually removed.
Further information about the MEG system has been provided, as well as a reference that fully describes the system we used, as well as the number of channels that were manually removed:
“The MEG system we used was equipped with 163 magnetometers, and was developed by the National Research Council of Italy at the Institute of Applied Sciences and Intelligent Systems (ISASI). All technical details regarding the MEG device are reported in (Rombetto et al., 2014).”
And
“(136 ± 4 sensors were kept)”
References
Abeyasinghe PM, Aiello M, Cavaliere C, Owen AM, Soddu A. 2021. A comparison of diffusion tractography techniques in simulating the generalized Ising model to predict the intrinsic activity of the brain. Brain Struct Funct 226:817–832. doi:10.1007/s00429020022116
Dehaene S, Kerszberg M, Changeux JP. 1998. A neuronal model of a global workspace in effortful cognitive tasks. PNAS 95:14529–14534. doi:10.1073/pnas.95.24.14529
Glomb K, Mullier E, Carboni M, Rubega M, Iannotti G, Tourbier S, Seeber M, Vulliemoz S, Hagmann P. 2020. Using structural connectivity to augment community structure in EEG functional connectivity. Network Neuroscience 4:761–787. doi:10.1162/netn_a_00147
Honey CJ, Sporns O, Cammoun L, Gigandet X, Thiran JP, Meuli R, Hagmann P. 2009. Predicting human restingstate functional connectivity from structural connectivity. Proceedings of the National Academy of Sciences of the United States of America 106:2035–2040. doi:10.1073/pnas.0811168106
Rombetto S, Granata C, Vettoliere A, Russo M. 2014. Multichannel System Based on a High Sensitivity Superconductive Sensor for Magnetoencephalography. Sensors 14:12114–12126. doi:10.3390/s140712114
Seguin C, Razi A, Zalesky A. 2019. Inferring neural signalling directionality from undirected structural connectomes. Nature Communications 10:1–13. doi:10.1038/s4146701912201w
Seguin C, Tian Y, Zalesky A. 2020. Network communication models improve the behavioral and functional predictive utility of the human structural connectome. Network Neuroscience 4:980–1006. doi:10.1162/netn_a_00161
Seguin C, Van Den Heuvel MP, Zalesky A. 2018. Navigation of brain networks. Proceedings of the National Academy of Sciences of the United States of America 115:6297–6302. doi:10.1073/pnas.1801351115
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. Journal of Neuroscience 33:7079–7090. doi:10.1523/JNEUROSCI.428612.2013
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. doi:10.1016/j.neuroimage.2018.10.079
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Reviewer #1 (Recommendations for the authors):
Many thanks to the authors for a clear and thorough response to my queries. Overall, the revised manuscript is greatly improved and the changes largely address my concerns.
I must ask for one final clarification, and possible addition, to the simulation scheme looking at leakage. The revision doesn't state whether any noise is added at the sensor level between projection through the leadfields and subsequent beamforming – is any noise added here? If so, please in include any details (and my apologies if I missed this detail somewhere)
If not, then this additional step should be added to the pipeline. As it stands, the simulation without sensor noise creates unusually favourable conditions for the beamformer. Environmental and sensor noise have their own dynamics and can be correlated in time and across sensors, these correlations mask neuronal signals and can be tricky for the beamformer to completely remove.
A gold standard would be to add sensor noise from an empty room recording, or if this is impractical then some modestly scaled noise whose sensor x sensor covariance matrix an example empty room recording.
The correlations of both the 'coupled' and 'uncoupled' simulations are substantially lower than the observed value in the 'correlation of lineally mixed surrogates' figure. I suspect this is partially due to favourable beamforming from missing sensor noise, as it stands the beamformer is able to remove nearly all spurious events. Including the sensor noise would likely induce some more spurious connections in both the coupled and uncoupled cases, it would be a very strong validation of the method to demonstrate that these are still distinguishable from the observed effects.
We thank the reviewer for their important suggestion. We have now repeated the same simulation as earlier, but this time adding white noise, correlated as 1/distanceamong sensors, to the sensorlevel signal, before performing the inversion. Each channel has SNR = 4. We report in Author response image 6, to the right, the correlations retrieved with the added noise, and to the left the correlations obtained without the noise, from Author response image 3, for comparison.
As correctly predicted by reviewer one, the noise structure introduces spurious correlations, that can be observed in both the coupled and uncoupled timeseries, given the shift of both distributions to the right, as compared to the correlations retrieved without noise. However, given that in both cases the observed correlations are significantly lower than the original ones is reassuring. We have now modified the manuscript, in the methods, field spread analysis section, as follows:“Noise, correlated as 1/distance among sensors, was then added to the sensorlevel time series, with a SNR = 4.”
We have also replaced our original simulation results (no noise) with the new results that include noise.
https://doi.org/10.7554/eLife.67400.sa2Article and author information
Author details
Funding
No external funding was received for this work.
Ethics
Human subjects: All participants gave written informed consent. The study complied with the declaration of Helsinki and was approved by the local Ethics Committee (Prot.n.93C.E./Reg. n.1417OSS).
Senior Editor
 Timothy E Behrens, University of Oxford, United Kingdom
Reviewing Editor
 Diego Vidaurre, University of Oxford, United Kingdom
Reviewers
 Andrew J Quinn, Oxford Centre for Human Brain Activity, Wellcome Centre for Integrative Neuroimaging, United Kingdom
 George O'Neill, University College London, United Kingdom
Publication history
 Preprint posted: November 25, 2020 (view preprint)
 Received: February 9, 2021
 Accepted: July 7, 2021
 Accepted Manuscript published: July 9, 2021 (version 1)
 Accepted Manuscript updated: July 13, 2021 (version 2)
 Version of Record published: July 21, 2021 (version 3)
Copyright
© 2021, Sorrentino et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,545
 Page views

 254
 Downloads

 9
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
Cycling of cosubstrates, whereby a metabolite is converted among alternate forms via different reactions, is ubiquitous in metabolism. Several cycled cosubstrates are well known as energy and electron carriers (e.g. ATP and NAD(P)H), but there are also other metabolites that act as cycled cosubstrates in different parts of central metabolism. Here, we develop a mathematical framework to analyse the effect of cosubstrate cycling on metabolic flux. In the cases of a single reaction and linear pathways, we find that cosubstrate cycling imposes an additional flux limit on a reaction, distinct to the limit imposed by the kinetics of the primary enzyme catalysing that reaction. Using analytical methods, we show that this additional limit is a function of the total pool size and turnover rate of the cycled cosubstrate. Expanding from this insight and using simulations, we show that regulation of these two parameters can allow regulation of flux dynamics in branched and coupled pathways. To support these theoretical insights, we analysed existing flux measurements and enzyme levels from the central carbon metabolism and identified several reactions that could be limited by the dynamics of cosubstrate cycling. We discuss how the limitations imposed by cosubstrate cycling provide experimentally testable hypotheses on specific metabolic phenotypes. We conclude that measuring and controlling cosubstrate dynamics is crucial for understanding and engineering metabolic fluxes in cells.

 Computational and Systems Biology
 Neuroscience
Seizure generation, propagation, and termination occur through spatiotemporal brain networks. In this paper, we demonstrate the significance of largescale brain interactions in highfrequency (80200 Hz) for identification of the epileptogenic zone (EZ) and seizure evolution. To incorporate the continuity of neural dynamics, here we have modeled brain connectivity constructed from stereoelectroencephalography (SEEG) data during seizures using multilayer networks. After introducing a new measure of brain connectivity for temporal networks, named multilayer eigenvector centrality (mlEVC), we applied a consensus hierarchical clustering on the developed model to identify the epileptogenic zone (EZ) as a cluster of nodes with distinctive brain connectivity in the ictal period. Our algorithm could successfully predict electrodes inside the resected volume as EZ for 88% of participants, who all were seizurefree for at least 12 months after surgery. Our findings illustrated significant and unique desynchronization between EZ and the rest of the brain in early to midseizure. We showed that aging and duration of epilepsy intensify this desynchronization, which can be the outcome of abnormal neuroplasticity. Additionally, we illustrated that seizures evolve with various network topologies, confirming the existence of different epileptogenic networks in each patient. Our findings suggest not only the importance of early intervention in epilepsy but the possible factor which correlates with disease severity. Moreover, by analyzing the propagation patterns of different seizures, we asserted the necessity of collecting sufficient data for identifying the epileptogenic networks.