PEAR: PEriodic And fixed Rank separation for fast fMRI

TL;DR

PEAR introduces a novel fMRI reconstruction method that models signals as a combination of periodic and fixed rank components, enabling faster undersampled data recovery with improved accuracy.

Contribution

The paper presents PEAR, a new approach that separates fMRI signals into periodic and fixed rank parts, enhancing reconstruction quality from undersampled measurements.

Findings

PEAR outperforms existing methods at high acceleration ratios.

Reconstruction fidelity is higher with PEAR than traditional low-rank or sparse models.

PEAR effectively estimates timecourses and activation maps from undersampled data.

Abstract

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. We decompose the fMRI signal into a component which a has 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. 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,16 (for simulated data) and R=6.66,10 (for real data). 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.

Tables

Table 1. Table 1: Area under ROC for 15 RSN maps. The bold values indicate the method with the highest value in each line, 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. Training dataset results for R=6.66 are also shown for completeness.

RSN Map #	Training dataset			Testing dataset
	R=6.66			R=6.66			R=10
	L+S	kt-FASTER	PEAR	L+S	kt-FASTER	PEAR	L+S	kt-FASTER	PEAR
1	0.96277	0.96577	0.96781	0.96877	0.96883	0.97277	0.93883	0.95361	0.95615
2	0.97591	0.9769	0.97924	0.9736	0.97222	0.97431	0.94021	0.95687	0.95948
3	0.97874	0.97789	0.98244	0.97933	0.97794	0.98208	0.96465	0.96612	0.96712
4	0.9754	0.97446	0.97754	0.97658	0.97452	0.97728	0.92735	0.95559	0.95503
5	0.97529	0.97461	0.97708	0.97526	0.97688	0.97755	0.94546	0.95622	0.94724
6	0.97462	0.97433	0.976	0.97239	0.97238	0.97091	0.94525	0.94523	0.93694
7	0.96876	0.97043	0.97082	0.96444	0.96823	0.96846	0.92203	0.9363	0.93836
8	0.97306	0.97881	0.97935	0.9784	0.97938	0.98079	0.9628	0.96357	0.96555
9	0.95426	0.95392	0.95535	0.95299	0.95289	0.95464	0.91514	0.93129	0.92877
10	0.97721	0.97639	0.97778	0.97244	0.97562	0.97369	0.94051	0.94603	0.9479
11	0.96271	0.97008	0.96684	0.97813	0.97901	0.97835	0.95839	0.95413	0.95635
12	0.97276	0.97506	0.9747	0.96785	0.9698	0.9693	0.90773	0.94283	0.94328
13	0.97282	0.97967	0.98012	0.96625	0.97007	0.97086	0.92664	0.93289	0.93045
14	0.97	0.97148	0.9722	0.96661	0.96741	0.96622	0.92057	0.94641	0.93858
15	0.97133	0.97555	0.97695	0.97087	0.97282	0.97465	0.95221	0.94615	0.95443

Equations

$$\mathbf{y}=\mathbf{E}\{\mathbf{X}\}+\mathbf{z}$$

$$\underset{\mathbf{X}\in {\mathbb{R}}^{N\times T}}{\min }$$

$$\underset{\mathbf{L},\mathbf{S}\in {\mathbb{R}}^{N\times T}}{\min }$$

$$\underset{\begin{array}{c}\mathbf{A}\in \mathbb{C}\\ \mathbf{P}\in {\mathbb{R}}^{N\times T}\end{array}}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{E}\{\mathbf{A}+\mathbf{P}\}{\Vert }_2^2+\lambda \Vert {\mathbf{F}}_t\{\mathbf{P}\}{\Vert }_1$$

$${\mathbf{A}}_n=\underset{\mathbf{A}\in \mathbb{C}}{\text{arg min}}D(\mathbf{A},{\mathbf{P}}_{n-1})=\underset{\mathbf{A}\in \mathbb{C}}{\text{arg min}}\frac{1}{2}\Vert \mathbf{y}-\mathbf{E}\{\mathbf{A}+{\mathbf{P}}_{n-1}\}{\Vert }_2^2$$

