Temporal Non-Local Means Filtering Reveals Real-Time Whole-Brain Cortical Interactions in Resting fMRI

Intensity variations over time in resting BOLD fMRI exhibit spatial correlation patterns consistent with a set of large scale cortical networks. However, visualizations of this data on the brain surface, even after extensive preprocessing, are dominated by local intensity fluctuations that obscure larger scale behavior. Our novel adaptation of non-local means (NLM) filtering, which we refer to as temporal NLM or tNLM, reduces these local fluctuations without the spatial blurring that occurs when using standard linear filtering methods. We show examples of tNLM filtering that allow direct visualization of spatio-temporal behavior on the cortical surface. These results reveal patterns of activity consistent with known networks as well as more complex dynamic changes within and between these networks. This ability to directly visualize brain activity may facilitate new insights into spontaneous brain dynamics. Further, temporal NLM can also be used as a preprocessor for resting fMRI for exploration of dynamic brain networks. We demonstrate its utility through application to graph-based functional cortical parcellation. Simulations with known ground truth functional regions demonstrate that tNLM filtering prior to parcellation avoids the formation of false parcels that can arise when using linear filtering. Application to resting fMRI data from the Human Connectome Project shows significant improvement, in comparison to linear filtering, in quantitative agreement with functional regions identified independently using task-based experiments as well as in test-retest reliability.


B Laplace-Beltrami filtering
Gaussian filtering is frequently used to isotropically smooth volumetric 3D data in functional and task based fMRI studies. Isotropic smoothing can also be performed on a 2D manifold such as the cortical surface. Because cortical surfaces have an intrinsic curvature that prevents them being mapped without distortion onto the plane, isotropic smoothing requires that curvature be taken into account. The Laplace-Beltrami (LB) operator, which generalizes the 2D Gaussian smoothing kernel from the plane to an arbitrary smooth manifold [4], can produce isotropic smoothing on cortex.
In this approach, a solution of the heat equation ∂f (p, t) ∂t = ∆f (p, t) with initial condition f (p, 0) = Y (p) is used to smooth the data Y (p), where ∆ is the LB operator on the cortical manifold with spatial parameters p [5,6,7]. The parameter t represents the diffusion time for the heat equation. The degree of smoothing increases with t. The solution is given by [8, chap. 1] f (p, t) = Y (q) K t (p, q) dq (S. 1) where, K t (p, q) is the heat kernel which can be expressed in terms of ordered eigen-functions φ 0 , φ 1 , φ 2 , · · · and corresponding eigen-values λ 0 , λ 1 , λ 2 , · · · of the Laplace-Beltrami operator ∆ as [8,6,7] K t (p, q) = Eq. (S.3) can be expressed in discrete form as a square matrix of dimension equal to the number of vertices on the cortical surface using a truncated eigen-function expansion of the LB operator (we use the first 800 eigenfunctions in all results with LB filtering in this paper). The level of smoothing of the final solution is parameterized by the time parameter t. The weights applied to neighboring vertices while filtering a particular vertex for different values of t can be intuitively interpreted as the smoothing kernel at the vertex. This approach allows efficient computation of filtered images as the eigen-decomposition needs to be computed only once irrespective of the number of rfMRI volumes or level of smoothing.

C Mapping and measures between parcellations
Assume we have two parcellations A and B for a set of vertices V :

Label-wise agreement
Agreement is a label-wise measure of the fraction of vertices that agree between two parcellations A and B. Specifically agreement of parcel A i ∈ A with the corresponding parcel in B is given by where | · | represents the cardinality of a set and S : Z M → Z N is a mapping of parcels in A to B as described later. The agreement measure ranges from 0 (no agreement) to 1 (perfect agreement).

Concordance: consistency between parcellations
Concordance is a global measure of consistency between two parcellations A and B. It is defined as the fraction of vertices between the two parcellations that agree: where | · | represents the cardinality of a set and S : Z M → Z N is mapping of parcels in A to B as described next. The concordance measure ranges from 0 (no agreement) to 1 (perfect agreement).

