PEAR: PEriodic And fixed Rank separation for fast fMRI.

PURPOSE: In functional MRI (fMRI), faster acquisition via undersampling of data can improve the spatial-temporal resolution trade-off and increase statistical robustness through increased degrees-of-freedom. High-quality reconstruction of fMRI data from undersampled measurements requires proper modeling of the data. We present an fMRI reconstruction approach based on modeling the fMRI signal as a sum of periodic and fixed rank components, for improved reconstruction from undersampled measurements.
METHODS: The proposed approach decomposes the fMRI signal into a component which has a fixed rank and a component consisting of a sum of periodic signals which is sparse in the temporal Fourier domain. Data reconstruction is performed by solving a constrained problem that enforces a fixed, moderate rank on one of the components, and a limited number of temporal frequencies on the other. Our approach is coined PEAR - PEriodic And fixed Rank separation for fast fMRI.
RESULTS: Experimental results include purely synthetic simulation, a simulation with real timecourses and retrospective undersampling of a real fMRI dataset. Evaluation was performed both quantitatively and visually versus ground truth, comparing PEAR to two additional recent methods for fMRI reconstruction from undersampled measurements. Results demonstrate PEAR's improvement in estimating the timecourses and activation maps versus the methods compared against at acceleration ratios of R = 8,10.66 (for simulated data) and R = 6.66,10 (for real data).
CONCLUSIONS: This paper presents PEAR, an undersampled fMRI reconstruction approach based on decomposing the fMRI signal to periodic and fixed rank components. PEAR results in reconstruction with higher fidelity than when using a fixed-rank based model or a conventional Low-rank + Sparse algorithm. We have shown that splitting the functional information between the components leads to better modeling of fMRI, over state-of-the-art methods.

Introduction

Accelerating acquisition in functional MRI (fMRI) has gained significant attention in neuroimaging. Accelerated fMRI can provide data at higher frame rates or sampling bandwidths, leading to higher temporal degrees of freedom.1 This enables the use of more powerful and sophisticated analysis techniques.2 Alternatively, accelerated fMRI may be used to increase spatial resolution without sacrificing temporal fidelity, enabling time‐resolved studies of the functional organization of the brain at finer scales, like cortical layers3 or columns.4 In resting state fMRI, where the goal is to estimate brain connectivity networks of the subject, accelerated data acquisition can improve the estimation of resting state networks (RSNs).5

Numerous methods for accelerating acquisition of MRI data by exploiting its intrinsic structure and redundancy have been published. In clinical dynamic MRI [e.g., Cardiac MRI and Dynamic Contrast Enhanced (DCE) MRI], many methods are based on undersampling in the k‐t space.6, 7, 8 Since the introduction of compressed sensing (CS),9, 10, 11 accelerated CS‐based methods for clinical MRI have also been developed.12, 13, 14, 15, 16, 17 While some of these methods have been adapted to fMRI,18, 19, 20, 21, 22, 23, 24 the different nature of the fMRI signal (when compared to clinical dynamic MRI), e.g., its low variance of signal of interest and its limited spatial compressibility, limits the adoption of those implementations.

Recently, we introduced an approach for reconstruction of fMRI from undersampled measurements that is based on modeling fMRI as fixed‐rank, i.e., k‐t FASTER (FMRI Accelerated in Space‐time by means of Truncation of Effective Rank).25 It aligns with the common analysis approaches in fMRI, which estimate limited numbers of spatial and temporal components from the data. It has been demonstrated that k‐t FASTER provides higher quality of activation and resting state maps when compared with other methods for fMRI reconstruction from undersampled data. This method has been extended recently to include multi‐channel coil sensitivity information and more flexible radial‐Cartesian sampling,26 providing additional encoding information and more incoherent sampling of k‐space, resulting in robust recovery of task‐based fMRI data at acceleration factors of 10 times or higher.

Several other approaches for reconstruction of fMRI from undersampled measurements have been recently suggested.27, 28, 29, 30, 31 Aggarwal et al. examined enforcing a low‐rank model and signal sparsity.27 Others explored exploiting low rank and sparsity after the separation of fMRI into two components. This approach, known as Low‐rank plus Sparse (or L + S),32 consists of modeling the data as a sum of two components, where low rank is enforced on one of them, and sparsity in some transform domain is enforced on the other. L + S has been examined for both clinical dynamic MRI33 and fMRI.28, 29, 30, 34

The common implementation of L + S for clinical dynamic MRI and fMRI consists of modeling the low rank component as a “background” image while the sparse component contains the dynamic information. This approach leads to satisfactory results for clinical dynamic MRI applications where the dynamic signal is significantly above the noise level [e.g., MR angiography (MRA) and DCE‐MRI] or periodic in the time domain (e.g., cardiac MRI). However, based on our experiments, its performance for fMRI, where the signal is often near the noise level and filtered by the hemodynamic response, is sub‐optimal.

In this study, we examine a different separation of the data, where both components contain functional information. While most previous methods that combine low rank and sparsity are based on solving an unconstrained minimization problem by singular value soft‐thresholding (SVT),35 in our approach we solve a constrained minimization problem based on truncating the singular values (a.k.a Truncated SVD or TSVD).36 Our approach forces one of the components to have a moderate fixed rank, and the other to be sparse in the temporal Fourier domain, leading to improved results compared to an SVT‐based approach. Reconstruction is performed via alternating minimization, that enforces the fixed‐rank requirement and sparsity iteratively.

