Hemodynamic cortical ripples through cyclicity analysis

The brain spontaneously generates neural activity even in the absence of any sensory input, motor output, or attention-demanding cognitive tasks. Such intrinsic activity “at rest” can be indirectly measured as infraslow (<0.1 Hz) fluctuations in blood-oxygen level-dependent (BOLD) signals using functional magnetic resonance imaging (fMRI). Despite its “noisy” appearance, many studies have shown that spontaneous BOLD activity carries meaningful information, having a dynamic structure that can interact with perception, cognition, and behavior (Allen et al., 2014; Cohen, 2018; Fox & Raichle, 2007; Gonzalez-Castillo et al., 2019; Sadaghiani, Poline, Kleinschmidt, & D’Esposito, 2015; Shahsavarani et al., 2023; Shine et al., 2016; Yu et al., 2015). Conventionally, using a measure known as resting-state functional connectivity (FC), this structure has been characterized as synchronous patterns of BOLD signals in bilaterally symmetric, distant brain regions, forming time-resolved functional neural networks distributed across the brain.

An emerging corpus of work has revealed that infraslow spontaneous fluctuations of BOLD signals are organized not only spatially into neural networks but also temporally in an orderly fashion (Mitra et al., 2018; Vidaurre, Smith, & Woolrich, 2017; Zhang, Watrous, Patel, & Jacobs, 2018). That is, there are time-resolved transitions among resting-state networks, which are not random but follow a specific sequential order. It is posited that a global structure guides the switching among the neural networks, which can be explained by cortical traveling waves, sequentially activating these networks (Matsui, Murakami, & Ohki, 2016; Munn, Müller, Wainstein, & Shine, 2021; Raut et al., 2021). Previously, spontaneous and stimulus-evoked traveling waves propagating across the cortex were reported using electroencephalography (EEG) and electrocorticography (ECoG) (Massimini, Huber, Ferrarelli, Hill, & Tononi, 2004; Muller, Chavane, Reynolds, & Sejnowski, 2018). In these cases, cortical waves have been primarily noted at higher frequencies (>0.1 Hz); however, cortical traveling waves can also be identified at lower frequencies. Infraslow neuronal activity (<0.1 Hz) was first observed in the rabbit cortex’s local field potential (Aladjalova, 1957). Later, full-band EEG was used to detect this activity in humans (Vanhatalo et al., 2004). These infraslow oscillations are a fundamental feature of brain function, which correlate with higher frequency oscillations and influence behavior (Lörincz, Geall, Bao, Crunelli, & Hughes, 2009; Palva & Palva, 2012; Sihn & Kim, 2022). In addition, it has been demonstrated that these oscillations are correlated with infraslow hemodynamic activity and underlie changes in fMRI resting-state networks (Fernandez & Lüthi, 2020; Grooms et al., 2017; Hiltunen et al., 2014). Notably, using wide-field voltage-sensitive dye imaging, Chan, Mohajerani, LeDue, Wang, and Murphy (2015) have shown that the infraslow spontaneous neural fluctuations recapitulate cortical motifs at higher frequencies. Using wide-field optical calcium imaging in transgenic mice, Matsui et al. (2016) detected infraslow global traveling waves across the cortex. Importantly, they showed that the phase of these cortical traveling waves regulated the resting-state FC patterns of hemodynamic data. Using concurrently recorded intracranial EEG and fMRI measures, Mostame et al. (2022) showed that while the two measures shared similar connectome configurations, they did so at asynchronous points in time, suggesting that hemodynamic and electrophysiological signals capture distinct aspects of connectivity rather than measuring the same phenomenon.

The primary finding of this paper is the detection of traveling waves propagating across the cortex by expanding the previously presented cyclicity analysis (CA) pipeline (Zimmerman, Abraham, Schmidt, Baryshnikov, & Husain, 2018).

We report three significant novel observations:

• A time-resolved temporal ordering among the resting-state activity of brain regions was identified with a wave propagating along a spatial axis, resembling cortical hierarchical organization.

• The dynamics of temporal ordering in the task paradigm were primarily characterized by distinct, short-term “bursts” in contrast to gradual changes.

• The direction of wave propagation repeatedly shifts during certain tasks, a phenomenon we term as “reversals” to contrast with the unidirectional steady dynamics observed in the resting state.

Our work also uses novel methodological approaches. While CA was conceptualized and tested in the context of cyclic phenomena—defined as those which follow a trajectory over and over again, in a contrast with the periodic phenomena, which repeat after time shift by a well-defined interval (Baryshnikov & Schlafly, 2016)—we discovered that in the systems we analyze here, the better interpretation is that of a ripple sweeping through the terrain, followed by an interval of calm.

