Neuronal Cascades Shape Whole-Brain Functional Dynamics at Rest

Connectome based modeling. A, The mean firing rate r and membrane potential V variables of the NMM are derived as the limit of infinite all-to-all coupled QIF neurons. Applying the Balloon–Windkessel model to V _n(t) we obtain the simulated BOLD signal B_n(t) at node n. B, The phase plane of each decoupled node (I_ext = 0) has a “down” stable fixed point and an “up” stable focus (full dots). These points are defined at the intersection of the nullclines r˙=0 (orange line) and V˙=0 (green line) where the dynamics freezes. The empty circle marks an unstable fixed point. As the external current I_ext is increased, the phase plane of the neural mass changes (see equations in Materials and Methods). In particular, the basin of attraction of the up state gradually becomes larger than that of the down state, while the fixed points move farther apart. C, D, The mouse connectome and structural connectivity W_nm were imported from the tracer experiments of the Allen Institute. The 104 cortical ROIs (corresponding to network nodes) are specified in Table 1. E, When the regions are coupled in a brain network, each node n receives an input current I_ext which is the sum of the other nodes’ firing rates, weighted by the structural connectivity. According to panel B, this input provokes a distortion of the local phase plane at node n.

Two regimes of dFC. A, Given two nodes n and m the edge CA signal E_nm(t) (orange box) is defined as the product of the z-scored BOLD signal B_n(t) and B_m(t). Averaging the BOLD-CA matrix over time we obtain the Pearson correlation across each pair of brain regions n and m (in black box, right), defining the static FC. Each column of the BOLD-CA matrix represents an istantaneous realization of the FC (iFC). B, The elements (t_i, t_j) of the dFC_e matrix are defined as the Pearson correlation between iFC(t_i) and iFC(t_j). Note in panel A the presence of transient bouts of strong BOLD-CA (e.g., in the blue boxes). During these events, the iFC remains relatively correlated for few consecutive time points, which gives rise to diagonal (yellow) blocks in the dFC_e matrix. The same CA burst (e.g., at t_i) can re-occur in time after long periods of time (e.g., at t_j), which gives rise to an off-diagonal dFC_e block (e.g., at the crossing of the dashed lines in panel B).

Two qualitatively distinct regimes of non-trivial functional dynamics. For every couple of global parameters (G, N) we calculated the dFC in a sliding window approach (dFC_w; as in Materials and Methods) and in an edge-centric approach (dFC_e; as in Fig. 2A,B). The “switching index” of each dFC_w matrix was evaluated as the variance of the respective upper triangular elements. We find two regimes of activity, named monostable and bistable, where qualitatively distinct neuroelectric organizations give rise to large-scale functional dynamics characterized by a non-vanishing switching index. In both regimes, the dFC_w and dFC_e display off-diagonal blocks, demonstrating a correlation between the functional activity at distinct times. The low global coupling G in the monostable regime (bottom left) does not guarantee a strong communication between the brain network regions, which most of the time populate the low firing rate (“down”) state. A strong noise N pushes the brain regions in the high firing rate (“up”) state for short transients. A higher value of the global coupling in the bistable regime (bottom right) promotes a subgroup of regions in the high firing rate (up) state. Low levels of noise perturb the equilibrium of the system provoking localized switching in both up →down and down →up directions (e.g., at t = 200 ms).

Mechanisms of cascade generation in the synthetic model. A, Different regions have a different fate depending on their location in the connectome. We classified the regions in five classes (D, U, J, D*, U*) according to their activity. B, Example exploration of the projected 2D phase space (top) and firing rates activity (bottom) of the “up-U” (light red), “down-D” (light blue), and “jumping-J” (green) regions. C, Distribution of the standardized firing rates in different classes. Class (J) regions have two modes but never cross the ±3 σ threshold (black lines). Class (U*) (dark blue) and class (D*) (dark red) regions dwell most of the time in the up and down states, respectively. Only in important rare occasions the *-regions cross the threshold to jump on the other side, substantially deviating from their baseline activity. The leading role of the *-regions as compared with the other classes is shown using PCA in Extended Data Figure 4-1A,B. D, Example of a cascade: when the (U*) node 63 jumps into the down state, it first drags down the node 62 (with which it shares the strongest structural link in the network). After them, other strongly connected nodes follow the trend.

