Hemodynamic Cortical Ripples through Cyclicity Analysis
Note: Submitted and under review.
1 Introduction
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 infra-slow (\(< 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 [1–7]. 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.
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.
1.1 Resting-state slow traveling waves
rs-fMRI: Time-series data obtained from fluctuations of the BOLD signal in the brain in the absence of any cognitive or mental tasks; i.e. at rest.
Investigating the dynamics of resting-state BOLD signals has been the target of significant research efforts. An emerging corpus of work has revealed that infra-slow spontaneous fluctuations of BOLD signals are organized not only spatially into neural networks but also temporally in an orderly fashion [8–10]. 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 [11–13]. Previously, spontaneous and stimulus-evoked traveling waves propagating across the cortex were reported using electroencephalography (EEG) and electrocorticography (ECoG) [14, 15]. 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. Infra-slow neuronal activity (\(<0.1\) Hz) was first observed in the rabbit cortex’s local field potential [16]. Later, full-band EEG was used to detect this activity in humans [17]. These infra-slow oscillations are a fundamental feature of brain function, which correlate with higher frequency oscillations and influence behavior [18–20]. In addition, it has been demonstrated that these oscillations are correlated with infra-slow hemodynamic activity and underlie changes in fMRI resting-state networks [21–23]. Notably, Chan and colleagues [24], using wide-field voltage-sensitive dye imaging, have shown that the infra-slow spontaneous neural fluctuations recapitulate cortical motifs at higher frequencies. Using wide-field optical calcium imaging in transgenic mice, Matsui and colleagues detected infra-slow global traveling waves across the cortex [11]. 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 [25], 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.
Functional connectivity (FC): Brain regions which show statsticially significantly correlated BOLD signals over a period of observation are said to be functionally connected.
1.2 Significance and innovation
The primary goal of this paper is to present the detection of traveling waves propagating across the cortex by expanding the previously presented Cyclicity Analysis (CA) pipeline [26]. To do so, we used both rest and task-based fMRI datasets from the Human Connectome Project. The voxel-level time series were aggregated to derive time series describing the hemodynamics of the ROIs, as described in the Section 6 on methods & materials.
We report three significant novel observations:
- A time-resolved temporal ordering among resting-state activity of brain regions was identified with the wave propagating along the anterior-posterior axis.
- Often, the dynamics of temporal ordering were characterized by time-localized `bursts’ rather than gradual changes.
- The direction of wave propagation during certain tasks exhibited reversals when compared with the resting state.
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).
Our work also uses novel methodological approaches. While CA was conceptualized and tested in the context of cyclic phenomena, we discovered that in 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 the reparametrization invariant data analysis, the length of the calm interval is immaterial. The key difference with the cyclic situation is the fact that at times most most of the network nodes are at rest. Exhibiting unexpected robustness, CA proved capable of detecting the presence and describing the structure of these ripples. This fact is one of the two methodological novelties of this work.
Reparametrization invariance: Features of time-series data that do not change with arbitrary monotonic transformations of the timeline are said to be reparametrization invariant.
Our second methodological innovation involves the systematic use of the oriented areas as functions of time. By design, the oriented areas (discussed in Section 6.1), or other iterated integrals, are numeric quantifiers: each of them associates a single number to a (multi-dimensional) time series. However, one can also consider these time intervals as expanding in time, making their upper bound variable, and evaluate the oriented areas on subintervals \([t_s,t]\) such that \(t_s < t \le 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 CA algorithms without any modifications to the time series extracted from Regions-Of-Interest (ROIs) - i.e., groups of voxels in the brain aggregated together.
Oriented area: The standard notion of area augmented with a \(\pm\) sign depending on whether the orientation attached to the area’s enclosing boundary is counter-clockwise or clockwise respectiviely. See Section 6.1.
2 Analytic approach
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 [27]. The applicability of the CA pipeline to analysis of fMRI traces was demonstrated in [26, 28, 29].
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, i.e. on the iterated integrals of order \(2\). These integrals, essentially, represent oriented areas encompassed by the \(2\)-dimensional 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 [30].
2.1 Cyclicity Analysis and Travelling Waves
Computationally, CA is a straightforward and efficient pipeline (we provide greater details on the algorithmic part of CA in Section 7.1). One begins with computing the oriented areas encircled by the parametric plots of the timeseries (appropriately interpolated; linearly in this study) for any pair of the time series \(x_k, x_l\) (see Figure 2). \[ A_{kl}=\frac12\int_{t_s}^{t_f} x_k dx_l -x_l dx_k, \tag{1}\]
It is clear that swapping the order of time series flips the sign of the oriented area: \(A_{kl}=-A_{lk}\). As we argue in Section 6.1, the positive oriented area \(A_{kl}>0\) is heuristically interpreted as \(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 anti-correlated signals vanishes, as well as the area spanned by signals with non-overlapping supports.)
To address pairwise leader-follower relationships synergistically, these oriented areas are arranged into the \(n\times n\) matrix \(A_{kl}\) for \(k,l=1,\ldots,n\) (which is therefore skew-symmetric). This matrix is referred to as the lead matrix (see Section 6.2 section for a precise mathematical definition & Figure 3 for a representative example), and its spectral decomposition allows one to assess the collective temporal interdependencies between the components \(\{x_k\}_1^n\).
Lead matrix: A skew-symmetric matrix whose \((k, l)\) entry records the oriented area computed between signals \(x_k\) and \(x_l\). The eigenvalues of this matrix are zero or purely imaginary complex conjugate pairs. See definition in Section 6.2.
The non-zero eigenvalues of the lead matrix consists of the pairs \(\pm i \lambda_s\), for some real \(\lambda_s\), which we assume are ordered by size: \(|\lambda_1|\geq |\lambda_2|\geq \ldots\). The eigenvectors corresponding to the conjugate eigenvalues \(\pm i \lambda_s\) are complex conjugate (and necessarily have non-vanishing 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 \(\pm i \lambda_s\) indicates a good approximation of the observations by a just few traveling waves.
2.2 Ripples
In some natural classes of cyclic, non-periodic signals - defined and referred to as the Chain of Offsets Model (COOM) in the following text - the components of the eigenvector corresponding to the eigenvalues with the largest (in absolute value) eigenvalue \(|\lambda_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.
COOM: Chain Of Ofsets Model; a model predicated on the assumption that the same unknown function \(\phi\) generates different signals \(x_k:= a_k \phi(t - \alpha_k)\) - not an uncommon situation. E.g. neuronal firing mechanisms in grey matter can be assumed to be the same whether in the frontal or parietal lobe. See Section 6.3.
Thus, we pick the eigenvalue with the largest absolute value \(|\lambda_k|\) and select one of the corresponding eigenvectors which we call the leading eigenvector. The collection of the components of the eigenvector (these components are complex numbers which 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 data (e.g. see Figure 4). 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.
Constellation: 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.
2.3 Examples
The plots below are examples of different multi-dimensional BOLD timeseries, the lead matrices obtained from them as well as the corresponding constellation from the leading eigenvector.
Under the reasonably mild 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 traveling wave. In Section 7.2 and Section 7.3 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, in some natural metric, parts of the constellation and focus on those.
2.4 Oriented areas as functions of time
While the coefficients of the lead matrix describes the overall, aggregated leader-follower relationship between a pair of network nodes (ROIs), the behavior of the oriented area as a function of time can be indicative of the time-heterogeneous nature of this relationship. To provide with a specific, baseline scenario of such behavior we consider the standard leader-follower pair of harmonic series, cosine and sine. In this case, if \(x=x_0+a\cos(t), y=y_0+b\sin(t), a,b>0\), one obtains
\[ A(t)=\frac12\int_0^t x dy-y dx= ab t+A\sin(t+\phi) \]
where \(A^2=x_0^2b^2+y_0^2a^2,\tan\phi=\dfrac{bx_0}{ay_0}\). In particular, the growth of the oriented areas will be smooth, with small oscillation if the averaged values of the time series, \(x_0, y_0\) are close to zero.
3 Results
As described above (and in greater detail in Section 6), we started by generating the timeseries describing the BOLD activities in 68 ROIs (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 recordings of duration \(t_f=864 s\). The time series were then centered (setting their time averages to zero) and normalized by their quadratic variation, i.e sum of squared differences. Oriented areas \(A_{kl}(\cdot)\) were then computed using a discretized version of Equation 1.
For the full lead matrix \(A\), the spectral decomposition was determined and the eigenvector \(v_*\) corresponding to the largest in absolute value eigenvalue \(i\lambda\) chosen. (There is an ambiguity as to which of the two complex conjugate eigenvectors corresponding to the opposite eigenvalues \(\pm i\lambda\) to chose. However, selecting a pronounced leader-follower pair of components \(x_k, x_l\) 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 a constellation; example of which is shown in Figure 4.
3.1 Dominant ROIs
In the case of perfect 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 [28]. 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.
In this study, we identified the CE or, rather, the positive-definite quadratic form \[ q(z)=ax^2+2bxy+cy^2, z=x+iy \] on the complex plane (for which the CE is the unit level set \(\{q(z)=1\}\)) using least square regression. The resulting norm \(|z|_q:=q(z)^{1/2}\) allows us to order the ROIs (points in the constellation) according to their distance defined in terms of \(q\): the point of the constellation with the largest norm has rank one, and so on - see Figure 5 for a representative example.
Naturally, these quadratic norms \(|\cdot|_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 (that is, of high typical \(|\cdot|_q\)-distance to the origin) are expected to be exposed to the presumed excitation ripples through the cortex, while those of higher average rank 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 6 summarizes the average rankings for each ROI corresponding to the first eigenvalue.
From here on we focus primarily on the leading eigenvalues (with the largest absolute value, i.e. \(\lambda_{1,2}\) in our convention), but one can extend the same procedures to the next largest in magnitude eigenvalue as well (see Figure 7) and so on.
3.2 Ordering the ROIs phases
Using the average ranking chart above we can choose a certain cutoff, to limit the number of ROIs affected by the leading ripple (i.e., corresponding to the leading eigenvalue), for further analysis.
Once such a subset is chosen, we can estimate the cyclic ordering between the ROIs by examining the arguments of the complex valued entries in the eigenvector corresponding to them. As noted in Section 2.2, CA is able to extract an approximate cyclic ordering if the nodes are excited by similar waveforms delayed by a collection of phase-shifts. It is cyclic because the ordering of the repeating sequence of events (such as the excitation of a node) has no intrinsic starting point: only the relative positions of the phases are important. Still, quite often there is a natural clustering, of the points in the constellation, corresponding intuitively, to a sequence of waves, coming through the network interspersed by the absence of the signal.
This phase shift extraction is especially intuitive when the constellation is clustered in one half-plane of the Argand diagram as in Figure 4. 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 the ordering of arguments of the (complex) components corresponding to the ROIs, normalized to vary from \(-\pi\) to \(\pi\). This is shown in Figure 8 which is obtained from Figure 5 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 level 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 particular position in the ordering of the phase shifts derived from the constellation. This is depicted in the violin plots below for the leading eigenvalue for \(k=20, 24, 28\) and \(54\) which correspond to natural cut-offs at L_supramarginal
, L_rostralmiddlefrontal
, L_superiorfrontal
and L_parsorbitalis
(see Figure 6) as well as all 68 ROIs. The ordering of the ROIs is by the average positions in the phase shifts order.
Keeping all ROIs, i.e. \(k=68\) results in the following visualization for wave propogation.
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 passage. Furthermore, changing the cut-off thresholds does not qualitatively change the observed pattern.
If one attempts to repeat the analysis of the excitations’ passage using the higher eigenvalues (i.e., \(\pm i\lambda_s, s\geq 2\)) one encounters loss of directionality: the violin plots become disordered, and lack any discernible structure, see Figure 14 and Figure 15 with the violin plots for the second and third largest eigenvalues.
3.3 Pairwise timeseries analysis
Given a subcollection of ROIs chosen as a dominant set, one can examine in further detail the time evolution of their leader-follower dynamics as previously mentioned. This results in the production of \(\binom{k}{2}\) new timeseries where \(k\) is the chosen cut-off applicable to Figure 6 (see Equation 2 for the precise equation). Within such a time series, a consistent increasing or decreasing trend over the period of observation shows that there is an average leader (or follower) relationship between the pair of time series - i.e. activity in one region precedes (or lags) activity in the other region. While the overall trend or increment is succinct indicator of this general relationship, the intermediate dynamics of the oriented area signal can provide insights into the nature between the pair of brain regions that constitute it. We noted primarily three kinds of dynamics that oriented area signals exhibited from the dataset.
- A near constant zero dynamics
- An overall increment or decrement punctuated by role-reversals in the leader-follower relationship
- An steady, nearly monotonic overall increment or decrement
We illustrate the three kinds of dynamics in the figures below:
We observed the first kind of dynamics mostly between bilateral region pairs while the last was a hallmark of the resting state dataset. More interestingly, the second kind of dynamics were observed in a larger fraction of the task paradigm 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.
The next two subsections deal with the detailed analysis of these observations, where we devise some tools to quantify these phenomena.
3.4 Reversals under the task paradigm
To further examine the phenomenon of dynamic reversals observed in the oriented area signals we accumulated over all pairs of ROIs the timestamps at which a change in the direction of the time series 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 Section 7.4. One does this by coupling the local minima and maxima in a well-defined way. The heights \(b,d\) of thus coupled minimum and maximum are thus referred as birth and death of a feature, and are visualized as points on the \(b,d\) plane, as shown in Figure 19. 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 [31] for more details).
In such a figure instances far from the diagonal line are said to be persistent features in the data (see [32] for more details). By choosing to keep some \(p\) most persistent features across the datasets, we can visualize the collection of their timestamps acrross all oriented area time-series using a carpet plot as follows:
In the above panel it is clear that reversal in the oriented area dynamics are a clear indicator of the difference between the resting state and task paradigm datasets. One prominent feature of the social task plots is the regular patterns occuring approximately every 52 frames (i.e., every \(.72s\times 52=37.4s\)). The documentation [33] reveals that the intervals between social task segments was \(38s\).
A similar analysis of the resting state sessions reveals that pairwise activities, as indicated by the above outlined analysis, are, indeed at rest, see Figure 21.
3.5 Burstiness in time series
To analyze burstiness in the time series, we extend 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 the oriented areas as functions of time was that activity between some regions tended to be much more consistent than others. These are exhibited by 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”. The below set of figures illustrate this phenomenon. Figure 22 shows a pairwise oriented area integral dominated by the jumps from deep local minima. Figure 23 however shows another pair of brain regions where the total accumulated oriented area is similar to that of the first display, but gets there in a more gradual fashion.
We have marked on either plot the ten longest persistent features observed in the time series calculated according 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 ten longest bars in the persistent diagram.
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 Figure 24, 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. (We 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.)
4 Discussion
We believe that Cyclicity Analysis holds the potential to provide new insights into dynamics of FC, building an improved understanding of brain networks and the interplay among them.
Our results provide supporting evidence for the notion of cortical waves propagating along the cortex between primary visual and transmodal areas of the brain. In addition, we found that the oriented areas between pairs of regions reflected in the lead matrix were driven mostly by short periods of sustained lead-lag relationship between the time series, suggesting a time dependent flow of information underlying the dynamic aspect of resting states hemodynamic measures.
Although the outcomes derived from both lagged correlations and CA support the idea of slow cortical ripples, they have different underlying mathematical apparatuses, assumptions, and levels of granularity. The lagged correlation method infers lag threads by deriving singular vectors of the time-delay matrix, whereas CA recovers inherent ordering among BOLD time series through eigenvectors of a lead matrix (a representation of the strength of temporal ordering between pairs of regions, see Section 2.1). Lagged correlation methods rely on interpolation and windowing to capture the dynamics of FC, which are sensitive to time delay estimation methods, autocorrelation [29], sampling variability [34], and parameters such as window length and window shift [35]. In contrast, CA offers a more robust approach with finer granularity resolution, and without assumptions regarding stationarity, state duration, or state transition, all following the principled reliance on reparametrization invariant tools.
Although cortical waves have been noted at much higher frequencies, detecting them in functional MRI data at sub-Hz frequencies 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 [9, 14, 36–40]. The temporal duration of these more typical cortical waves tends to be on the order of tens to hundreds of milliseconds rather than the sub-Hz frequencies which BOLD signals are able to resolve. Nevertheless, over the past decade, a growing body of studies has observed BOLD signal flow across cortex [12, 41–44] .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 non-linear reparameterizations of the timeline, a feature lacking in the correlation-based methods.
One of the more striking patterns our analysis uncovered was that the overall lead-lag relationships between pairs of regions in the brain were often the result of short bursts of strong interactions, rather than long consistent stretches of moderate ones. This should be compared to the recent paper by [45], where the authors investigated the contributions of moment-to-moment BOLD activity to the overall pattern of functional connectivity. 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 corresponding to fluctuations in the level of default mode and attentional control network activity, which are often viewed as in opposition to each other [46].
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 [47, 48]. According to this framework, a topography of connectivity patterns is reflected in the spatial distances between so-called “higher” areas of the brain, where more multi-modal, 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, rather than as independent networks. The functions associated with various cortical networks were correlated with the level of the hierarchy of sensory or motor processing. Also, [12] 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, that should be investigated in future studies.
In the context of the resting state scans, we see ordering in the direction that moves in the reverse hierarchy, through association areas toward the visual regions (see Figure 13). The precentral gyrus activity leads and the lateral occipital cortex follows. 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 20), we observed a strikingly different pattern. This task consisted of 5 blocks, each consisting of a period of watching a video and then judging whether or not the objects in the video appear to perceive the other object’s feelings and thoughts. In this task, we still observed strong directional constraints, but which 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, that may correspond with periods of watching the videos and periods of making judgements of the objects’ intentions within the video. The direction of the leader-follower relationships are also much burstier compared to 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 21 vs Figure 20). Here the participants were presented with three-second cues to perform twelve motions. There are 13 motor blocks and three 15-second 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 and motor movement events 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 activity belonging to hierarchical visual processing and motor processing.
Our analysis identified a temporal ordering that primarily propagated along the anterior-posterior axis of the brain. 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 \(|\lambda_1/\lambda_3|\) 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 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 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 the anterior-posterior axis. 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 (MEG). 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 whole-brain function.
5 Materials
From the Human Connectome Project 1200 Release [33, 49], we considered 889 de-noised minimally preprocessed participants who completed all the structural, resting state, and task fMRI sessions using a customized Siemens Skyra 3T 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 regions of interest (ROI) information in the fMRI parcellations available for download.
ROI | ROI full name | ROI | ROI full name |
---|---|---|---|
1 | Banks of superior temporal sulcus | 18 | Pars orbitalis cortex |
2 | Caudal anterior cingulate cortex | 19 | Pars triangularis cortex |
3 | Caudal middle frontal cortex | 20 | Pericalcarine cortex |
4 | Cuneus | 21 | Postcentral cortex |
5 | Entorhinal cortex | 22 | Posterior cingulate cortex |
6 | Fusiform cortex | 23 | Precentral cortex |
7 | Inferior parietal cortex | 24 | Precuneus |
8 | Inferior temporal cortex | 25 | Rostral anterior cingulate cortex |
9 | Isthmus cingulate cortex | 26 | Rostral middle frontal cortex |
10 | Lateral occipital cortex | 27 | Superior frontal cortex |
11 | Lateral orbitofrontal cortex | 28 | Superior parietal cortex |
12 | Lingual cortex | 29 | Superior temporal cortex |
13 | Medial orbitofrontal cortex | 30 | Supramarginal cortex |
14 | Middle temporal cortex | 31 | Frontal pole |
15 | Parahippocampal cortex | 32 | Temporal pole |
16 | Paracentral cortex | 33 | Traverse temporal cortex |
17 | Pars opercularis cortex | 34 | Insula |
The ROIs extracted were 34 in number (bilateral, therefore total N=68) and are listed in the 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(TR) of 720 ms, for the resting state fMRI scans, this resulted in a set of arrays each with elements of dimension 68 \(\times\) 1200 that were fed into the CA pipeline. (Here 1200 is the number of the samples; the total session durations therefore being 14.4 minutes). Further details of the scan protocol are available here.
The HCP dataset includes fMRI scans recorded under different task paradigms. Such scans recorded during progression of certain motor and cognitve tasks were also included in part of our analysis. The two task protocols are as follows - for the motor task [50]:
Participants are presented with visual cues that ask them to tap their left or right fingers, squeeze their left or right toes, or move their tongue to map motor areas. Each block of a movement type lasts 12 s (10 movements), and is preceded by a 3 s cue. In each of the two runs, there are 13 blocks, with 2 of tongue movements, 4 of hand movements (2 right and 2 left), 4 of foot movements (2 right and 2 left) and three 15 s fixation blocks per run.}
and for the social task [50]:
Participants are presented with short video clips (20 s) of objects (squares, circles, triangles) either interacting in some way, or moving randomly. These videos were developed by either Castelli and colleagues (Castelli et al., 2000) or Martin and colleagues (Wheatley et al., 2007). After each video clip, participants chose between 3 possibilities: whether the objects had a social interaction (an interaction that appears as if the shapes are taking into account each others feelings and thoughts),Not Sure, or No interaction (i.e., there is no obvious interaction between the shapes and the movement appears random). Each of the two task runs has 5 video blocks (2 Mental and 3 Random in one run, 3 Mental and 2 Random in the other run) and 5 fixation blocks (15 s each).}
6 Methods
Recall that CA pipeline builds on the iterated integrals of the second order [26].
6.1 Oriented areas
Such iterated integrals, not expressed in terms of the iterated integrals of 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 25(a)) and its time shifted copy, \(y(t)\) (orange curve). We interpret the interrelationship between \(x\) and \(y\) as a leader-follower one: \(x\) leads; \(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 region \(R\) of the plane of positive area (as shown in Figure 25(b)).
Consider a closed parametric curve \((x(s), y(s))\in \mathbb{R}^2, s\in [t_1,t_2]\) with piece-wise 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. Alternatively, the algebraic area can be calculated using the representation
\[ \operatorname{area}\left(R\right) = \dfrac{1}{2} \oint \limits _C xdy - y dx= \dfrac{1}{2} \int_{t_1}^{t_2} x\dot{y} - y \dot{x} \; ds, \tag{2}\]
implying that the oriented area is an iterated integral of order 2 (here \(C\) is the curve, serving as the contour of integration, oriented by the time and \(R\) is the region it encloses). Note that this area is signed; in particular, reversing the orientation of the curve or interchanging the variables flips the sign of the oriented area.
Note that the quantity in Equation 2 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 two signals do not overlap in time, their oriented area measure is also zero.
6.2 Lead matrices
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 \(\mathbb{R}^n\), we arrange the oriented areas into the square \(n\times n\) matrix (here \(n\) is the number of time series observed). This matrix, whose \(kl\)-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 [28]. We remark that this matrix is obviously skew-symmetric, implying in particular, that its eigenvalues form pairs of purely imaginary numbers, \(\pm i\lambda,\lambda\in\mathbb{R}\), and the corresponding eigenvectors form complex-conjugated pairs with necessarily complex-valued components.
6.3 Chain-Of-Offsets-Model
Consider now the situation where the coordinates \(x_j, j=1,\ldots,n\) of the trajectory correspond to the same function (which we will interpret here as the function on the internal clock space, a circle) \(\phi:S^1\to\mathbb{R}\), just offset by a different phase. In other words,
\[ x_k(t)=a_k\phi(t-\alpha_k); \quad \alpha_j\in S^1, k=1,\ldots,n \tag{3}\]
This is the Chain-Of-Offsets-Model (COOM) alluded to in Section 2. In this case, the lead matrix can be readily computed in terms of the Fourier coefficients \(c_m\), of \(\phi\): it is given by the summation
\[ A^\phi_{k,l}=2\pi a_ka_l \sum_{m\geq 1} m|c_m|^2 \sin(m(\alpha_k-\alpha_l)) \tag{4}\]
In particular, the lead matrix given by Equation 4 decomposes into the sum of rank two skew symmetric matrices with entries
\[ A^{(m)}_{k,l}=m|c_m|^2a_ka_l\sin\left(m(\alpha_k-\alpha_l)\right) \]
If one of the coefficients in the Fourier series for \(\phi\) dominates, the skew symmetric matrix \(A^\phi\) is well approximated (in Frobenius norm) by the rank two matrix \(A^{(m)}\) with each \((k,l)\)-entry given by:
\[ \left(|c_m|^2a_ka_l\sin(\alpha_k-\alpha_l)\right)_{kl}=|c_m|^2\left(v_k u_l -v_l u_k\right)_{kl} \tag{5}\]
where \(u_k=a_k\cos(\alpha_k), v_k=a_k\sin(\alpha_k)\).
In general, if a skew-symmetric operator of rank 2 is represented as \(Q=u\otimes v'-v\otimes u'\) (where \(u,v\) are linearly independent, and \(v'\) denotes the conjugate vector to \(v\)), then one can see that \[w=-e^{-i\theta} u/|u|+v/|v|\] is an eigenvector of \(Q\) with the eigenvalue \(i\sin\theta |u||v|\) (here \(\theta\) is the angle between \(u,v\)).
This implies that the real and imaginary components of the eigenvector \(w_k=p_k+iq_k\) are obtained from the real and imaginary components \(u_k\) and \(v_k\) defining the matrix \(A^{(m)}\) via the linear transformation of the plane
\[ \left( \begin{array}{c} p_k\\ q_k \end{array}\right)= \left( \begin{array}{c} Lu_k+Mv_k\\ Nu_k \end{array}\right) \]
where \(L=-\cos\theta/|u|, M=1/|v|\) and \(N=\sin\theta/|u|\).
Therefore, if \(p_k+iq_k\) are the components of the eigenvector corresponding to the rank \(2\) matrix stemming from a purely harmonic COOM with the offsets \(\alpha_k, k=1,\ldots,n\), they are just obtained from the complex numbers \(a_k\exp \left( {i\alpha_k} \right), k=1,\ldots,n\) by a linear transformation of the complex plane, and therefore the cyclic order defined by the arguments of these components is the same as the cyclic order of the collection of points on the unit circle: \[ \{\cos(\alpha_k)+i \sin(\alpha_k)\}_{k=1,\ldots,n} \]
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 \(\mathbf{x(t)}\).
7 Appendix
7.1 Theoretical considerations
The problem of finding reparameterization invariants (RI) of trajectories \[\mathbf{x}(t)=(x_1(t),x_2(t),\ldots,x_d(t))\] in a multidimensional space \(\mathbb{R}^d\) was addressed by [27], where the author established that essentially, all RI quantifiers are given by the iterated integrals, a family of functionals of trajectories that can be defined iteratively and computed efficiently.
Chen’s iterated integrals form a countable family: the degree of an iterated integral is one more than the number of iterative integrations needed to compute it. For each degree, there are several distinct (functionally independent) iterated integrals, and their number, and the overall structure are described by combinatorial apparatus relevant in several disjoint areas of mathematics. From empirical perspective, however, the iterated integrals of higher orders become excessively prone to the noise, and we still don’t have convincing examples of their practical usefulness.
7.2 Model example: signal propagation in networks
To illustrate how the CA pipeline works, we introduce a model example, where the signals are propagating through a network by gossiping. This material is dedicated to the description of a simulated model, with which we illustrate the computational pipeline of tools we applied to the fMRI data. Conceptually, this model, of the waves propagating in a network, broadcast from an unknown source, is a highly stylized caricature of the cortex waves. While we presented the results on the (tentative) recovery of those waves from the fMRI readings in the main text, it is enlightening to see what the situation is in a model example, where the ground truth is known.
7.2.1 Networks & signals
Consider a weighted undirected graph \(G\) with \(n\) vertices, interpreted as a network, with edge lengths proportional to internode communication delays. Consider the following broadcasting protocol: if a node emits a signal, it goes to its neighbors, who will receive it after the edge-specific delay, and immediately rebroadcast it. If the signal was received by the node earlier, it ignores it. This implies that each signal reaches any given node along the shortest path connecting it to the source.
Such peer-to-peer propagation models are often referred to as gossip, epidemic or first passage percolation networks [51, 52], and are relevant in studies of social networks, peer-to-peer networks, distributed resource location, etc. [53–55]. Assuming such a model, one can easily recover the structure of the underlying network, if the shortest path lengths between the pairs of the nodes are known.
Suppose a signal \(s_i(t)= f(t)\) propagates from a source node \(i\). If the underlying graph \(G\) is a connected one, then a node \(j \neq i\) of the network will eventually receive the signal from \(i\) as a delayed shapeform so that,
\[ s_j= f(t-d(i,j)) \tag{6}\]
where \(d(i,j)\) is the shortest path distance between nodes \(i\) and \(j\). For connected graphs, this distance is well-defined and finite even if there is no direct connection between nodes \(i\) and \(j\). Quite often, one is only interested only in the ordinal information about the edgelengths (for example, the minimal spanning trees depend only on the ordering of the edge lengths).
Can this ordinal information be reconstructed using only the observations of the (perhaps, noisy) signals at the nodes?
7.2.2 Cyclicity Analysis
We will apply the CA pipeline to recover the network structure from the observations of the signals \(s_k, k=1, \dots, n\), as it manifestly fits the Chain Of Offsets Model (see Section 6 section of the paper). For each of the \(\binom{n}{2}\) pairs of signals given by Equation 6 observed at the nodes of the network, we form the oriented area of the corresponding \(2\)-dimensional projection as,
\[ \operatorname{area} = \dfrac{1}{2} \oint \limits _C xdy - y dx= \dfrac{1}{2} \int_{t_1}^{t_2} x(s)\dot{y}(s) - y(s) \dot{x}(s) ds, \tag{7}\]
with \(x, y = s_i, s_j\) for \(i, j = 1, 2, \dots, n\) and \(C\) being the parametric contour linearly closed up. As explained in the article, the orientation is indicated by attaching a sign to the computed area value, with a positive sign corresponding to counter-clockwise integration; interpreted as \(j(t)\) following \(i(t)\) in time. The spectral decomposition of square matrix \(A\) whose entries are given by Equation 7, i.e., the lead matrix, allows us to determine the relative distances of the nodes of the network to the seed node, where the signal originates.
Indeed, as we argued previously, the order of the arguments of the components of the (complex) eigenvector corresponding to the leading (in absolute value) eigenvalue of the lead matrix will reflect the cyclic order of the propagating wave, if the rank two skew-symmetric lead matrix corresponding to the purely harmonic signal propagating through the network is a good enough approximation of the sampled lead matrix \(A\). In this case, spectral analysis of the lead matrix recovers the lag-structure between the source node and every other node.
Eigenvalue ratio: From Section 6.3, recall the rank-2 decomposition; in particular for such matrices, the \(|\lambda_1\left( A \right)|/|\lambda_3\left( A \right)|\) is undefined because \(\lambda_3 = 0\). Thus, larger \(|\lambda_1\left( A \right)|/|\lambda_3\left( A \right)|\) ratios indicate better approximations.
Note that one measure of how well the assumptions and results hold is obtained by computing the eigenvalue ratio \(|\lambda_1\left( A \right)|/|\lambda_3\left( A \right)|\) of the matrix \(A\), with larger ratios leading to better results. Repeating this analysis for various source nodes would then enable one to reconstruct the structure of the network.
7.3 Simulation study
In Figure 26 we introduce a network on \(n=12\) nodes. The mathematical graph of Figure 26(a) is constructed by starting with the \(C_3\) graph and adding nodes with \(k=3\) edges at a time. The edges are attached to vertices at random following a distribution proportional to the vertex degree . The shown matrices represent the internode distances, both as given by the edge lengths in Figure 26(b), and by the shortest path distances between pairs of nodes in Figure Figure 26(c). Note that the shortest distance using relays can be shorter than the direct link.
To generate the signal in Figure 26(d) emitted by the seed note (node \(6\) in this case), we use a convex combination of randomly modulated Gaussian functions, randomly displaced,
\[ f(t) = \sum \limits _{k=1}^{s} r_k g(t,k), \quad g\left(t,k\right) = \exp \left(-\left(t-k\pi + c\mathbf{R}_k \right)^2 \right). \tag{8}\]
Here the displacements \(\mathbf{R}_k\) and amplitudes \(r_k\) were drawn uniformly from the appropriate intervals. The signal propagates through the network according to the equation Equation 6, to which we add a small noise, realized as the scaled Brownian motion, independently at each node.
Details of the model example are as follows:
- Panel (a) shows the graph \(G\) with \(N=\) 12 nodes and source node highlighted. The edge thickness are in proportion to the internode signal delay.
- In panel (b) the adjacency matrix with random edge weights is visualized. Zeros indicate the absence of a direct link between the nodes.
- Panel (c) visualizes the shortest-path length matrix. Note that the connectivity of the graph results in all non-diagonal elements being nonzero.
- Finally panel (d) shows the randomly generated signal at the source node as per Equation 8.
This grid shows the slightly delayed and noisy waveforms that is recieved at each node. The source node itself is assumed to be subject to some process noise and each recieving node also has independent sensor noise entering additively.
This grid shows the results of CA for the model example.
- In panel (a) we show the skew-symmetric lead matrix.
- In pabel (b) we have the (absolute values of) the eigenvalues of the lead matrix. The rank \(2\) matrix corresponding to the first two (complex-conjugated) eigenvectors approximates the lead matrix well, as the first pair of eigenvalues dominates:
\[ \left| \lambda_1 \right| / \left| \lambda_3 \right|= 5.45 \]
- In panel (c) the ``constellation” of the components of the leading eigenvector winds around the origin; the circular order of their arguments (phases) indicates the propagation of the gossip wave through the network. More precisely, as shown …
- In panel (d) the scatter plot of these arguments against the known distances from the source node shows a near perfect linear dependence, implying that lag-structure is encoded with eigenvector.
Typically the distances between the nodes and the source will be unknown, however repeated analysis with multiple sources and careful consideration of obtained complex arguments from the eigenvectors enables inference of the underlying pairwise relative distances.
The results of the cyclicity processing pipeline for the time series generated in the example of Figure 26 are presented in Figure 28. Figure 28(a) shows the sampled lead matrix. The interpretation of the entries is quite intuitive; the vanishing oriented area between a pair of nodes is consistent with the fact that the network distance from the seed node to either of them is approximately the same.
Moreover, the first pair of the conjugate eigenvalues dominates in the spectrum of the lead matrix, indicative of the high resolution of the analysis as shown in Figure 28(b). The eigenvector corresponding to (one of) the leading eigenvectors has components shown in Figure Figure 28(c). The arguments of these components, that is the angles formed by the rays pointing towards them and the ray of positive real numbers are shown on the display Figure 28(d) as the scatter plot, against the shortest distance to the seed node. The strong, essentially linear dependence between these phases and the distances to the seed shows that cyclicity can be used to detect the latter from the former, in a reparameterization invariant way.
7.4 A dash of persistent homology
We will describe here a small subset of the persistent homology theory, dealing with the basic case of \(0\)-dimensional homology, - essentially, tracking the connecting components of the sublevel sets for a function \(f\).
Consider a (continuous) function \(f\) of a topological space \(X\) (in our situation, the topological space is just a subinterval \(X=[t_s,t_f]\) of the real line, but the definitions are valid in the generality we formulate them).
Consider a point \(t_*\in X\) which is a _local minimum of \(f\): this means in some small vicinity \(N_\epsilon(t_*)\) (in our situation, again, the vicinity is just a subinterval \((t_*-\epsilon,t_*+\epsilon)\)), the function does not dip below the level \(f(t_*)\): \(f(t)\geq f(t_*), t\in N_\epsilon(t_*)\), and any path from \(t_*\) to outside of \(N_\epsilon(t_*)\) has to pass through a value that is strictly higher than \(f(t_*)\).
Assume that the global minimum of \(f\) lies below \(f(t_*)\). In this case, there is a path connecting \(t_*\) to a point where \(f\) attains a value below \(f(t_*)\). What is the minimum of the maximal elevation \(f(\cdot)\) along all such paths? It is easy to see that it is achieved at some point \(t^*\) with \(f(t^*)>f(t_*)\) which is a saddle in the dimensions \(>1\), or a local maximum in the situation of \(d=1\), which is of primary interest to us. This point is well-defined in generic case (where the critical values, - the values of the function at its critical points, such as minima, saddles etc) are all different.
In this situation one says that the critical values \(f(t_*)\) and \(f(t^*)\) are coupled, and the pair \((b,d)=(f(t_*),f(t^*))\) forms a \(0\)- persistence bar. The endpoints of this bar are referred as the birth and _death values \(b<d\), represented often as a point on the \((b,d)\) plane. The collection of these points is referred to as the (\(0\)-dimensional) persistence diagram.
The global minimum is, by convention, always coupled to \(+\infty\).
Informally, it is useful to think of the process of coupling of local minima \(t_*\) with their counterparts \(t^*\) as following: imagine \(f\) as the landscape that one fills with water, starting with a basin at \(t_*\). Eventually, the water spills out to other basins, - the elevation at which it can be connecting to some {} basin is coupled to the height of the local minimum \(t_*\), see below.
Longer bars correspond to larger scale features of the function \(f\), the intuition we use in Section 3.4.