We call our approach PEAR: PEriodic And fixed Rank separation for fast fMRI. We examine reconstructions from undersampled data acquired using golden‐angle radial sampling,37 and correspondence to both time‐course information [using General Linear modeling (GLM)] and spatial information (resting state network maps estimated via dual regression38). Our experiments consist of a purely synthetic simulation, to show the concept of separation between the components, a synthetic simulation using real timecourses to examine correspondence of results to real and known time‐courses, and retrospective sampling of a real resting state fMRI dataset, to examine resting state network recovery. We compare PEAR to k‐t FASTER that uses a fixed‐rank model only, and to a conventional, SVT‐based L + S implementation. We also explore the contributions of the different components in PEAR. Based on our experiments, PEAR exhibits better estimation of the timecourses and resting state networks from undersampled data, compared to k‐t FASTER and to L + S, using only 9.38% of the data in the simulations and 10% of the data in the retrospective sampling experiment.

The paper is organized as follows. Section 2 presents the proposed method for faster fMRI via separation of signals into periodic and fixed rank components. Section 3 describes experimental results. Section 4 discusses theoretical aspects and implementation details of our method and Section 5 concludes by highlighting the key results.

Method

In MRI, data is acquired in the spatial Fourier domain (k‐space). In dynamic MRI applications, such as cardiac MRI, MRA and DCE MRI, as well as in fMRI, the k‐space of each temporal frame is acquired. By undersampling k‐space (i.e., taking only partial k‐space measurements for each temporal frame), one can obtain higher frame rate, or alternatively, cover a greater extent in k‐space, thereby increasing spatial resolution without decreasing temporal resolution.

In the problem of fMRI reconstruction from undersampled k‐space, our goal is to recover the time series of acquired images. For simplicity, the time series is represented as a space‐time matrix, where each column is a 3D temporal frame concatenated as a vector, N denotes the number of pixels in a single frame, and T denotes the number of frames in the time series. The measurement model, which takes into account that in most cases data is acquired using multiple coils is:where y is a vector of undersampled measurements and E is a general linear operator that maps a matrix to a vector. For acquisition with multiple receiver coils, E consists of multiplication by coil sensitivities followed by an undersampled Fourier transform. The vector z represents the measurement noise, modeled as complex Gaussian with zero mean.

Since y is generated via undersampling, proper reconstruction of X from y requires assumptions on X. Relying on framework of CS, many methods that are based on sparsity of X in some transform domain were examined for dynamic MRI in general and for fMRI in particular. In our recent work, we considered modeling X as a fixed rank matrix, which aligns with the theory that X is composed of a relatively small number of spatially coherent temporal processes. The fixed‐rank based approach for fMRI (i.e., k‐t FASTER) solves the following minimization problem25:where r is a fixed, moderate rank that ranges between 20 and 50 in our fMRI approach (but may be much lower in other MRI modalities). Unlike some other types of dynamic MRI that exhibit high variance of signals of interest, in fMRI valuable information may also be embedded in low variance components, slightly above the noise level. Consequently, method to retrieve information from higher dimensions is expected to provide better results for fMRI.

An approach that was applied initially for clinical dynamic MRI,33 and has been examined recently for fMRI,28, 29, 34 consists of modeling the dynamic sequence as a sum of two components. A low‐rank component that represents mainly the background (denoted as the L component), and a component that contains the valuable signal. The latter is modeled as sparse in some transform domain, and denoted as the S component. This approach, known as L + S,32 is based on solving the following unconstrained problem:where L and S denote the low‐rank and sparse components, ‖·‖* and ‖·‖1 are the nuclear norm and the ℓ1 norm, λ 1,2 are tuning parameters that control the weight given to each term in the optimization problem and the reconstructed space‐time matrix is X = L+S. The linear transformation Ψ is a sparsifying transformation applied on S, depending on the specific dynamic MRI application. For MRA, Ψ may be chosen as an identity transform, whereas for cardiac MRI, which consists of periodic temporal structure, Ψ may be a temporal Fourier transform (i.e., applying a Fourier transform row‐wise, independently on each of the rows of S).33 For fMRI, both types of transformations were examined: Otazo et al.34 considered the identity transform (although their decomposition is used for analysis rather than acceleration) and Singh et al.28 examined the temporal Fourier transform. We note that L + S solves an unconstrained problem that does not explicitly enforce a fixed rank, and the solution is often based on SVT. By viewing the results of current implementations of L + S for fMRI28, 34 we found that in practice the resulting L component tends to have a very low rank and typically contains only background information, while the important functional information is in the S component.

Some fMRI analysis models suggest that while neural signals have strong band limited components, they also exist across the frequency spectrum.39 Therefore, we consider including sparsity in the temporal spectrum, to capture the bandlimited assumption, in addition to a fixed rank representation that would be more suitable for broader‐band signals. In particular, we propose modeling the fMRI signal as a sum of a fixed rank component, which contains the high variance information, and a periodic component that captures the periodicity that is not captured in the fixed rank component. Thus, we model the fMRI data X as X = A+P, where A and P are the fixed rank and periodic components, respectively. We enforce a limited number of temporal periodic signals for P, by demanding sparsity in the temporal Fourier domain, and a fixed rank for A which contains the high variance signal.