From the perspective of reparametrization-invariant data analysis, where CA belongs, the length of the calm interval is immaterial. The key difference with the cyclic situation is the fact that, at times, most of the network nodes are at rest. Exhibiting unexpected robustness, CA proved to be capable of detecting the presence and describing the structure of these ripples.

Our second methodological innovation involves the systematic use of oriented areas (discussed in the Oriented Areas section) as functions of time. By design, the oriented areas, or other iterated integrals, are numeric quantifiers. Each of them associates a single number to a (multidimensional) time series, defined over a time interval. However, one can also consider these time intervals as expanding in time by making their upper bound variable and evaluate the oriented areas on subintervals [t_s, t], t_s < t ≤ t_f of the total observational time interval [t_s, t_f]. This turns the oriented areas into the functions of t. A careful investigation of how the oriented areas evolve can, as it turned out, provide significant insights into the nature of the process. In what follows, we apply the CA computational pipe without any modifications to the time series extracted from the regions of interest (ROIs).

The computational procedure of CA, which underlies our findings, exploits pairwise temporal relations between the time series corresponding to different ROIs. These temporal relations are quantified through the so-called iterated path integrals as the baseline quantifiers of our observed time series (Chen, 1958). The applicability of the CA pipeline to the analysis of fMRI traces was demonstrated in various studies (Baryshnikov & Schlafly, 2016; Shahsavarani, Abraham, Zimmerman, Baryshnikov, & Husain, 2020; Zimmerman et al., 2018).

Our approach relies on the iterated integrals of the lowest degree, where the outputs are not reducible to the functions of the increments of the time series, that is, on the iterated integrals of order two. These integrals, essentially, represent the oriented areas encompassed by the two-dimensional (2D) projections of the trajectory. One might also consider iterated integrals of higher orders. While they tend to become increasingly corrupted by noise, the idea of using iterated integrals of all orders for data analysis is advanced by T. Lyons and his school (Bonnier, Kidger, Arribas, Salvi, & Lyons, 2019).

As we have shown previously in Zimmerman et al. (2018), a positive oriented area a_kl > 0 is heuristically indicative of x_k being the leading indicator to the follower x_l, and conversely, the negative a_kl < 0 indicates that x_l is the leader and x_k is the follower. Note that the oriented areas generated by perfectly correlated or anticorrelated signals vanish, as well as the area spanned by signals with nonoverlapping supports.

The nonzero eigenvalues of the lead matrix consist of the purely imaginary pairs ±iλ_s—for some, real λ_s—which we assume are ordered by size: ∣λ₁∣ ≥ ∣λ₂∣ ≥ …. The eigenvectors corresponding to the conjugate eigenvalues ±iλ_s are complex conjugate (and necessarily have nonvanishing imaginary components).

How fast the absolute values of the spectrum of the lead matrix decrease is a heuristic indicator of the dominance of the temporal interaction network by a small number of dynamic patterns. Observing just a few pairs of large (in absolute value) eigenvalues ±iλ_s indicates a good approximation of the observations by just a few propagating ripples.

In some classes of cyclic (not necessarily periodic) signals—defined and referred to as the Chain of Offsets Model (COOM) below—the components of the eigenvector corresponding to the eigenvalues with the largest (in absolute value) eigenvalue ∣λ_k∣ line up along an ellipse in the complex plane in cyclic order corresponding to the one according to which the time series follow each other in the underlying system (we provide further details in the Chain of Offsets Model section).

Thus, we pick the eigenvalue with the largest absolute value ∣λ_k∣ and select one of the corresponding eigenvectors. We will be referring to such an eigenvector as the leading eigenvector. The collection of the components of the eigenvector (these components are complex numbers that we consider as a collection of indexed points in the complex plane) is referred to as the constellation. Each of the points in the constellation corresponds to one of the ROIs in our sample. Note that there is some ambiguity as those eigenvectors come in complex conjugate pairs, but in practice, this does not complicate matters, as the signs of the entries of the lead matrix allow one to easily distinguish between two possible orientations.

Under the COOM assumption, we interpret the angular order of the points in the constellation as indicative of the relative order of the nodes—or ROIs—in the wake of the ripple sweeping through the network. (In Section 1 of the Supporting Information, we show, using simulated wave propagation, that in an excitable network with lags, this heuristic remains valid.) It should be noted that, in general, not all nodes would be receiving a particular signal, and for those nodes, the corresponding components of the leading eigenvector will cluster near the origin. This implies that meaningful analysis should first single out the largest parts of the constellation (as defined by an intrinsic metric) and focus on those. In the following sections, we show that this leads to a metric defined by quadratic forms.