$${\mathbf{P}}_n=\underset{\mathbf{P}\in R^{N\times K}}{\text{arg min}}D({\mathbf{A}}_n,\mathbf{P})=\underset{\mathbf{P}\in R^{N\times K}}{\text{arg min}}\frac{1}{2}\Vert \mathbf{y}-\mathbf{E}\{{\mathbf{A}}_n+\mathbf{P}\}{\Vert }_2^2+\lambda \Vert {\mathbf{F}}_t\{\mathbf{P}\}{\Vert }_1.$$

$${\mathbf{A}}_n=R_r(S_{\mu }({\mathbf{A}}_{n-1}-\alpha {\mathbf{E}}^H\{\mathbf{E}\{{\mathbf{A}}_{n-1}+{\mathbf{P}}_{n-1}\}-\mathbf{y}\}))$$

$${\mathbf{P}}_n={\mathbf{F}}_t^H\{{\Lambda }_{\lambda }({\mathbf{F}}_t\{{\mathbf{P}}_{n-1}-\alpha {\mathbf{E}}^H\{\mathbf{E}\{{\mathbf{A}}_{n-1}+{\mathbf{P}}_{n-1}\}-\mathbf{y}\}\})\}$$

$${\mathbf{A}}_n=R_r(S_{\mu }({\mathbf{X}}_{n-1}-{\mathbf{P}}_{n-1}))$$

$${\mathbf{P}}_n={\mathbf{F}}_t^H\{{\Lambda }_{\lambda }({\mathbf{F}}_t\{{\mathbf{X}}_{n-1}-{\mathbf{A}}_{n-1}\})\}.$$

$${\mathbf{P}}_n={\mathbf{F}}_t^H\{D({\mathbf{A}}_n,\mathbf{Q})\}$$

$$D({\mathbf{A}}_n,\mathbf{Q})=\underset{\mathbf{Q}\in R^{N\times K}}{\text{arg min}}\frac{1}{2}\Vert \mathbf{y}-\mathbf{E}\{{\mathbf{A}}_n+{\mathbf{F}}_t^H\{\mathbf{Q}\}\}{\Vert }_2^2+\lambda \Vert \mathbf{Q}{\Vert }_1.$$

$${\mathbf{Q}}_{k+1}={\Lambda }_{\lambda }({\mathbf{Q}}_k-\alpha {\mathbf{F}}_t\{{\mathbf{E}}^H\{\mathbf{E}\{{\mathbf{A}}_n+{\mathbf{F}}_t^H\{{\mathbf{Q}}_k\}\}-\mathbf{y}\}\})$$