To understand the rationale behind this modeling, Fig. 1 shows 4 timecourses of arbitrarily selected pixels in a resting state fMRI dataset, that exhibited high correspondence (|Z| > 6) with regressors that represent a Default Mode Network (DMN) map or a visual network map after a dual‐regression against those network maps. The signals in the figure were extracted from the fMRI sequence after conventional fMRI pre‐processing (including skull stripping, motion correction, and slice timing correction), and the spatial locations of the selected pixels are shown at the bottom left of the figure. The timecourses are shown in the time domain (top left) and in the temporal Fourier domain, after removing the DC component (top right). It can be seen that while three of the timecourses exhibit peaks in the spectral domain and are suitable for sparse representation in the temporal Fourier domain, one of the timecourses exhibits a broad spectrum in the temporal Fourier domain. As a result, separation into fixed rank and periodic components allows better representation of signals (compared to fixed‐rank only or periodic component only) as it is expected to capture both broad‐band and band‐limited temporal spectra.

Figure 1

Top: timecourses (left) and amplitude spectrum (right) after removing DC of four selected pixels from a preprocessed resting‐state fMRI dataset. Bottom left: spatial locations of selected pixels. Pixels #1‐#3 showed high correspondence (|Z| > 6) with Default Mode Network and pixel #4 showed high correspondence with a visual network. The amplitude spectrum of the pixels shows that while pixels #1, #2 and #4 exhibit a limited number of peaks in the spectral domain (and therefore may be suitable for sparse modeling in the temporal Fourier domain) pixel #3 involves a wide range of spectral components. As a result, separation into a fixed‐rank and periodic components as proposed by our method would allow better representation of those signals, compared to using a fixed rank component (k‐t FASTER) or a periodic component only (L + S).

To obtain the separation into A and P components, we propose the following minimization problem:where is the set of matrices with a fixed rank r (ranges between 20 and 50 in our fMRI approach) and F is the temporal Fourier transform that applies a Fourier transform row‐wise on each of the rows of P (where ).

We solve Eq. (4) using alternating minimization (AM).40 In this approach, in each iteration we perform minimization with respect to one variable while keeping the other one fixed, and then switch between the variables. In our case, we start with an arbitrary initial point P 0. For n≥1 we iteratively compute: We solve each sub‐problem Eqs.(5), (6) via gradient projection41 where the proximal gradient is used for the non‐differentiable ℓ1 function in Eq. (6). A solution for Eq. (5) has been proposed by Goldfarb et al.42, a.k.a Iterative Hard Thresholding with Matrix Shrinkage (IHT‐MS). It consists of a gradient step for data consistency followed by a projection step onto the subspace . The general step is:where R (Q) is the projection onto the subspace , defined as: where σ 1 ≥ σ 2 ≥ … ≥ σ are the singular values of Q, and u and v are the singular vectors associated with σ . The operator S (Q) is the singular value soft‐thresholding operator, defined as: S (Q) = U[Σ−μ I]+ V where Σ = diag(σ 1,σ 2,…,σ ) and U and V are the left and right singular vectors associated with Σ. The parameter α is a step size, which controls the rate of convergence.

The rationale behind applying S (·) before applying R (·) can be explained by modeling the lower singular values, as representing noise, while represent data contaminated with additive noise. The subtraction of μ from is explained as removing noise from the higher singular values. Based on our experiments,25 the selection of μ = c·σ where c is a fixed parameter, leads to good results.

To solve Eq. (6) we use the Iterative Shrinkage‐Thresholding Algorithm (ISTA),43, 44 whose details are given in Appendix A. The general step for the solution of Eq. (6) is:where Λ indicates the soft‐thresholding operator with parameter λ, applied element‐wise. Finally, by defining X  = A  + P −α E {E{A  + P }−y}} we get that Eqs. (7) and (8) become: The convergence of an iterative solution based on Eq. (7) has been proven by Goldfarb et al. and the convergence of an iterative solution based on Eq. (8) is well studied in the literature.43, 45 Therefore, based on Csiszar and Tusndy,46 the convergence of our proposed AM framework is guaranteed.

The proposed algorithm is coined PEriodic And fixed Rank separation for fast fMRI (PEAR) and is summarized in Algorithm 1, where SVD represents the singular value decomposition and are the singular values of the matrix X −P in descending order. The operator is the conjugate temporal Fourier transform, and E is the conjugate transpose of E.

To summarize, the major differences between L + S and the PEAR approach for fMRI are outlined below:

Thresholding mechanism and a fixed rank solution: The solution of the L + S problem given in Eq. (2) using the same AM approach used for solving Eq. (4) results in an algorithm that is different from Algorithm 1 in the thresholding mechanism of the singular values. While in Algorithm 1 we perform singular value soft‐tresholding (SVT) with the value c·σ followed by truncating the r+1…m singular values (TSVD), for the solution of Eq. (2) we perform SVT with value λ 1, with no rank constraint. The issue of an SVT‐based solution versus a TSVD‐based solution has been studied in the literature previously,47 and it has been shown that TSVD performs better for fixed rank problems.47 Based on our experience, fMRI can be considered as a fixed rank problem, since the number activation networks is relatively small and in many cases kept fixed for analysis. This statement also aligns with our experimental results for fMRI hereinafter.

