A Novel Spatial-Spectra Dynamics-Based Ranking Model for Sorting Time-Varying Functional Networks from Single Subject FMRI Data

. Accumulating evidence suggests that the brain state has time-varying transitions, potentially implying that the brain functional networks (BFNs) have spatial variability and power-spectra dynamics over time. Recently, ICA-based BFNs tracking models, i.e., SliTICA, real-time ICA, Quasi-GICA, etc, have been gained wide attention. However, how to distinguish the neurobiological BFNs from those representing noise and artifacts is not trivial in tracking process due to the random order of components generated by ICA. In this study, combining with our previous BFNs tracking model, i.e., Quasi-GICA, we proposed a novel spatial-spectra dynamics-based ranking method for sorting time-varying BFNs, called weighted BFNs ranking, which was based on the dynamical properties in both spatial and spectral domains of each BFN. This proposed weighted BFNs ranking model mainly consisted of two steps: first, the dynamic spatial reproducibility (DSR) and dynamic fraction of amplitude low-frequency fluctuations (DFALFF) for each BFN were calculated; then a weighted coeffi-cients-based ranking strategy for merging the DSR and DFALFF of each BFN was proposed, to make the meaningful dynamic BFNs rank ahead. We showed the effective results by this ranking model on the simulated and real data, suggesting that the meaningful dynamical BFNs with both strong properties of DSR and DFALFF across the tracking process were ranked at the top.


Introduction
Blood oxygen level dependent (BOLD) functional magnetic resonance imaging (fMRI) is powerful modality to discover functional connectivity (FC) among discrete brain regions.Recently, a growing number of reports have suggested that the FC un-der rest or task condition is not static but exhibits complex spatiotemporal dynamics [1].For example, besides the well-known temporal variability of FC [2][3][4], the brain activation regions also show considerable spatial variability and dynamical power spectrum over time [5][6], which has potentials for diagnosing brain diseases [7].
Independent component analysis (ICA), which could extract the brain functional networks (BFNs) from fMRI data at the single subject level or group level under rest or task conditions [1][2][3][8][9][10], has turned out to be a promising tool to decode the brain activity.With respect to tracking the dynamic BFNs at the single subject level, some variants of ICA such as SliTICA, real-time ICA and Quasi-GICA, have been proposed to capture the spatial dynamics over time [5,[11][12].However, after the BFNs identification using ICA-based tracking methods, how to distinguish the neurobiological and meaningful time-varying BFNs from noise and artifacts is not trivial due to the following reasons.To the beginning, most ICA algorithms do not automatically rank BFNs, which leads to manually examine all BFNs one by one.Further, the BFNs at single subject level always have limited information in contrast to ones from multiple subjects, which leads to that sorting BFNs at single subject level is more challenging task.In addition, to our knowledge, most of BFNs ranking methods such as percent variance ranking [2], power spectrum ranking [13], support vector machine (SVM) based ranking [14], multiple ICA estimations based ranking [15][16] and MMC ranking [17], are designed for sorting static BFNs, rather than dynamic BFNs.Thus, in this study, a novel spatial-spectra dynamics-based ranking method was proposed to sort the time-varying BFNs, taking advantage of the dynamic features in both spatial and spectral domains of each BFN.

2
Theory and Methods

Brief Review of Quasi-GICA
It is well-known that the BFNs identification has always been modulated as a blind source separation (BSS) problem, assuming that there exists functional integration of activity in multiple macroscopic loci [2-3, 8-10, 18].This BSS problem aims at retrieving the underlying sources, namely BFNs, denoted in vector notation as S having Q rows of underlying sources i s with size of Q M × , from the observed mixture X having P time points with size of P M × , which can be written as: where each column of A is called time course (TC) with the size P Q × , the row vector i s and i x are with the same size 1 M × .Founding on formula (1), the dynamical BFNs identification in Quasi-GICA [5] has the following four steps: firstly, the sliding window technique [12] is applied to generate the sliding window subset i , namely , where the window length is equal to L , and ; secondly, the data compression for each sliding window subset ˆi X is formulated as , the compressed data of X can be formulated as , , , , ; thirdly, ICA estimation is implemented on Y to obtain BFNs, and formulated as follows: ˆ( ) ) pinv denotes the pseudo-inverse operation, M is a ( ) and Ŝ is a N-by-M matrix with each row representing a BFN.To obtain the TCs (i.e., Â ), the multivariate linear regression is applied to the mixture X , depending on the separated Ŝ , i.e., ˆ( ) (4) Finally, to obtain the dynamical BFNs and dynamical TCs from each sliding window subset (i.e., ˆi X ), the dual regression is applied, and sequentially expressed as: ˆ( ) ˆˆ( )