Mapping between parcellations
Given the two parcellations A and B, we aim to match each parcel in A to a unique parcel in B. Let g(A i , B j ) be a measure of the goodness of the match of A i to B j such as the Dice coefficient or Jacard index. We want a map S : Z M → Z N such that we maximize the goodness of match across all parcels: The exact solution of eq. (S.3) is combinatorial and scales approximately as n! where n = max(M, N )). We use an approximate solution of eq. (S.3) by noting its similarity to the famous stable matching problem [9]. Stable matching finds a match between elements of two sets of equal size when a preference order of matching is specified for each element. A match is stable if there does not exists a pair (a, b) in the match in which both a and b have higher preference elements which also prefer a and b respectively over their current match. We use the Gale-Shapley algorithm [9] after transforming our parcel mapping to a stable matching problem as described below. We will match each parcel in A to a unique parcel in B: In the language of the Gale-Shapley algorithm, elements in A are suitors and the elements in B are reviewers.
2. The Gale-Shapley algorithm works with sets of the same size. Hence, we define an n × n matrix G by appending an appropriate number of row or columns to G. All the elements of appended rows/columns are set to δ = min(G) − , where is a small positive constant. This modifies our problem by adding dummy suitors/reviewers which can be easily ignored.
3. Next, we compute the preference order of each element from G. There are a total of n suitors and n reviewers. We define the preference order for each suitor by arranging the indices of the elements in the corresponding row in G in the descending order of magnitude of their entries. Similarly, we define the preference order of each reviewer as the indices of the elements of the columns arranged in descending order of magnitude of their entries. If two elements have same preference, then we break the tie by randomly assigning a preference order.
4. We use the resulting preference order from the last step to find a stable match S using the Gale-Shapley algorithm [9]. The map S : Z n → Z n is modified to getŜ : Z M → Z N as follows: If M = N thenŜ ≡ S; If M < N then all the appended dummy suitors are ignored in S; If M > N then all suitors which are matched to an appended dummy reviewer in S are modified to match to an empty set so that BŜ (i) = ∅.
The matching solution obtained by the above procedure is the suitor-optimal solution (in the sense of preference order). If a reviewer-optimal solution is required A and B should be swapped [9]. Further, the solution obtained is only an approximation of eq. (S.3) as the Gale-Shapley algorithm is blind to the absolute values of g A i , B S(i) and only uses relative preference order, which may not always maximize the total cost. Nonetheless, in our experience, this approximate solution produces reasonable matching and is more computationally tractable than combinatorial approaches.      Figure H: An example of task labels for a single subject obtained by thresholding the task activation statistical maps at Z-score of 3.0 followed by ignoring all connected components with less than 40-vertices (see section 2.4.3 for more details). A total of 17 task-pair activation maps were obtained for each subject from HCP's processed data, which were computed by processing the task fMRI data with three different level of smoothing: (a) 2mm, (b) 4mm, and (c) 8 mm. The task-pairs corresponding to each task label is reported in (d). The 4mm smoothing result is also shown in Fig. 4 in the body of the paper and used for subsequent quantitative analysis in Fig. 10 and Table 1.  Only a subset of these results are included in the body of the paper. Specifically, overall performance results for 4mm smoothing (Fig. 10(b)) and mean agreement fraction for left foot task and 4mm smoothing ( Fig. 10(a)).  Detailed performance of all filtering approaches across several parameters for task-labels obtained by 2mm smoothing of task fMRI data. Line plots show the mean agreement of N-cuts clusters with 2mm-smoothed task labels across 40 subjects with 4 sessions each (40 × 4 = 160 N-cuts clustering results). For each task, we also show the spatial histogram of the task label across the 40 subjects. Only seven tasks, shown above, could survive the statistical and contiguity thresholding (see section 2.4.3 for more details), and can also be seen in Fig. H. The tongue task was also studied after breaking the activation maps into two parts, each lying on one hemisphere of the brain.

Quantitative comparison with probabilistic Brodmann areas (BAs)
Since, it is believed that the cytoarchitectonic areas reflect functional specialization in cerebral cortex [10], it is also interesting to study how the functional parcellations obtained from rfMRI data compare with probabilistic Brodmann areas (BA). Probabilistic BA were obtained using histology studies of 10 postmortem human brains and were transferred to the subject's 32K Conte-69 cortical mesh, which are available from HCP as discrete labels (after thresholding) for each subject [11,12]. An example of the discrete Brodmann label map is shown in Fig. M(a). Similarly to our task based analysis, we computed label-wise agreement between BA labels and N-cuts parcellations to study the effect of different filtering approaches.
tNLM > LB 0.146 Table 1: Statistical tests for higher agreement with probabilistic Brodmann areas: Table of (uncorrected) p-values for one-sided tests for 'best' performance of LB and tNLM filtering (parameters of 'best' parcellation are reported in Fig. M(b)). The procedure and description of the non-parametric tests are same as that of task labels.
We also computed quantitative agreement of N-cuts parcellations with probabilistic BAs. Fig. M(b) shows the summary of the best performances for different filtering approaches for each BA. Consistent with our previous observation for task fMRI, tNLM achieves the highest agreement fraction across several BA (for 18 out of 24 BAs). We also performed the non-parametric statistical tests with the same null hypothesis as the task analysis. The results and alternate hypothesis of the statistical tests are reported in table 1 which reveal that tNLM filtering significantly improves the agreement of several BA labels with the parcellation results as compared with LB filtering. There are few BA labels, such as BA 3a (L), for which LB parcellations shows higher median agreement than tNLM. More complete data are included in Fig. N in which mean agreement is plotted for each method and parameter for each of the BAs.