Separation of functional information between components: In conventional implementation of L + S for fMRI,28, 29, 34 the L component is modeled as very low rank, and therefore contains background information and no functional information. In PEAR, functional information is split between the components. Consequently, recovered functional information is not limited to periodic signals only, and results are improved compared to L + S, as will be shown in the next section. Indeed, the solution of Eq. (4) can be approximated by solving Eq. (2) with appropriate selection of λ 1,2 values. However, in PEAR we do not aim at separating between background information (which is typically modeled as a very low rank of 1 or 2) and functional information. Instead, we enforce a fixed rank (based on a priori knowledge of typical fMRI dimensionalities, often between 20−50), which allows functional information to be captured in both A and P components, and taking advantage of both fixed‐rank and sparsity modeling. The use of fixed rank simplifies the implementation of PEAR for a variety of datasets, as it obviates the need to examine a range of λ 1,2 values for each separate dataset. In addition, our experiments show that solving (4) provides better results for fMRI, when compared to solving Eq. (2) also for the case where λ 1,2 were carefully chosen for optimality.

To demonstrate the performance of PEAR compared to a well defined ground truth and additional algorithms for the reconstruction of fMRI from undersmpled measurements, we performed 3 types of experiments. In all experiments PEAR is compared to k‐t FASTER25 and L + S28 [where sparsity is enforced in the temporal Fourier domain, as shown in Eq. (2)].

The first experiment is a simulation based on synthetically generated mixtures of periodic and aperiodic signals, and aims to compare the results of PEAR to the aforementioned methods and to examine how PEAR separates the signals into periodic and fixed rank components. The second experiment is an extension of the first experiment using realistic time‐courses instead of purely synthetic ones. In the third experiment we examine the performance of PEAR for a retrospectively undersampled realistic 3D fMRI data sequence.

In all experiments, data is undersampled retrospectively, through a radial sampling approach using the golden‐angle37, 48 and the NUFFT49 package was used for forward and adjoint spatial Fourier transforms. The output of each experiment is provided as z‐statistics maps that reflect the degree to which each timecourse (in the case of experiments 1 and 2) or spatial regressor (in case of experiment 3) is expressed with a unique time‐course in the data. Output maps were null‐corrected using a Gaussian and Gamma mixture model.5

The parameters c = 0.7, α = 1 (for k‐t FASTER) and α = 0.5 (for PEAR and L + S) were selected experimentally. In all cases, all time‐points were initialized to the mean image calculated from all projections. All algorithms were run 100 iterations or until the minimum update between consecutive iterations was below 10−4. For k‐t FASTER and PEAR we examined the parameter r in the range between 1 and 50, and experimentally used r = 32 for k‐t FASTER, r = 27 for PEAR in experiments 1 and 2, and r = 20 for PEAR in experiment 3. In experiments 1 and 2 the parameters were tuned for optimal performance for each method (where optimal performance is evaluated by examining the z‐stat maps), and in experiment 3 parameters were tuned for optimality on a training set and results were obtained using the same parameters for an unseen fMRI sequence.

Experimental results

Experiment 1: purely synthetic simulation

In this experiment, we simulated a phantom consisting of 5 Regions of Interest (ROIs). Each ROI is spatially formed as a single letter from the letters “FMRIB”, and contains one of 5 purely synthetic timecourses, generated as follows. The letters “F” and “I” were purely periodic timecourses where “F” contains a single frequency and “I” contains a mixture of three frequencies, “R” was a purely aperdioc timecourse, and “M” and “B” were a superposition of periodic and aperiodic timecourses. The phantom was added to a realistic background fMRI dataset, to form a 2D fMRI sequence with known functional timecourses, of size 64 × 64, with 512 time points. The timecourses and their spatial locations in the simulated image are shown in Fig. 2 (left and top).

Figure 2

Left: spatial locations of the ROIs used in experiments 1 and 2, each ROI is formed as single letter and contains a single timecourse. Top: the timecourses used for each ROI in experiment 1. “F” and “I” are purely periodic timecourses where “F” contains a single frequency and “I” contains a mixture of three frequencies, “R” is a purely aperdioc timecourse, and “M” and “B” are superposition of periodic and aperiodic timecourses. Bottom: Timecourses used in experiment 2 for each simulated ROI (letter). [Color figure can be viewed at wileyonlinelibrary.com]

Undersampling was carried out in the k‐t space. We examined two undersampling ratios, first by taking 8 radial projections at each timepoint (corresponding to acceleration ratio of R = 8 relative to a fully‐sampled, maximally efficient Cartesian acquisition). We then repeated the experiment using only 6 radial projections at each timepoint (corresponding to R = 10.66). This simulates one slice of a hybrid radial‐Cartesian trajectory, which rotates an EPI trajectory within 3D k‐space.26 An additive white Gaussian noise with zero mean was added to the samples in the k‐space domain to obtain SNR of 25dB. For PEAR, λ was examined in the range of 0.45–3.4 and was selected as λ = 0.91 experimentally. For L + S, λ 1 was examined in the range of 1.1–3.4 and λ 2 was examined in the range of 0.45–3.4. These parameters were selected as λ 1 = 1.6 and λ 2 = 0.91 experimentally (the values for λ, λ 1,2 are provided after normalization with respect to the standard deviation of the data).

To examine the correspondence of the reconstruction with the ground‐truth, we performed regression against the original timecourses using General Linear Model (GLM).50 Figure 3 shows the F‐test results as null‐corrected z‐statistics maps,5 for the ground‐truth data (fully sampled image without the addition of noise), L + S, k‐t FASTER and PEAR, for both R = 8 and R = 10.66. All maps are thresholded at |Z| > 4.4 and shown with color scale mapped between 4.4 < |Z| < 15.

Figure 3