$${\mathbf{P}}_n={\mathbf{F}}_t^H\{{\Lambda }_{\lambda }({\mathbf{F}}_t\{{\mathbf{P}}_{n-1}\}-\alpha {\mathbf{F}}_t\{{\mathbf{E}}^H\{\mathbf{E}\{{\mathbf{A}}_n+{\mathbf{P}}_{n-1}\}-\mathbf{y}\}\})={\mathbf{F}}_t^H\{{\Lambda }_{\lambda }({\mathbf{F}}_t\{{\mathbf{P}}_{n-1}-\alpha {\mathbf{E}}^H\{\mathbf{E}\{{\mathbf{A}}_{n-1}+{\mathbf{P}}_{n-1}\}-\mathbf{y}\}\})\},$$

Full text

Methods: The proposed approach decomposes the fMRI signal into a component which a has 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,16 (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.

fMRI, Low rank, Compressed Sensing

I 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 Feinberg et al. [2010]. This enables the use of more powerful and sophisticated analysis techniques Smith et al. [2012]. 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 layers Koopmans et al. [2011] or columns Yacoub et al. [2008]. 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) Beckmann et al. [2005].

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 Madore et al. [1999], Pruessmann et al. [1999], J et al. [2003]. Since the introduction of compressed sensing (CS) YC and G [2012], M et al. [2007], accelerated CS-based methods for clinical dynamic MRI have also been developed U et al. [2008], H et al. [2009]. While some of these methods have been adapted to fMRI Jung and Ye [2009], Zong et al. [2014], Jeromin et al. [2012], Holland et al. [2013], Chavarrías et al. [2015], Fang et al. [2015], Zhao et al. [2010], 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) Chiew et al. [2015]. 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 Chiew et al. [2016], 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 Aggarwal and Gupta [2016], Singh et al. [2015], AY et al. [2016], Aggarwal et al. [2017], Lam et al. [2013]. Aggarwal et al. examined enforcing a low-rank model and signal sparsity Aggarwal and Gupta [2016]. 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) Candès et al. [2011], 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 MRI Otazo et al. [2015] and fMRI Singh et al. [2015], R et al. [2015], AY et al. [2016], Aggarwal et al. [2017].

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) Cai et al. [2010], in our approach we solve a constrained minimization problem based on truncating the singular values (a.k.a Truncated SVD or TSVD) Gavish and Donoho [2014]. 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 samplingWinkelmann et al. [2007], and correspondence to both time-course information (using General Linear modeling (GLM)) and spatial information (resting state network maps estimated via dual regression F et al. [2009]). 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 6.25% of the data in the simulations and 10% of the data in the retrospective sampling experiment.

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

II 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, $\mathbf{X}\in\mathbb{R}^{N\times T}$ 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:

[TABLE]

where $\mathbf{y}$ is a vector of undersampled measurements and $\mathbf{E}$ is a general linear operator that maps a matrix to a vector. For acquisition with multiple receiver coils, $\mathbf{E}$ consists of multiplication by coil sensitivities followed by an undersampled Fourier transform. The vector $\mathbf{z}$ represents the measurement noise, modeled as complex Gaussian with zero mean.