RSN formation. A, top, The standardized firing rate activity in the (U*) and (D*) classes (class-specific average; dark red and dark blue, respectively) is characterized by peaks (the strongest are marked as I, II, III, IV) occurring in correspondence of cascades similar to Figure 4D. A, middle, During a cascade, we also observe a peak of BOLD-CAs, appearing as vertical strips. Many, but not all, edges are recruited. A, bottom, The blocks in the dFC_e matrix appear in correspondence of CA events, showing that these bursts generate stable epochs of FC correlated in time. B, In each selected epoch (I, II and III, IV), the large firing rate cascades trigger the jump of other nodes away from baseline activity (circled in black) and promote specific functional hubs at the BOLD level, represented by colored nodes in the network plots. A functional hub is defined by the components of the first leading eigenvector (linear combination of brain regions explaining most of the variance in the data; eigenvalue λ > 0.41) associated to the iFCs at times tI,tII,tIII , and t_IV, respectively. The most representative hub regions are depicted in yellow. Gray regions have been excluded as they do not contribute substantially. Only the edges with the highest CAs are displayed. Importantly, CA events generated from neuronal cascades at specific sites support distinct functional networks which are not correlated among themselves (e.g., no off-diagonal dFC_e block between I and III).

Neuronal cascades and neuronal avalanches. A, Standardized EEG activity extracted from a resting-state human EEG/fMRI dataset (top). The activity is binarized assigning a unitary/null value every time the activity in a region is above/below a certain threshold (e.g., ±3 σ; black lines). The obtained binary raster plot (bottom) is characterized by intermittent epochs of deviations from baseline activity. Neuronal avalanches are defined as consecutive deviations from baseline activity (e.g., red box). B, We extract the global magnitude of the deviations from baseline (top, gray signal) by summing the binary EEG raster plot over the ROIs. This signal is convoluted with a Gaussian kernel [width = 1 BOLD TR] and downsampled to obtain the same resolution of the BOLD activity, which defines the neuronal cascades signal (blue). Neuronal cascades can be thought of as clustering of high magnitude avalanches, whose occurrence in time is not homogeneous (bottom).

Neuronal cascades drive the functional dynamics. A, Example of neuronal cascades in the monostable and bistable synthetic regime and for a representative subject of the empiric EEG/fMRI human dataset. B, The BOLD-CA in simulated and empiric mouse and human datasets (middle panels) are characterized by sudden collective events involving large network parts (vertical stripes). The root sum square of BOLD-CA across all edges (RSS, green lines, top panels) defines the global CA amplitude signal for each dataset. Concurrently, the dFC_e matrices (bottom panels) display both diagonal and off-diagonal blocks, remarking the non-trivial re-occurrence of the same stable functional network at distinct times (see Fig. 2B). At a visual inspection, the BOLD-CA events happen in coincidence with dFC_e blocks and, most notably, the neuronal cascades and RSS signals (blue and green lines in panels B, C, respectively) co-fluctuate in most instances.

A, The largest BOLD-CAs events (CA, above the 98th percentile of the RSS; Fig. 4C, green line) are distinguished from non-events (nCA, below threshold). We report the synthetic and empiric correlations between iFCs at times within CA events (left in every panel), between CA events and non-events nCA (center of panels), within non-events nCA (right of panels). These correlations are by definition the off-diagonal values of the dFC_e matrix (see Fig. 2A,B). The distribution of the correlations within events is wider and explains the greatest off-diagonal correlation values of the dFC_e across all the synthetic and empiric datasets. This principle is explicitly shown in panel B, where the original dFC_e extracted from an empiric human trial (top) was sorted according to increasing RSS (bottom), leading to the clustering of high correlations toward high CA times. This shows that most of the non-trivial temporal correlations involve CA times falling in the last quartile of the RSS (above the 75th percentile, central green line). Thus, the strongest CA events drive the dynamics of FC.