Experiment 1: GLM F‐test results of L + S, k‐t FASTER and PEAR for purely synthetic simulation, for R = 8 (top) and R = 10.66 (bottom). The z‐stat map of the ground truth is also shown (left). All maps are thresholded at |Z| > 4.4 and with color scale mapped between 4.4 < |Z| < 15. It can be seen that L + S is unable to recover the letter “R”, due to its purely aperiodicity. While PEAR exhibits some improvement compared to the other methods for R = 8, its superiority is clearly demonstrated for R = 10.66, as PEAR provides almost perfect recovery of the letters “R” and “B”, at minimal ratio of false positive errors. In addition, PEAR provides better statistical significance (compared to k‐t FASTER) for the recovery of “F”, “M”, “I” and “B”. [Color figure can be viewed at wileyonlinelibrary.com]

It can be seen that for R = 8, PEAR and k‐t FASTER provide similar results that outperform L + S, where PEAR provides some improvement compared k‐t FASTER. However, for R = 10.66 the superiority of PEAR is clearly demonstrated, as PEAR provides the most reliable result, being the only method that almost perfectly recovers both “R” and “B” with minimum false positive errors at R = 10.66. In addition, we see that L + S is unable to recover the aperiodic timecourse “R”, as opposed to both k‐t FASTER and PEAR thanks to their fixed‐rank component. Note that when comparing the results of PEAR to k‐t FASTER, the greater z‐statistics produced by PEAR should also be taken into account. We see that while both k‐t FASTER and PEAR show similar recovery patterns for “F”, “M”, and “I”, PEAR provides higher z‐stat values, demonstrating improved signal recovery and significance of activation detection.

An interesting analysis is the contribution of each component in PEAR. Figure 4 shows the GLM results for the A and P components of PEAR separately (for R = 10.66), where the z‐statistics maps are thresholded at |Z| > 4.9 and shown with color scale mapped between 4.9 < |Z| < 15. Note that the sum of the null‐corrected z‐statistics maps of A and P is not equal to the z‐statistics map of PEAR, due to the null‐correction applied for each map, that depends with each map's noise level. However, the separation of PEAR into periodic and fixed rank components is clearly demonstrated. The A component highly corresponds to the letter “R” which is a purely aperiodic timecourse, and with the letters “M” and “B” that include an aperiodic part. The P component highly corresponds to the letters “F” and “I” which are purely periodic timecourses, and to the letters “M” and “B” that include an periodic part. As demonstrated, this separation allows better modeling and leads to better recovery compared to k‐t FASTER and L + S. Another analysis is presented in Fig. 5, where example portions of the mean timecourses from the five letter ROIs are shown for the ground truth, L + S, k‐t FASTER and PEAR reconstruction results, including the timecourses for the A and P components of PEAR separately. The timecourses are shown in arbitrary units, to allow proper examination of their structure. It can be seen that as expected, L + S is limited in its ability to track the rapid changes that appear in the letter “R”. In addition, the P component of PEAR indeed contains the periodic part of the signal, and therefore exhibits high correspondence with letters that are fully periodic (“F” and “I”).

Figure 4

Experiment 1: GLM F‐test results of A and P components of PEAR for R = 10.66. Maps are thresholded at |Z| > 4.9 and with color scale mapped between 4.9 < |Z| < 15. It can be seen that the A component captures the letter “R” which is purely aperiodic, and the letters “M” and “B” that include an aperiodic part. The P component captures the letters “F” and “I” which represent periodic timecourses, and the letters “M” and “B” that include an periodic part. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 5

Example portions of the mean timecourses of the five letter ROIs are shown for the ground truth, L + S, k‐t FASTER and PEAR reconstruction results, including the timecourses for the A and P components of PEAR separately (Experiment 1, R = 10.66). The timecourses are shown in arbitrary units, to allow proper examination of their structure. It can be seen that L + S is limited in its ability to track the rapid changes that appear in the letter “R”. In addition, the separation of PEAR into A and P component is clearly seen, as the P component exhibits high correspondence with letters that are fully periodic (“F” and “I”), and A exhibits high correspondence with the aperiodic letter, “R”. [Color figure can be viewed at wileyonlinelibrary.com]

These simulations clearly demonstrate the expected behavior of our proposed approach under conditions where signals include pure periodicity; however, these are not a realistic depiction of fMRI data, even in task conditions, since these signals are rarely strongly periodic. The following experiment examines the performance of the various algorithms for realistic timecourses.

Experiment 2: simulation with realistic timecourses

In this experiment, we generated 5 realistic timecourses. First, we used a regression of a realistic fMRI dataset taken from the Human Connectome Project (HCP) database51 against 15 canonical resting state network maps (RSNs),38 derived from high‐dimensional group‐level Independent Component Analysis (ICA) of resting state fMRI datasets. The regression result provided 15 timecourses, each one corresponds to a single RSN regressor. We pulled 5 timecourses and used them instead of the purely synthetic timecourses used in experiment 1. These timecourses were used for the same simulation of the letters FMRIB, rather than retained in the original network spatial maps, due to the ease of visually evaluating the output parameter maps.

We repeated the same setting of experiment 1 (including the same SNR, sampling ratios, and selected parameters for each algorithm) where the only difference is the use of realistic timecourses instead of simulated ones. The realistic timecourses, including the corresponding regressors used for generation of the timecourses, are shown in Fig. 2 (bottom and left).