Since $\mathbf{y}$ is generated via undersampling, proper reconstruction of $\mathbf{X}$ from $\mathbf{y}$ requires assumptions on $\mathbf{X}$. Relying on framework of CS, many methods that are based on sparsity of $\mathbf{X}$ in some transform domain were examined for dynamic MRI in general and for fMRI in particular. In our recent work, we considered modeling $\mathbf{X}$ as a fixed rank matrix, which aligns with the theory that $\mathbf{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 problem Chiew et al. [2015]:

[TABLE]

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 Otazo et al. [2015], and has been examined recently for fMRI AY et al. [2016], Singh et al. [2015], R et al. [2015], 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 Candès et al. [2011], is based on solving the following unconstrained problem:

[TABLE]

where $\mathbf{L}$ and $\mathbf{S}$ denote the low-rank and sparse components, $\|\cdot\|_{*}$ and $\|\cdot\|_{1}$ are the nuclear norm and the $\ell_{1}$ norm, $\lambda_{1,2}$ are tuning parameters that control the weight given to each term in the optimization problem and the reconstructed space-time matrix is $\mathbf{X}=\mathbf{L}+\mathbf{S}$. The linear transformation $\mathbf{\Psi}$ is a sparsifying transformation applied on $\mathbf{S}$, depending on the specific dynamic MRI application. For MRA, $\mathbf{\Psi}$ may be chosen as an identity transform, whereas for cardiac MRI, which consists of periodic temporal structure, $\mathbf{\Psi}$ may be a temporal Fourier transform (i.e., applying a Fourier transform row-wise, independently on each of the rows of $\mathbf{S}$) Otazo et al. [2015]. For fMRI, both types of transformations were examined: Otazo et el. R et al. [2015] considered the identity transform (although their decomposition is used for analysis rather than acceleration) and Singh et al. Singh et al. [2015] 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 fMRIR et al. [2015], Singh et al. [2015] 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 Thompson et al. [2014]. 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 $\mathbf{X}$ as $\mathbf{X}=\mathbf{A}+\mathbf{P}$, where $\mathbf{A}$ and $\mathbf{P}$ are the fixed rank and periodic components, respectively. We enforce a limited number of temporal periodic signals for $\mathbf{P}$, by demanding sparsity in the temporal Fourier domain, and a fixed rank for $\mathbf{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.

To obtain the separation into A and P components, we propose the following minimization problem:

[TABLE]

where $\mathbb{C}$ is the set of matrices with a fixed rank $r$ (ranges between 20 and 50 in our fMRI approach) and $\mathbf{F}_{t}$ is the temporal Fourier transform that applies a Fourier transform row-wise on each of the rows of $\mathbf{P}$ (where $\hat{\mathbf{X}}=\hat{\mathbf{A}}+\hat{\mathbf{P}}$).

We solve (4) using alternating minimization (AM) Niesen et al. [2009]. 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 $\mathbf{P}_{0}$. For $n\geq 1$ we iteratively compute:

[TABLE]

[TABLE]

We solve each sub-problem (5,6) via gradient projection Bertsekas [1999] where the proximal gradient is used for the non-differentiable $\ell_{1}$ function in (6). A solution for (5) has been proposed by Goldfarb et al. Goldfarb and Ma [2011], 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 $\mathbb{C}$. The general step is:

[TABLE]

where $R_{r}(\mathbf{Q})$ is the projection onto the subspace $\mathbb{C}$, defined as: $R_{r}(\mathbf{Q})=\sum_{i=1}^{r}\sigma_{i}\mathbf{u}_{i}\mathbf{v}_{i}^{H}$ where $\sigma_{1}\geq\sigma_{2}\geq...\geq\sigma_{m}$ are the singular values of $\mathbf{Q}$, and $\mathbf{u}_{i}$ and $\mathbf{v}_{i}$ are the singular vectors associated with $\sigma_{i}$. The operator $S_{\mu}(\mathbf{Q})$ is the singular value soft-thresholding operator, defined as: $S_{\mu}(\mathbf{Q})=\mathbf{U}[\mathbf{\Sigma}-\mu\mathbf{I}]_{+}\mathbf{V}^{H}$ where $\mathbf{\Sigma}=\text{diag}(\sigma_{1},\sigma_{2},...,\sigma_{m})$ and $\mathbf{U}$ and $\mathbf{V}$ are the left and right singular vectors associated with $\mathbf{\Sigma}$. The parameter $\alpha$ is a step size, which controls the rate of convergence.

The rationale behind applying $S_{\mu}(\cdot)$ before applying $R_{r}(\cdot)$ can be explained by modelling the lower singular values, $\{\sigma_{i}\}_{i=r+1}^{m}$ as representing noise, while $\{\sigma_{i}\}_{i=1}^{r}$ represent data contaminated with additive noise. The subtraction of $\mu$ from $\{\sigma_{i}\}_{i=1}^{r}$ is explained as removing noise from the higher singular values. Based on our experiments Chiew et al. [2015], the selection of $\mu=c\cdot\sigma_{r+1}$ where $c$ is a fixed parameter, leads to good results.

To solve (6) we use the Iterative Shrinkage-Thresholding Algorithm (ISTA) Daubechies et al. [2004], A and M [2009], whose details are given in Appendix A. The general step for the solution of (6) is:

[TABLE]

where $\Lambda_{\lambda}$ indicates the soft-thresholding operator with parameter $\lambda$, applied element-wise. Finally, by defining $\mathbf{X}_{n-1}=\mathbf{A}_{n-1}+\mathbf{P}_{n-1}-\alpha\mathbf{E}^{H}\{\mathbf{E}\{\mathbf{A}_{n-1}+\mathbf{P}_{n-1}\}-\mathbf{y}\}\}$ we get that (7) and (8) become:

[TABLE]

[TABLE]

The convergence of an iterative solution based on (7) has been proven by Goldfarb et al. and the convergence of an iterative solution based on (8) is well studied in the literature Daubechies et al. [2004], Combettes and Wajs [2005]. Therefore, based on Csiszar and Tusnády Csisz and Tusnády [1984], 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 $\{\sigma_{i}\}_{i=1}^{r+1}$ are the singular values of the matrix $\mathbf{X}_{n-1}-\mathbf{P}_{n-1}$ in descending order. The operator $\mathbf{F}^{H}_{t}$ is the conjugate temporal Fourier transform, and $\mathbf{E}^{H}$ is the conjugate transpose of $\mathbf{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 (3) using the same AM approach used for solving (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\cdot\sigma_{r+1}$ followed by truncating the $r+1...m$ singular values (TSVD), for the solution of (3) we perform SVT with value $\lambda_{1}$, with no rank constraint. The issue of an SVT-based solution versus a TSVD-based solution has been studied in the literature previously Josse and Sardy [2013], and it has been shown that TSVD performs better for fixed rank problems Josse and Sardy [2013]. 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 Singh et al. [2015], AY et al. [2016], R et al. [2015], 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 (4) can be approximated by solving (3) with appropriate selection of $\lambda_{1,2}$ values. However, PEAR enables enforcing a fixed rank (based on a priori knowledge of typical fMRI dimensionalities, often between 20-50) for a variety of datasets, obviating the need to examine a range of $\lambda_{1,2}$ values for each separate dataset. In addition, our experiments show that solving (4) provides better results for fMRI, when compared to solving (3) also for the case where $\lambda_{1,2}$ were carefully chosen for optimality.

III Experimental Results

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 FASTER Chiew et al. [2015] and L+S Singh et al. [2015] (where sparsity is enforced in the temporal Fourier domain, as shown in (3)).

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-angle Winkelmann et al. [2007], Graedel et al. [2016] and the NUFFT Fessler and Sutton [2003] 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 Beckmann et al. [2005].

The parameters $c=0.7$, $\alpha=1$ (for k-t FASTER) and $\alpha=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.

III.1 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\times 64$, with $512$ time points. The timecourses and their spatial locations in the simulated image are shown in Fig. 2 (left and top).

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 4 radial projections at each timepoint (corresponding to R=16). This simulates one slice of a hybrid radial-Cartesian trajectory, which rotates an EPI trajectory within 3D k-space Chiew et al. [2016]. An additive white Gaussian noise with zero mean was added to the samples in the k-space domain to obtain SNR of $25$dB. For PEAR, $\lambda$ was examined in the range of 0.45-3.4 and was selected as $\lambda=0.91$ experimentally. For L+S, $\lambda_{1}$ was examined in the range of 1.1-3.4 and $\lambda_{2}$ was examined in the range of 0.45-3.4. These parameters were selected as $\lambda_{1}=1.6$ and $\lambda_{2}=0.91$ experimentally (the values for $\lambda,\lambda_{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)Jenkinson et al. [2012]. Figure 3 shows the F-test results as null-corrected z-statistics mapsBeckmann et al. [2005], 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=16. All maps are thresholded at $|Z|>4.3$ and shown with color scale mapped between $4.3<|Z|<15$.

It can be seen that PEAR provides the most reliable result, being the only method that almost perfectly recovers both “M” and “B” with minimum false positive errors at R=16. 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.

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=16), where the z-statistics maps are thresholded at $|Z|>4.3$ and shown with color scale mapped between $4.3<|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 with 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 modelling 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”).

These simulations clearly demonstrate the expected behaviour 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.

III.2 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) database Van Essen et al. [2013] against 15 canonical resting state network maps (RSNs) F et al. [2009], 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-test Jenkinson et al. [2012] 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=16. All maps are thresholded at $|Z|>4$ and shown with color scale mapped between $4<|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=16, we see that PEAR provides cleaner results, as can be seen mainly when comparing the recovery of the letters “F” and “B”. It can also be seen that all methods are unable to recover “R” and “I” due to the high undersampling ratio. It can also be seen that “R” and “I” are the letters with the lowest Z values in the fully-sampled, ground truth data, due to their low energy in this experiment.

III.3 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, taken from the HCP database Van Essen et al. [2013]. Data was registered to an MNI standard space with dimensions $91\times 109\times 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 91x109 2D axial slice with 512 timepoints. After reconstruction, all 2D reconstructed slices were stacked together to form a 3D image with 512 time point for further analysis.

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 Van Essen et al. [2013]. Evaluation was done using dual regression Beckmann et al. [2009] 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.

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, $\lambda$ was examined in the range between 0.25 and 2, and was selected as $\lambda=1.75$ experimentally. For L+S, $\lambda_{1}$ was examined in the range between 0.25 and 1.75, and $\lambda_{2}$ was examined in the range between 0.13 and 1.75. These parameters were selected as $\lambda_{1}=1.25$ and $\lambda_{2}=0.25$ experimentally (the values for $\lambda,\lambda_{1,2}$ are provided after normalization with respect to the standard deviation of the data).

Figure 7 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 results for R=10, are shown in Fig. 8. The green ellipses in the images show regions where PEAR’s activation pattern is the most similar one to the ground truth. It can also be seen that while both k-t FASTER and PEAR provide reliable results for both R=6.66 and R=10, L+S is does not provide satisfactory results at the higher acceleration ratio of R=10.

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. 7 (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 9 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 summary of the results is shown in Table 1, where the method that provides the highest AUC for each RSN in marked in bold (for R=6.66 and R=10 separately). 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. The results of the training dataset for R=6.66 are also shown in the table, for completeness.

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. 10 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.5$). Figure 11 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, with some strong peaks in the spectrum for frequencies that have strong total energy. 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.

IV Discussion

IV.1 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 Lam et al. [2013], Singh et al. [2015], Aggarwal et al. [2017], R et al. [2015]. 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.

IV.2 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 works Chiew et al. [2015]), 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.