As described in detail in the Materials and Methods section, we generated time series describing the BOLD activities in the 68 ROIs (i.e., from 34 bilateral pairs). These time series were obtained by averaging time series over voxels within those ROIs (defined in Table 1) in each resting-state fMRI (rs-fMRI) recording with a duration of t_f = 864 s. The time series were then centered (setting their time averages to 0) and normalized by their quadratic variation, that is, the sum of squared differences. Oriented areas a_kl(·) were then computed using a discretized version of Equation 1 (see comments in the Discretization section).

For the full lead matrix A, the spectral decomposition was determined and the eigenvector v_* corresponding to the largest in absolute value eigenvalue iλ was chosen. There may be an ambiguity as to which of the two complex conjugate eigenvectors corresponding to the opposite eigenvalues ±iλ to chose. However, selecting a leader-follower pair of components x_k, x_l with a well-defined phase shift, and ensuring that the orientation of the pair of complex numbers v_k, v_l—the k-th and l-th components of the eigenvector v—is positive, resolves it.

The components of this eigenvector are complex numbers (indexed by the 68 ROIs). The scatterplot of these complex numbers are referred to as constellation; examples of them are shown in Figure 1.

In the case of a perfectly harmonic time series, such a constellation would necessarily align with an ellipse (referred to as the constellation ellipse, CE), except for the components corresponding to the nodes not activated during the passage of the excitation through the network. Those components would collapse to the origin (Baryshnikov & Schlafly, 2016). For a noisy signal, it is beneficial to separate the activated nodes from irrelevant ones by selecting the components farthest away from the origin. It should be emphasized, however, that the distance should be measured with respect to the quadratic form defining the CE, not the standard Euclidean metric.

Naturally, these quadratic norms ∣·∣_q differ from recording to recording. One can, however, construct an average ranking over all recordings, derived from the constellations associated with the leading eigenvalue of the lead matrices. The points of the low average rank (i.e., of large ∣·∣_q distance to the origin) are expected to be exposed to the presumed excitation rippling through the cortex, while those of high average rank (i.e., closer in ∣·∣_q distance to the origin) are heuristically not affected by the ripple.

Again, these ranks are computed in terms of the constellations corresponding to the leading eigenvalues in each recording; one can also consider the ranks corresponding to the second highest eigenvalue and so on, recovering potentially independent ripples. Figure 2 summarizes the average rankings for each ROI corresponding to the first eigenvalue (the results for the second eigenvalue is provided in the Supporting Information).

From hereon, we focus primarily on the leading eigenvalues (with the largest absolute value, i.e., λ_1,2 in our convention), but one can extend the same procedures to the next largest in magnitude eigenvalue as well.

Using the average ranking chart in Figure 2, we can choose a cutoff, akin to choosing a knee point from an eigenvalue plot in principal component analysis, to limit the number of ROIs affected by the leading (i.e., corresponding to the leading eigenvalue) ripple for further analysis.

Once a subset of ROIs is chosen in this way, one can estimate the ordering of those significant ROIs by examining the arguments of the complex valued entries in the eigenvector corresponding to them. As explicated in the Ripples section, if the nodes are excited by similar waveforms delayed by a collection of the phase shifts, CA is able to extract those phase shifts, and, correspondingly, an approximate order in which the ROIs are excited as the ripple passes through the network (Baryshnikov & Schlafly, 2016; Zimmerman et al., 2018).

This phase shift extraction is especially intuitive when the constellation is clustered in a part of the complex plane, as on the bottom left display of Figure 1. In this situation, one can reposition the axis corresponding to the zero phase so that it points away from the center of mass of the constellation. In this way, the ordering of the shifts coincides with the ordering of arguments of the (complex) components corresponding to the ROIs, normalized to vary from −π to π. This is shown on the bottom right display in Figure 1, which is obtained from the bottom left display by the rotation aligning the center of mass of the constellation with the positive real ray.

Therefore, for each eigenvalue, and each threshold separating the ROIs with high excitation levels with respect to this eigenvector, one derives the order in which the ripple captured by that eigenvalue excited the nodes.

Again, these obtained orderings change across recordings, being impacted by the subjects’ and experimental conditions’ variability and noise. However, the statistical distributions of those phase shifts are indicative of a well-defined ripple, propagating through the cortex. To this end, we derived the histograms of the frequencies with which any given ROI appears in a particular position in the ordering of the phase shifts derived from the constellation. Those histograms, depicted as the violin plots in Figure 3, are shown for the ranking thresholds of k = 20 on the top left display of Figure 4. The ordering of the ROIs is by the average positions in the order of phase shifts.

It is immediate that these frequencies are indicative of a ripple passing through the selected significant ROIs. Among remarkable features of this plot is the symmetry: The bilateral pairs of ROIs are excited at the same stages of the ripple’s passage. Anatomically, the direction of the imputed ripple is shown on the top right display of Figure 3.