Figure 6 shows the General Linear Model (GLM) F‐test50 results as null‐corrected z‐statistics maps that were computed against the realistic time courses of all letters, for the ground‐truth data (fully sampled image without the addition of noise), L + S, k‐t FASTER and PEAR, for R = 8 and R = 10.66. All maps are thresholded at |Z| > 4.7 and shown with color scale mapped between 4.7 < |Z| < 15. In correspondence with experiment 1, for realistic timecourses we see that for R = 8, both k‐t FASTER and PEAR provide similar results that outperform L + S. For R = 10.66, we see that PEAR provides cleaner results, as can be seen mainly when comparing the recovery of the letters “R”, “I”, and “B”. Similarly to experiment 1, we need to examine the z‐statistics when comparing the results of PEAR to k‐t FASTER. It can be seen that although both k‐t FASTER and PEAR show similar recovery pattern for “F”, and “M”, PEAR results provide higher z‐stat values, demonstrating improved signal recovery.

Figure 6

Experiment 2: GLM F‐test results of ground truth, L + S, k‐t FASTER (for ranks of 27 and 32) and PEAR (for rank 27) for simulation with realistic timecoureses, for R = 8 (top)and R = 10.66 (bottom). The z‐stat map of the ground truth is also shown (left). All maps are thresholded at |Z| > 4.7 and with color scale mapped between 4.7 < |Z| < 15. It can be seen that while PEAR and k‐t FASTER provide similar results for R = 8, PEAR provides better results for R = 10.66 with minimal ratio of false positive errors. In addition, it can be seen that choosing the same rank (27) for kt‐FASTER and PEAR leads to different results, a fact that indicates the contribution of the P component on top of the fixed‐rank components. [Color figure can be viewed at wileyonlinelibrary.com]

In Fig. 6 we show the results of k‐t FASTER for two different ranks (r = 27,32). It can be seen that choosing the same rank (27) for kt‐FASTER and PEAR leads to different results, a fact that indicates the contribution of the P component on top of the fixed‐rank components.

To examine the ability of PEAR to recover nonstationarities in resting‐state patterns (we relate to the functional definition of “non‐stationary”, in which we consider that the connectivity between regions is non‐stationary, i.e., correlations between regions varying over time), we generated a new timecourse for the letter “B” by concatenating the first 256 timepoints of the “F” timecourse and the last 256 timepoints of the “R” timecourse. Together with the original ROI's and timecourses of “F” and “R” we repeated the experiment for R = 8 and R = 10.66. The timecourses used in this experiment and the results of PEAR are seen in Fig. 7. It can be seen that although “B” is non‐stationary, PEAR recovers “B” reliably.

Figure 7

PEAR results for non‐stationary timecourse: Left: timecourses used for “F”, “R” and “B”. “F” and “R” represent realistic timecourses of brain networks and “B” is composed of concatenating the first 256 timepoints of the “F” timecourse and the last 256 timepoints of the “R” timecourse. Right: GLM F‐test results of ground truth, and PEAR at R = 8 and R = 10.66. All maps are thresholded at |Z| > 5.6 and with color scale mapped between 5.6 < |Z| < 15. It can be seen that although “B” has a non‐stationary timecourse, PEAR provides reliable results in recovering it, for both R = 8 and R = 10.66. [Color figure can be viewed at wileyonlinelibrary.com]

Experiment 3: retrospective undersampling of real fMRI dataset

In this experiment, we examined kt‐FASTER, L + S and PEAR for the scenario of undersampled 3D fMRI data. We used simultaneous multi‐slice EPI data acquired on a 3T MRI system, with parameters based on HCP protocols: TR = 836 ms, multi‐band factor(MB) = 8, and isotropic resolution of 2 mm. Data was registered to an MNI standard space with dimensions 91×109×91 and included 512 timepoints. We examined two sampling ratios. First, we used 15% (R = 6.66) of the data by taking only 14 radial blades at each timepoint for each axial slice. In addition, we examined the scenario of using only 10% of the data (R = 10). To keep memory requirements under control, reconstructions were performed independently for each 91 × 109 2D axial slice with 512 timepoints. An additive white Gaussian noise with zero mean was added to the samples in the k‐space domain to obtain SNR of 25 dB. After reconstruction, all 2D reconstructed slices were stacked together to form a 3D image with 512 time point for further analysis. We applied a brain mask to focus only on within‐brain voxels.

After data reconstruction, we evaluated correspondence of the data to 15 canonical Resting State Networks (RSN) derived from high‐dimensional group‐level ICA of resting fMRI datasets from the HCP database.51 Evaluation was done using dual regression52 as follows. First, we performed spatial regression of the reconstructed dataset against the canonical maps (regressors) to extract the timecourses corresponding to each map. Then, we performed temporal regression of the dataset against the time series. The output is a set of z‐statistic maps (one for each regressor) that reflect the degree to which each spatial regressor is expressed with a unique time‐course in the data. More details about the dual regression process appear in Appendix B.

We compared the z‐statistics maps from PEAR to those computed from the fully sampled data (ground truth) and from reconstructions using L + S and k‐t FASTER methods. In this experiment, the parameters for each algorithm were tuned for a training sequence, and the results are evaluated using the same parameters for an unseen data sequence. For PEAR, λ was examined in the range between 0.25 and 2, and was selected as λ = 1.75 experimentally. For L + S, λ 1 was examined in the range between 0.25 and 1.75, and λ 2 was examined in the range between 0.13 and 1.75. These parameters were selected as λ 1 = 1.25 and λ 2 = 0.25 experimentally (the values for λ,λ 1,2 are provided after normalization with respect to the standard deviation of the data).