Weighted BFNs Ranking Model
The random order of BFNs generated by ICA has negative effect on seeking biologically meaningful BFNs, which brings extra burden for manual selection of BFNs.To our knowledge, the meaningful BFNs firstly should have relatively high spatial reproducibility across different sessions and/or subsampled datasets [17,19].Also, the TCs corresponding to meaningful BFNs exhibit strong low-frequency fluctuations within the frequency band [0.01Hz, 0.1Hz] [6,8,20].Thus, by assuming that the BFNs and TCs have different properties of dynamic spatial reproducibility (DSR) and dynamic fraction of amplitude low-frequency fluctuations (DFALFF), respectively, a novel weighted BFNs ranking model is proposed, where the dynamical BFNs with both high values of DSR and DFALFF will be ranked in front of those who have relatively low values.The DSR for the spatial map as to each BFN can be defined as follows: , , where ˆk i S , ˆkj S and ˆk S denote the spatial maps of the kth BFNs from the sliding window subset ˆi X , ˆj X and X , respectively.The 1 k DSR in equation ( 7) emphasizes the mutually spatial reproducibility among the dynamical BFNs from the different slidingwindows, while 2 k DSR in equation ( 8) emphasizes the spatial reproducibility between the dynamical BFNs and the corresponding BFNs from the whole sliding-windows.
To calculate DFALFF values, the TCs are first extracted from Â and ˆi ), and then manipulated by Fourier transformation according to equation ( 9) to obtain the power spectra F and ˆi F , respectively.
( ) For a continuous frequency band [a, b], the energy at the specific interval is calculated as .
According to the low-frequency fluctuations within band [0.01Hz, 0.1Hz] of BFNs [6,8,20], the dynamic FALFF values as to the power spectra F and ˆi F for the kth BFN can be formulated as follows: The where . In this equation, each 3 Experimental Tests