Changing the cutoff thresholds does not qualitatively change the observed pattern. The lower part of Figure 3 shows the result for the cutoff at k = 24 ROIs. Again, the excitation is symmetric, and the patterns of propagation are anatomically similar to what was observed for k = 20 ROIs.

If one attempts to repeat the analysis of the excitations’ passage using the higher eigenvalues (i.e., ±iλ_s, s ≥ 2), one encounters loss of directionality. The violin plots become disordered and lacked any discernible structure; see Figure 4 with the violin plots for the second and third largest eigenvalues.

The overall increment of a_kl(t) over the interval of observation is the a_kl entry of the lead matrix A, but how this increment is reached supplies us with additional information. While the overall trend or increment is a succinct indicator of the average leader-follower relationship, the dynamics of the oriented area signal can provide insights into the dynamics of the interaction between the corresponding brain regions.

We noted primarily three kinds of dynamics that oriented area signals exhibited from the dataset.

1. Neutral dynamics with low-level fluctuations near zero.

2. Steady, nearly monotonic overall increment or decrement.

3. Steady increment or decrement punctuated by role reversals in the leader-follower relationship.

The three kinds of dynamics are illustrated in Figure 5 using examples from the data.

One would expect the first kind of dynamics for pairs of time series that are strongly correlated (or anticorrelated), and not unexpectedly, we observed them between the pairs of bilateral ROIs. The steady dynamics of the second type manifested itself frequently in the recordings of the resting-state sessions. More interestingly, the third intermittent kind of dynamics was observed in a large fraction of the task datasets. One can note, by observation, the timings of the reversals of oriented area accumulations, and it becomes readily apparent that they appear in a regular fashion and are quite localized.

(We should remark that the behavior exhibited on the bottom row of Figure 5 is unlikely to be reproduced by a low-complexity autonomous 2D system and is highly indicative of external excitations, causing the reversals of the oriented areas.)

The next two subsections deal with the detailed analysis of these observations, where we devise some tools to quantify these phenomena.

To further examine the phenomenon of dynamic fluctuations observed in the oriented area signals, we accumulated over all pairs of ROIs the time stamps at which a change in the direction of the oriented area occurred. While, in general, this happens at the (local) maxima and minima, one can filter out the small-scale jitter using the toolbox of Persistent Homology, briefly outlined in the Section 2 of the Supporting Information. One does this by coupling the local minima and maxima in a well-defined way. The heights b, d of such coupled minimum and maximum are thus referred to as the birth and death of a feature and are visualized as points on the (b, d) plane, as shown in Figure 6. We interpret the features with the largest span (i.e., the points furthest from the diagonal {b = d}) as indicative of the reversals. The coupled pairs whose span d − b is high are said to be persistent features in the data (see Ghrist, 2008 for more details).

We mark the time stamps of the minima of these persistent features. Combining all time stamps of such persistent features across all oriented area time series (corresponding to pairs of ROIs), we visualize them using a carpet plot as shown in Figure 7 (for task paradigm data).

One prominent feature of the social task plots is the regular patterns occurring approximately every 52 frames (i.e., every 0.72 s × 52 = 37.4 s). The documentation (Barch et al., 2013) revealed that the intervals between social task segments were 38 s. A similar analysis of the motor task data showed approximately 11–13 repeating epochs within the duration of the observation. Once again, the documentation recorded that for motor tasks, the participants were presented with stimulus in 13 blocks (10 movements, each were 12 s long, and with three 15-s fixation blocks). The strongest indications of a reversal consistently appeared at the onset and offset of the task blocks, which may reflect a switching between task-positive and task-negative systems at these transitions between task and rest.

A similar analysis of the resting-state sessions revealed that pairwise activities, as indicated by the above outlined analysis, were, indeed at rest; see Figure 8.

To analyze burstiness in the time series, we extended the concepts outlined above to arrive at a mechanism for quantifying the difference between the resting and task paradigms, as evidenced by the activity of pairs of brain ROIs. One salient observation that arose from the study of oriented areas as functions of time was that activity between some regions tended to be much more consistent than others. These are exhibited by a gradual increase in the overall oriented area integral. The less-consistent pairs, on the other hand, contributed to the total increment of the oriented area in jagged “bursts.” Figure 9 illustrates this phenomenon in its top two panels. The top left image shows a pairwise oriented area integral dominated by the jumps from the deep local minima. The top right image, however, shows another pair of brain regions, where the total accumulated oriented area is similar to that of the left display, but does that accumulation in a more gradual fashion.

We have marked on either plot the 10 longest persistent features observed in the time series calculated according to the technique described in the previous section. Using this heuristic, we define the burstiness of the oriented area function as the total sum of the lengths of the 10 longest bars in the persistence diagram. (We remark that the results are robust with respect to the number of the longest bars chosen for analysis: Selecting 8 or 20 longest bars leads to essentially identical outputs.)