Figure 8 shows the Default Mode Network (DMN) regressor used for dual regression overlaid on the MNI atlas, as well as the z‐stat dual regression outputs of the ground truth, kt‐FASTER, L + S and PEAR for R = 6.66. All maps are thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8. The green ellipses show regions where differences between PEAR's activation patterns and L + S's and/or k‐t FASTER's activation patterns can clearly be seen. We note that PEAR provides better correspondence to ground truth in most activation regions compared to k‐t FASTER, and less false positive errors compared to L + S. While L + S shows the best activation pattern in the frontal lobe, evaluation of the results needs to take into account both correspondence of activation regions and false positive errors ratio, which are higher in L + S results. The trade‐off between the true and false positive errors is provided by the ROC curves, later in this paper.

Figure 8

Experiment 3: retrospective sampling of realistic fMRI dataset. Left: the regressor used for dual regression overlaid on the MNI template, and the dual regression results of the ground truth. Right: dual regression results of L + S, k‐t FASTER and PEAR obtained at undersampling ratio of 15% (R = 6.66). All maps are thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8. The ellipses illustrate regions where differences between PEAR to k‐t FASTER and/or L + S are clearly demonstrated. [Color figure can be viewed at wileyonlinelibrary.com]

The results for R = 10, are shown in Fig. 9. It can be seen that while L + S provides noisy results and k‐t FASTER underestimates the activation regions, PEAR maps are more similar to the ground truth. In addition, the phenomenon that has been demonstrated in experiments 1 and 2 repeats itself in real data – PEAR provides some improvement compared to the other methods for moderate acceleration ratio (R = 6), and significant improvement for high acceleration ratio (R = 10). Reconstruction results (R = 6.66)

Figure 9

Experiment 3: retrospective sampling of realistic fMRI dataset. Dual regression results of L + S, k‐t FASTER and PEAR obtained at undersampling ratio of 10% (R = 10) . All maps are thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8. It can be seen that while L + S provides noisy results and k‐t FASTER underestimates the activation regions (compared to the ground truth shown in Fig. 8), PEAR maps are more similar to the ground truth. [Color figure can be viewed at wileyonlinelibrary.com]

To provide measure for comparison between the methods, we computed the receiver operating characteristic (ROC) curves. This curve shows the performance of each method when compared to the ground truth (in terms of true positive ratio (TPR) against false positive ratio (FPR)) as the discrimination threshold varies. As a reference, we used the ground truth DMN z‐stat map shown in Fig. 8 (thresholded at |Z| > 3.3). For the generation of ROC we computed the TPR and FPR for each DMN z‐stat map for each method, as the threshold Z value ranges between 0 and 10.

A common scalar measure for the performance of the algorithm is the area under the curve (often referred to as AUC). Figure 10 shows the ROC curves for kt‐FASTER, L + S and PEAR, including the AUC computed for each curve, for R = 6.66 and R = 10. It can be seen that PEAR provides the most convex shape with the highest AUC in both cases. In addition, the degraded performance of L + S for R = 10 can clearly be seen by examining its ROC curve for R = 10. To check the validity of this result for different RSN maps, we computed the AUC for all the 15 canonical RSN maps used in our dual regression process. The results are shown in Fig. 11. We see that while there are cases in which k‐t FASTER or L + S outperform PEAR for some of the maps, PEAR provides the best AUC for the majority of the maps for both R = 6.66 and R = 10.

Figure 10

Receiver operating characteristic (ROC) curve for experiment 3. Performance of L + S, kt‐FASTER and PEAR are shown. The numbers in brackets indicate the area under the curve (AUC) for each method. The ROC results are generated for the ground truth map thresholded at |Z| > 3.3 serving as a reference. It can be seen that PEAR provides the highest AUC for DMN, for both R = 6.66 and R = 10. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 11

Area under ROC for 15 RSN maps, for R = 6.66 and R = 10 separately. It can be seen that PEAR provides the best performance for most of the maps, for both R = 6.66 and R = 10. [Color figure can be viewed at wileyonlinelibrary.com]

Finally, we examined the separation of PEAR into into A and P components, in terms of both time courses and spatial z‐stat maps (for R = 6.66). For this purpose, we first performed dual regression analysis for the A component and for the P component separately, to generate a z‐stat map for each. Those maps, thresholded at |Z| > 3.3 and with color scale mapped between 3.3 < |Z| < 8, are shown in Fig. 12 and demonstrate that both A and P components contain functional activity. Then, we arbitrary selected a single pixel that exhibited high correspondence with DMN for both A and P (z‐stat value of |Z| > 4). Figure 13 shows the timecourses and amplitude spectra of the A and P components, for the selected pixel, where the spatial location of the pixel is shown at the bottom of the Figure. The timecourses and amplitude spectra show that the A component contains a wide range of frequencies. The P component contains a limited number of temporal frequencies and captures the low‐energy periodicity that is not captured in A. This separation shows that a selection of a fixed, moderate rank for the A component, in addition to a demand for a limited number of temporal frequencies for the P component leads to the desired separation, which results in improved z‐stat results shown earlier.

Figure 12

Z‐stat maps of A and P components of PEAR for R = 6.66. Maps are thresholded at |Z| > 3.3 and with color scale mapped between 3.3 <| Z| < 8. It can clearly be seen that both components contain functional information and present activation regions in locations that correspond to similar activation regions in the ground truth. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 13