IV.3 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\times 91$ with 512 time points) is approximately 15 minutes using MATLAB running on a single machine with 3.2GHz 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 complexityKiperwas et al. [2017] and implementing algorithms in a real-time programming environment.

V 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.

Acknowledgements

This work was supported by the Israeli Ministry of Science, by the ISF I-CORE joint research center of the Technion and the Weizmann Institute, Israel, by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 646804-ERC-COG-BNYQ, by the British Technion Society and the Coleman-Cohen Fellowship for post-doctoral studies in the UK, and by the EPSRC and the Wellcome Trust.

The authors have no relevant conflicts of interest to disclose.

Appendix A Solution of (6) using ISTA

Let $\mathbf{Q}=\mathbf{F}_{t}\{\mathbf{P}\}$. Using the fact that $\mathbf{F}_{t}$ is unitary, (6) can be written as:

[TABLE]

where

[TABLE]

An iterative solution to (12) can be obtained using the iterative shrinkage-thresholding algorithm (ISTA)Daubechies et al. [2004], Combettes and Wajs [2005], A and M [2009] , whose general step is:

[TABLE]

Here $\Lambda_{\lambda}$ is the soft-thresholding operator with parameter $\lambda$, $\alpha$ is the step size and $\mathbf{Q}_{0}=\mathbf{P}_{n-1}$. Using the fact that $\mathbf{P}_{n}=\mathbf{F}_{t}^{H}\{\mathbf{Q}_{K}\}$, setting $K=1$ and using (11) we get that the general step for the solution of (6) is:

[TABLE]

where the last equality is (8) in the paper.