Remarkably, the distributions of the burstiness, computed for all pairs of ROIs and scans, are dramatically different for the resting-state and the tasked sessions. As one can see in the bottom panel of Figure 9, the burstiness for the resting state is far lower than that of either motor or social tasked sessions. On the other hand, the distributions of burstiness for the motor and social tasks are very close. Note that all time series were energy-normalized, so that the low burstiness of the resting state cannot be attributed to different scales of the cortical activity.

We believe that CA stands as a promising avenue for gaining novel insights into the dynamics of FC, thereby enhancing our understanding of brain networks and their interactions. Our research findings substantiate this perspective by revealing a distinctive temporal ordering in the resting-state activity of brain regions, particularly along a spatial axis suggestive of hierarchical cortical organization. The results provide empirical support for the existence of cortical waves propagating between transmodal areas and the primary visual cortex, highlighting the relevance of CA in uncovering these patterns. Additionally, the identification of short periods of sustained lead-lag relationships in oriented areas further emphasizes the time-dependent flow of information underlying the dynamic aspects of hemodynamic measures. This temporal organization that is characterized by abrupt “bursts” rather than gradual changes underscores the dynamics of the brain’s functional architecture in the task paradigm. Importantly, we observed noteworthy shifts in the direction of wave propagation during specific tasks, presenting repeated reversals that contrast with the observed steady patterns in the resting state.

There are several other well-regarded computational techniques aimed at the discovery of the temporal interrelationships between brain regions. The lagged correlations method of Mitra, Snyder, Blazey, and Raichle (2015) infers, as its first step, the relative lags between pairs of time series, defined as the time shift that maximizes the correlation coefficient between them. (The resulting skew-symmetric matrix can be considered a nonreparameterization-invariant version of our lead matrix.) Next, they adjust this matrix by removing the mean lag for each time series. Then, they extract the leading singular vector so as to determine the sequence of lags best compatible with their modified lag matrix.

The approach of extracting the leading singular vector from a delay matrix is used broadly, with differently defined delay matrices. For a different example, see Gu et al. (2021), where the time series are sliced into segments by the local minima of the overall brain activity, and then within each such segment, the delays are detected, again by analyzing the local peaks of activity for each ROI.

In all these methods, the inferred time lags are the key element used in the subsequent numerical procedures. As such, they are inherently dependent on the consistent rhythm of the observed processes. Lagged correlation methods rely on interpolation and windowing to capture the dynamics of FC, which are sensitive to time delay estimation methods, autocorrelation (Shahsavarani et al., 2020), sampling variability (Raut et al., 2020), and parameters such as window length and window shift (Shakil, Lee, & Keilholz, 2016). In contrast, CA offers an approach with finer granularity resolution, and without the aforementioned sensitivities or assumptions regarding stationarity, state duration, or state transition, all following the principled reliance on reparametrization-invariant tools.

Further, the spectral analysis under CA is critically different from that of the extant methods. To extract the temporal ordering of the signals, we rely on the inherently complex valued eigenvalues of the skew-symmetric lead matrix (see the Cyclicity Analysis section), while other methods preprocess it to extract its singular vectors, which are then interpreted as the delay vector. From these perspectives, the CA is a method of extraction of the temporal ordering of the observed signals. which is radically different from the existing ones.

Although cortical waves have been noted at much higher frequencies, detecting them in fMRI data at subhertz frequencies with CA tools is a novel observation. At higher frequencies, cortical waves of oscillatory activity are thought to contribute to a spatiotemporal framework for neuronal communication by coordinating a range of synchronizing patterns (Das et al., 2022; Muller et al., 2016, 2018; Roberts et al., 2019; Rubino, Robbins, & Hatsopoulos, 2006; Townsend et al., 2015; Zhang et al., 2018). Although spike trains of individual neurons have shown to be coupled to both slow and fast time scales (Okun, Steinmetz, Lak, Dervinis, & Harris, 2019), the temporal duration of these more typical cortical waves tends to be on the order of tens to hundreds of milliseconds rather than the subhertz frequencies, which BOLD signals are able to resolve. Nevertheless, over the past decade, a growing body of studies has observed BOLD signal flow across cortices (Gu et al., 2021; Hindriks, Mantini, Gravel, & Deco, 2019; Ma & Zhang, 2018; Majeed et al., 2011; Mitra, Snyder, & Raichle, 2020; Raut et al., 2021). Taken together, these studies provide evidence that there is likely a propagation of the BOLD signal throughout the brain that dynamically carries or reflects longer periods of stable information flow between brain regions. The method that we have outlined here should be a useful tool for investigating these dynamics due to its invariance to the arbitrary nonlinear reparameterizations of the timeline, a feature lacking in the correlation-based methods.