Top: Timecourses and amplitude spectra of a pixel resulted with |Z| > 4 for DMN z‐stat map. Timecourses and amplitude spectra of ground truth (top row) and A and P components of PEAR result (for R = 6.66) (second row) are shown. Values are shown in arbitrary units for better view. The spatial location of the selected pixel on an axial slice is also shown (left). It can be seen that the A component contains a wide range of frequencies. The P component contains a limited number of temporal frequencies and captures the periodicity that is not captured in the fixed‐rank component, A. This separation shows that a selection of a fixed, moderate rank for the A component, in addition to a demand for a limited number of temporal frequencies for the P component leads to the desired separation, which results in improved z‐stat results shown earlier. [Color figure can be viewed at wileyonlinelibrary.com]

Discussion

Relation to previous works

As described in the Introduction, a few methods that have been published recently also focus on separation of fMRI into two components.28, 30, 31, 34 The main differences between our work and these prior works lie in the methodology and experiments. First, we enforce both components to contain functional information explicitly, by solving an unconstrained minimization problem that enforces the A component to have a fixed moderate rank, using a TSVD based solution. This is in contrast to L + S‐based methods that practically enforce all the functional information to be contained in a single component, required to be sparse in some transform domain. As a result, the L + S solution might be sub‐optimal in reconstruction of signals that are neither periodic, nor strong enough to be captured in the low‐rank component.

Second, the experimental part of this paper presents a thorough validation of the suggested approach (compared to other existing methods), based on realistic nonuniform undersamping, and examines the spatial activation of resting state networks with broad spectrum characteristics. This analysis, which also includes the examination of the functional information that is contained in each of the components, is expanded compared to previous papers that deal with separation of fMRI into components; in particular, some of them show task‐based MRI where the design is periodic or basic RSN analysis.

Error measures and reproducibility

It is important to note that in the case of fMRI, conventional measures between the reconstructed datasets (e.g., MSE or correlations) may be misleading, as lower MSE does not necessarily mean that low variance functional information is preserved. Therefore, evaluation of results in this work is based the z‐stat analysis of activation maps, which is the common tool used today for resting state fMRI analysis.

In addition, while prospective undersampling is possible (and has been carried out in our previous works25), its evaluation has to be performed against connectivity maps taken from the literature or group‐averaged RSN spatial maps, leading to uncertainty in ensuring that subject‐specific detail is retained. Therefore, we focus on retrospectve undersampling, which allows accurate comparison against a well defined ground truth.

Limitations

This work focuses on resting state fMRI, which is a branch in fMRI research that deals with mapping brain connectivity based on an fMRI experiment that does not involve stimulation. In contrast, task‐based fMRI involves known manipulations of brain activity (a.k.a task‐based fMRI), although some of the properties of task‐based fMRI data are very similar to resting fMRI (i.e., the low variance signal based on the BOLD effect). Since many other reconstruction techniques place strong assumptions on brain activity, e.g., that it is periodic, highly reproducible or smooth in time, PEAR does not place these assumptions and therefore is expected to provide reliable results also for task‐base fMRI. However, the analysis of task‐based fMRI with PEAR is reserved for future research.

In addition, the reconstruction time of PEAR for a single realistic 2D axial slice in the dimensions described in our experimental results (109 × 91 with 512 time points) is approximately 15 min using MATLAB running on a single machine with 3.2 GHz CPU. To allow recovery of a multi‐slice image, we used cluster‐based computing that processed all 91 slices in parallel and provided a 3D volume approximately at the same running time. While the computation time and the need for cluster‐based computing are drawbacks of all iterative methods used to reconstruct fMRI sequences examined in this paper (k‐t FASTER and L + S), it is performed offline, while the subject is no longer in the scanner and does not require expensive MRI resources. In addition, we are examining approches to accelerate the reconstruction process, mainly by improving the NUFFT, using methods with lower computational complexity53 and implementing algorithms in a real‐time programming environment.

Conclusions

This paper presents PEAR, an under‐sampled fMRI reconstruction approach based on separating the fMRI signal to periodic and fixed‐rank components. The higher accelaration ratio offered by PEAR results in reconstruction with higher fidelity than when using a fixed‐rank based model or a conventional L + S algorithm. We have shown that splitting the functional information between the A and P components, by solving a constrained problem that enforces a fixed, moderate rank for the A component, leads to better modeling for fMRI, due to the unique nature of the fMRI signal. Future work will focus on extending this work to task‐based fMRI using both retrospective and prospective sampling.

Algorithm 1 PEAR: PEriodic And fixed Rank separation for fast fMRI
Input
Multicoil undersampled k‐t data: y
Space‐time multicoil encoding operator: E
Temporal Fourier transform: F t
Predefined rank: r, Soft‐shrinkage parameter: c
Tuning constant: λ, Step size: α
Iteration limit: N
Output: Estimated fMRI time‐series: X^
Initialize:
P 0 = 0, X 0 = E H{y}
Iterations for n = 1..N
UΣV H = SVD(X n−1−P n−1)
Σ(j,j)=Σ(j,j)−c·σ_r+1,j<randΣ(j,j)>c·σ_r+10,otherwise.
A n = UΣV H
P_n=FtH{Λ_λ(F_t{X_n−1−A_n−1})}
X n = A n+P n−α E H{E{A n+P n}−y}}