Experimental Datasets
(1) Simulation Dataset One subject dataset was generated by SimTB toolbox, with V = 148×148 voxels, 12 spatial sources and 120 time points at time of repetition (TR) = 2s.Two sources shared the task-related block modulation in addition to having unique fluctuations.Activation for the other ten sources was simulated based on solely unique hemodynamic fluctuations with no task-related variation.Additive noise was also included to reach a specified contrast-to-noise ratio of 1.0.This dataset was also used in the studies [5,[9][10].
(2) Visual Task Dataset Three healthy subjects took part in a visual task.The visual stimulus paradigm was (OFF-ON) × 3-OFF in 20-second blocks.The visual stimulus was a radial blue/yellow checkerboard, reversing at 7 Hz, corresponding to the "ON" state.At the "OFF" state, the participants were required to focus on the cross at the center of the screen.The BOLD fMRI dataset was acquired on a Philips 3.0 Tesla scanner with a single-shot SENSE gradient echo EPI with TR of 2.0s.There were 40 slices providing the whole-brain coverage, with a SENSE acceleration factor of 3.0 and scan resolution of 80×80.The in-plane resolution was 3mm×3mm.The slice thickness and slice gap were 3 mm and 1 mm, respectively.

Data Processing
For the simulated data, no preprocessing step was involved; as to the real data, the standard preprocessing procedure was performed in SPM8, including slice-timing, motion correction, spatial normalization and spatial smoothing with FWHM kernel equal to 8mm.The MRIcroN software was used to display BFNs.
The Laplace approximation [21] was used to estimate the component number in Quasi-GICA [5].The sliding window length L was empirically set to 70 and 40 time points for the simulation data and visual task data, respectively.For the visualization of BFNs from the real data, the z-score normalization was performed, and then the threshed maps was obtained with thresh value set to 2.0.

4
Results and Analysis

Results of Dynamic Low-frequency Fluctuations
Fig. 1A showed the different mean FALFF values for the twelve components estimated by Quasi-GICA across the sliding windows on the simulated data, implying that different BFNs had variable properties of the dynamic low-frequency fluctuations; also, the non-zero standard deviation (std) FALFF value for each BFN showed variant properties of the dynamic low-frequency fluctuations across different sliding windows.Fig. 1B and 1C demonstrated the power spectra of the estimated BFNs with lowest (Comp7) and largest (Comp5) mean FALFF values, where the power spectra were estimated from the first sliding window and the last sliding window, respectively.Based on the results in Fig. 1A-C, we could draw a conclusion that the simulated BFNs had different properties of the dynamic low-frequency fluctuations; also, the block-task modulated BFNs, e.g., Comp 5 and Comp 3, had more consistently dynamic low-frequency fluctuations across all the sliding windows than those without task modulation.
Fig. 1D depicted the different mean FALFF values for the estimated BFNs by Quasi-GICA across the sliding windows on the visual task data of Subject 1. Fig. 1E and 1F displayed the power spectra of the estimated BFNs with lowest (Comp12) and largest (Comp5) mean FALFF values, where the power spectra were estimated from the first sliding window data and the last sliding window data, respectively.According to the results in Fig. 1D-F, it can be concluded that the estimated BFNs had different properties of the dynamic low-frequency fluctuations, and the visual network modulated by visual stimulus, i.e., Comp5, was more consistent on dynamic lowfrequency fluctuations across all the sliding windows than the others.Meanwhile, the results of the dynamic low-frequency fluctuations for Subject 2 and 3 also showed the similar phenomenon, which were not depicted for saving space.

Results of Dynamic Spatial Reproducibility
Fig. 2A showed the mean and std values of DSR as to BFNs across the sliding windows on the simulated data, and Fig. 2B and 2C orderly displayed the normalized spatial maps of BFNs with lowest and highest DSR values, which were estimated from the first, middle and last sliding window data, respectively.According to the results in Fig. 2A-C, twelve estimated BFNs from simulation data showed different properties of dynamic spatial reproducibility, while the task-modulated BFN, i.e., Comp5, had the least spatial variance as the time went by, and Comp8 had the largest spatial variance across the sliding windows.This finding implied that the designed stimulus could strongly modulate the consistent patterns of brain activity in spatial domain.
Fig. 2D depicted the mean and std values of DSR as to BFNs across the sliding windows on the visual task data of subject 1, and Fig. 2E and 2F orderly displayed the spatial maps of BFNs with lowest and largest DSR values, which were estimated from the first and last sliding window data, respectively.Similarly, the visual stimulus modulated BFN, i.e., Comp5, had the strongest spatial reproducibility across the sliding windows, while Comp10 had the largest spatial variance as the time went by.This phenomenon implied that each BFN had different degrees of the time-varying spatial variance, and the task-modulated BFNs had strong spatial reproducibility across the sliding windows in contrast to the ones without modulation.Meanwhile, the similar performance in dynamic spatial reproducibility for Subject 2 and 3 was also found, whose results were not displayed for saving space.

Results of Weighted BFNs Ranking
Fig. 3 depicted the sorted results of the estimated BFNs from simulation data, where two estimated BFNs sharing the same task stimulus (Comp5 and Comp3) were ranked in the top, with the largest WRC values, while Comp7 ranked at the bottom with least WRC value, mainly due to its DFALFF value.Fig. 4 showed the sorted results of the estimated BFNs from visual task data of Subject 1, where the estimated BFNs corresponding to the visual stimulus (Comp5) ranked at the top, with the largest WRC value, while Comp12 ranked at the bottom with least WRC value, mainly due to its lowest DFALFF value (see Fig. 2D).Based on the results from the simulated and real data, we could draw a conclusion that the proposed BFNs weighted ranking method could effectively sort the meaningful dynamical BFNs at the top.

Discussion
At beginning, taking the results of dynamic low-frequency fluctuations (Fig. 1) and dynamic spatial reproducibility (Fig. 2) of the estimated BFNs together, we could find that the BFNs were not always with consistently dynamic properties of spatial reproducibility and low-frequency fluctuations.For example, based on the results of the visual task data of Subject 1, Comp12 was with the lowest FALFF value, while Comp10 had the lowest DSR value.On one hand, this phenomenon demonstrated the complexity of dynamical BFNs in spatial and spectral domain.On the other hand, the good ranking results depicted in Figs. 3 and 4 showed the effectiveness of the weighted BFNs ranking model based on the dynamic properties of power spectrum and spatial variance, while most of ranking methods focused on the static BFNs ranking only, e.g., power spectrum ranking, SVM-based ranking, MMC ranking, etc.Moreover, in contrast to power spectrum ranking [13], the weighted BFNs ranking is not subject to task-related BFNs ranking, and is also useful for the resting-state data.Compared with SVM-based ranking [14], the weighted BFNs ranking is simpler, having no prerequisite training process.Compared to the multiple ICA-estimationbased ranking [15][16], the only one-time ICA estimation is needed in Quasi-GICA.With respect to MMC ranking [17], the weighted BFNs ranking takes advantage of dynamical properties of spatial variance and power spectrum from the sliding windows, while MMC ranking uses the spatiotemporal reproducibility from the odd and even sampled dataset and needs at least two independent runs of ICA on sub-sampled dataset, which also suffers from some deficiencies, i.e., inaccurate estimation of order number on sub-sampled data, etc.

Conclusion
In this study, based on the dynamical BFNs tracking model, i.e., Quasi-GICA, a novel weighted BFNs ranking method was proposed, founding on the dynamic properties of spatial reproducibility in spatial domain and low-frequency fluctuations in spectral domain, aiming at determining the effective sorting orders of dynamical BFNs.Our results depicted and verified that different BFNs indeed had distinct dynamics in spatial variance and power spectrum in time-varying process, and the sorted results also demonstrated the effectiveness of this proposed model, which could sort BFNs with strongly dynamic spatial reproducibility and dynamic low-frequency fluctuations at the top.

F
is the N-by-L reducing matrix (determined by PCA), and N is the size of the re-tained time dimension.Considering the data compression on each ˆi X ,

1 kDPR
in equation(10) emphasizes the DFALFF property of BFNs from the different sliding-windows, while 2 k DPR in equation (11) emphasizes FALFF property of BFNs from the whole sliding-windows.For sorting the kth BFN, the dynamic spatial reproducibility, i.e., first calculated, and then the weighted ranking coefficient (WRC) is expressed as contribution to the time-varying properties of DSR and DFALFF.It is worth noting that the larger WRC implies that the corresponding BFN has stronger properties of the dynamic spatial reproducibility and dynamic low-frequency fluctuations.

Fig. 1 .
Fig. 1. (A)-(C) Dynamic properties of the low-frequency fluctuations regarding BFNs estimated from the simulated data; (D)-(F) Dynamic properties of the low-frequency fluctuations regarding BFNs estimated from the visual task data of Subject 1.

Fig. 2 .
Fig. 2. (A)-(C) Dynamic properties of the spatial reproducibility regarding BFNs from the simulated data; (D)-(F) Dynamic properties of the spatial reproducibility regarding BFNs from the visual task data of Subject 1.

Fig. 3 .
Fig. 3.The sorted order of the BFNs estimated from the simulated data.

Fig. 4 .
Fig. 4. The sorted order of the BFNs estimated from the visual task data of Subject 1.