One striking pattern revealed by our analysis is the predominant role of distinct, short-term “bursts” in lead-lag relationships between pairs of brain regions, as opposed to prolonged, consistent stretches of moderate associations. This should be compared with the recent paper by Esfahlani et al. (2020), where the authors investigated the contributions of moment-to-moment BOLD activity to the overall pattern of FC. Similar to our observation, they saw that only a small fraction of frames in the time series explained a significant amount of the variance in the network connectivity. These fluctuations corresponded with activity patterns aligned with fluctuations in the level of default mode and attentional control network activity, which are often viewed as in opposition to each other (Anderson, Ferguson, Lopez-Larson, & Yurgelun-Todd, 2011).

This opposition between the default mode and attention networks has a strong overlap with the idea of principal gradients of macroscopic cortical organization in the brain (Margulies et al., 2016; Samara, Eilbott, Marguilies, Xu, & Vanderwal, 2023). According to this framework, a topography of connectivity patterns is reflected in the spatial distances between the so-called “higher” areas of the brain, where the more multimodal, abstract, predictive information is encoded, and “lower” areas, such as the primary sensory/motor regions. This primary connectivity gradient predicts the positions of canonical resting-state networks, which are viewed in this framework as reflecting representational hierarchies rather than distinct modules. In other words, resting-state networks are reflective of the temporal phase of propagating patterns. The functions associated with various cortical networks were correlated with the level of the hierarchy of sensory or motor processing. Also, Raut et al. (2021) observed slowly propagating waves of activity measured with fMRI in humans, which were associated with cyclic fluctuations in arousal. They then replicated the result in macaque monkeys using hemisphere-wide ECoG. Critically, they found that functional networks maintained phase shifts relative to one another in their relation to autonomic arousal, rather than specific time delays. This result validates the usefulness of a time reparametrization-invariant analytical method like the one we present here.

The results of CA, as applied to the HCP dataset, revealed some unexpected findings, which should be investigated in future studies.

In the context of the resting-state scans, we see ordering in the direction moving along the cortical hierarchy, through association areas toward the visual regions (see Figure 3). The precentral gyrus activity leads, and the lateral occipital cortex follows. Compared with the intermittent bursts in the task state, wave here follows a reverse hierarchy, moving from association areas toward the visual processing regions. Interestingly, our results echo those by Gu et al. (2021), which found that while the global propagation pattern follows the cortical hierarchy at the network level, there may be a reverse hierarchy within the sensorimotor networks (e.g., ending at the lateral occipital cortex rather than the pericalcarine cortex). The temporal ordering between these regions is steady and nearly monotonic. It is possible that during resting state, brain activity is mostly internally generated as in mind-wandering. Possibly because there is no stimulus to attend to, activity mostly occurs in the higher-to-lower direction in the hierarchy, reflecting endogenous processing.

In the context of the social task (see Figure 7), we observed a strikingly different pattern. This task consisted of five blocks, each consisting of a period of watching a video and then judging whether the objects in the video appear to perceive the other objects’ feelings and thoughts. In this task, we still observed strong directional constraints, but which brain area is the leader and which is the follower is shifting over time. This may be interpreted as periods of time that switch between internally generated (endogenous) activity and externally generated (exogenous) activity, which may correspond with periods of watching the videos and periods of making judgments of the objects’ intentions within the video. The direction of the leader-follower relationships is also much burstier compared with the direction during the resting-state data.

The pattern of activity related to the motor task was also bursty and more similar to the social task than the resting-state task (see Figure 7). Here, the participants were presented with 3-s cues to perform 12 motions. There are 10 motor blocks and three 15-s fixation blocks in each run. Like in the social task, we observed switches between the leaders and the followers during the run. However, it is more difficult to match the activity observed to events in the task. This may be because the cuing events (3-s long) and motor movement events (12-s long) are too short to be adequately captured by the low-frequency movement across the brain. We may be observing some aliasing in the directionality of the lead-lag signal. It is also possible that in the motor task, there is more separability between activities belonging to hierarchical visual processing and motor processing.

Our analysis identified a temporal ordering that primarily propagated in opposite directions along an axis resembling the principal gradient of the brain’s spontaneous activity identified by Margulies et al. (2016). However, it should be noted that our analysis focused on analyzing the first eigenvector of the lead matrix, and this favors patterns of activity from subjects having high ∣λ₁/λ₃∣ ratios in the lead matrices. Thus, our choice likely selected for subjects with strong leader-follower activity along a single direction. It is possible that the analysis of the other eigenvectors would reveal a pattern of temporal ordering moving in other directions, reflecting other kinds of neural dynamics, such as the signaling across hemispheres.

In this paper, we deployed CA to detect transient processes in the brain. This new approach complements recent advancements in effective or directional connectivity research. The method outlined here has intrinsic advantages due to its lack of assumptions regarding the temporal properties of the BOLD signal dynamics. It does not assume the stationarity of the time series or specific properties of latency, state duration, or state transitions, which could bias correlational, spectral, or lag-based approaches. We have shown in group data a primary propagating wave of BOLD activity during resting state along a spatial axis related to cortical hierarchical organization. We further observed that these propagating waves appear to switch directions in a task-dependent manner.

The findings reported here open a wide range of future research directions. Future work should apply CA techniques to other measures of human brain activity that are more reflective of direct neural activity, including EEG, ECoG, the fast optical signal, and magnetoencephalography. Applying CA to multiple sensing modalities may help clarify the relationship between patterns observed in the fast temporal domain of neuronal activation and the longer duration patterns observed in the BOLD signal, which may be more reflective of broader stable states of the whole-brain function.

From the Human Connectome Project 1200 Release (Glasser et al., 2013; Van Essen et al., 2012), we considered 889 denoised minimally preprocessed participants who completed all the structural, resting-state, and task fMRI sessions using a customized Siemens Skyra 3 T scanner. Of those, 27 were excluded from the analysis for the segmentation issues noted in the HCP Quality Control process or functional preprocessing errors reported in the HCP Data Release Updates. After exclusion, 862 participants remained for analysis who ranged in age from 22–45 years and included 464 females. We used the Connectome Workbench software to extract time series from ROI information in the fMRI parcellations available for download.

The ROIs extracted were 34 in number (bilateral; therefore, total N = 68) and are listed in Table 1. The parcellated output files were further processed in Python/R to create time courses for the regions listed in Table 1. At a time repetition of 720 ms, for the rs-fMRI scans (numbering four per participant), this resulted in a set of arrays D, with elements of dimension 68 × 1, 200 that were fed into the CA pipeline. (Here, 1, 200 is the number of fMRI frames; the total session’s duration, therefore, being 14.4 min). Further details of the scan protocol are available online from the HCP Project Consortium. The top left panel in Figure 1 shows a representative of BOLD signals in D.

The HCP dataset includes fMRI scans recorded under different task paradigms. Such scans recorded during the progression of certain motor and cognitive tasks were also included in part of our analysis. The two task protocols are as follows—for the motor task (Barch et al., 2013):

For CA, lead matrices were generated from the time courses and their eigenstructure were analyzed. These matrices have the dimension 68 × 68, with each (i, j) entry denoting the average leader-follower relationship between ROI i and ROI j. See the right panel in Figure 1 for a representative example of a generic lead matrix. The ∣λ₁∣/∣λ₃∣ ratio (here, λ_k are eigenvalues of the lead matrix), a measure of the dominance of the rank 2 approximation (described under the Chain of Offsets Model section), was computed for all lead matrices. See Figure 10 for relevant statistics. (Note that λ₁ = $λ¯$₂ since the lead matrix is real and skew-symmetric by construction.)

Recall that the CA pipeline builds on the iterated integrals of second order (Zimmerman et al., 2018).

Such iterated integrals, not expressed in terms of the iterated integrals of the first order, are known as the oriented areas. It relies on an intuitive interpretation of the oriented area encompassed by a curve in the plane, parametrically represented by a pair of trajectories, x and y. Consider a function of time x(t) shaped as a pulse (blue curve in Figure 11), as well as its time-shifted copy, y(t) (red curve). We interpret the interrelationship between x and y as a leader-follower one: x leads and y follows. Now, to express this relationship in a reparameterization-invariant fashion, we render it using just the parametric plot (taking into account the orientation of the resulting curve in the (x, y) plane). The key observation now is that when plotted against each other, these two trajectories enclose a regionRof the plane of positive area (as shown in Figure 11).

Consider a closed parametric curve (x(t), y(t)) ∈ ℝ², t ∈ [t₁, t₂] with piecewise differentiable functions x, y as the components. This curve partitions the plane into the open connected domains. Within each of these domains, the winding number of the curve around a point is well-defined: One can imagine an observer located within such a domain, tracing the point of the trajectory; the winding number counts the total number of turns the observer makes (if the curve winds around the point counterclockwise, the point’s winding number is positive; if clockwise, negative). The sum of the areas of those domains, weighted by their winding numbers, is called the algebraic or oriented area encircled by the curve.

Note that the quantity in Equation 3 behaves differently from the correlation coefficients. While the correlation of two signals is largest when they are proportional, the oriented areas are intrinsically antisymmetric, so that it necessarily vanishes on a pair of proportional functions. On the other hand, if the time supports of the two signals do not overlap in time, their oriented area measure is also zero.

In practice, the integrals are replaced by the discretized sums requiring an estimate of the resulting error. From the analytical viewpoint, one needs to evaluate the dependence of the integrals (Equation 1) expressing the oriented areas on the variations of the integrands x_k, x_l: how they change when passing to the functions $xk′$⁠, $xl′$⁠. One can a priori assume that the error terms e_k ≔ $xk′$ − x_k, e_l ≔ $xl′$ − x_l have bounded C¹ norms C_k, C_l, respectively, and vanish at the discretization points Δn, n = 0, …, N.

In isolation, oriented areas just provide information about the pairwise relationships between two time series. To capture the collective information a multidimensional trajectory represents, we need to aggregate these data in a certain way.

For a trajectory in an n-dimensional Euclidean space ℝⁿ, we arrange the oriented areas into a square n × n matrix (where n is the number of time series observed). This matrix, whose (k, l)-th entry is the oriented area spanned by the pair x_k, x_l of the time series, is referred to as the lead matrix (Baryshnikov & Schlafly, 2016). We remark that this matrix is obviously skew-symmetric, implying particularly that its eigenvalues form pairs of purely imaginary numbers, ±iλ_s, λ ∈ ℝ, s = 1, 2, …, n, and the corresponding eigenvectors form complex conjugated pairs with necessarily complex valued components.

We will be referring to this model as the COOM.

One should remark that the reparameterization-invariant nature of our tools makes COOM applicable when the waves coming through the substrate parameterized by ξ need not be necessarily overlapping. The whole formalism works for a sequence of solitary waves, between which the system comes to rest.

Therefore, the spectral decomposition of the lead matrix can lead to the recovery of the order in which the signals are represented by the components of the time series x(t).

Yuliy Baryshnikov was partially supported by the ARO, through MURI SLICE. Fatima T. Husain was partially supported by a grant from DOD W81XWH-15-2-0032. Benjamin Zimmerman was partially supported by the National Institutes of Health grant R90AT008924-08. Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research, and by the McDonnell Center for Systems Neuroscience at Washington University.

A very early preprint was published on the bioRxiv preprint server: https://doi.org/10.1101/2021.05.16.444387. A supporting information to this manuscript provides a greater detail regarding the CA pipeline as well as the toolbox of persistent homology used here. Supporting information for this article is available at https://doi.org/10.1162/netn_a_00392.

Ivan Abraham: Formal analysis; Investigation; Methodology; Software; Visualization; Writing – original draft. Somayeh Shahsavarani: Conceptualization; Validation; Writing – original draft; Writing – review & editing. Benjamin Zimmerman: Conceptualization; Validation; Writing – original draft; Writing – review & editing. Fatima Husain: Conceptualization; Validation; Writing – review & editing. Yuliy Baryshnikov: Conceptualization; Formal analysis; Methodology; Supervision; Writing – original draft; Writing – review & editing.

Yuliy Baryshnikov, Army Research Office (https://dx.doi.org/10.13039/100000183), Award ID: MURI SLICE. Fatima Husain, U.S. Department of Defense (https://dx.doi.org/10.13039/100000005), Award ID: W81XWH-15-2-0032. Benjamin Zimmerman, Foundation for the National Institutes of Health (https://dx.doi.org/10.13039/100000009), Award ID: R90AT008924-08.

The Human Connectome Project S1200 dataset is available for download from the project website with its own license agreement.

BOLD:

Acronym for blood-oxygen level-dependent signal, a measure of oxygenation levels of different regions of the brain used as a proxy for neuronal activity.

Functional connectivity:

Brain regions that show statsticially significantly correlated BOLD signals over a period of observation are said to be functionally connected.

Reparametrization invariant:

Features of time series data that do not change with arbitrary monotonic transformations of the timeline are said to be reparametrization invariant.

Oriented area:

The standard notion of area augmented with a ± sign depending on whether the orientation attached to the area’s enclosing boundary is counterclockwise or clockwise, respectiviely.

Lead matrix:

A skew-symmetric matrix whose (k, l)-th entry records the oriented area computed between signals x_k and x_l. The eigenvalues of this matrix are 0 or purely imaginary complex conjugate pairs.

Cyclic signals:

Signals that repeat in time often absent any fixed or predictable period. Periodic signals are cyclic but cyclic signals need not be periodic. Analysis of such signals is called cyclicity analysis (CA).

Chain of offsets model:

A model predicated on the assumption that the same unknown function ϕ generates different signals: x_k = α_kϕ(t − α_k), possibly in the presence of corrupting noise.

Constellations:

Each eigenvector of the lead matrix, with its complex valued entries considered as a collection of indexed points on the Argand plane, constitutes a constellation.

rs-fMRI:

Time series data obtained from fluctuations of the BOLD signal in the brain in the absence of any cognitive or mental tasks, that is, at rest.

Persistence diagram:

A scatterplot encoding the generation and annihilation of features (for the purposes of this paper, local maxima and minima) in a signal.