Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2659356

Formats

Article sections

Authors

Related links

Neuroimage. Author manuscript; available in PMC 2009 August 1.

Published in final edited form as:

Published online 2008 April 23. doi: 10.1016/j.neuroimage.2008.04.179

PMCID: PMC2659356

NIHMSID: NIHMS63385

Fa-Hsuan Lin,^{1,}^{2} Thomas Witzel,^{3} Joseph B. Mandeville,^{1} Jonathan R. Polimeni,^{1} Thomas A. Zeffiro,^{4} Douglas N. Greve,^{1} Graham Wiggins,^{1} Lawrence L. Wald,^{1} and John W. Belliveau^{1}

Correspondence: Fa-Hsuan Lin, MGH-HST Athinoula A. Martinos Center Rm 2301, Bldg. 149 13^{th} St., Charlestown MA 02129, Email: ude.dravrah.hgm.rmn@nilhf, Tel: 617-726-8954, Fax: 617-726-7422

The publisher's final edited version of this article is available at Neuroimage

See other articles in PMC that cite the published article.

Developments in multi-channel radio-frequency (RF) coil array technology have enabled functional magnetic resonance imaging (fMRI) with higher degrees of spatial and temporal resolution. While modest improvement in temporal acceleration has been achieved by increasing the number of RF coils, in parallel data acquisition techniques, the maximum attainable acceleration is intrinsically limited only by the amount of independent spatial information in the combined array channels. Since the geometric configuration of a large-*n* MRI head coil array is similar to that used in EEG electrode or MEG SQUID sensor arrays, the source localization algorithms used in MEG or EEG source imaging can be extended to also process MRI coil array data, resulting in greatly improved temporal resolution by minimizing k-space traversal during signal acquisition. Using a novel approach, we acquire multi-channel MRI head coil array data and then apply inverse reconstruction methods to obtain volumetric fMRI estimates of blood oxygenation level dependent (BOLD) contrast at unprecedented whole-brain acquisition rates of 100 ms per sample. We call this combination of techniques magnetic resonance Inverse Imaging (InI), a method that provides estimates of dynamic spatially-resolved signal change that can be used to construct statistical maps of task-related brain activity. We demonstrate the sensitivity and inter-subject reliability of volumetric InI using an event-related design to probe the hemodynamic signal modulations in primary visual cortex. Robust results from both single subject and group analyses demonstrate the sensitivity and feasibility of using volumetric InI in high temporal resolution investigations of human brain function.

Hemodynamically based fMRI (Belliveau et al., 1991; Belliveau et al., 1990; Kwong et al., 1992; Ogawa et al., 1990) is typically limited to a temporal sampling period of two to four seconds if whole brain coverage is desired. Most fMRI data acquisition methods employ an EPI technique that utilizes many phase encoding steps and multiple read-out gradients. Consequently, this reliance on gradient encoding results in long image acquisition times and relatively loud acoustic noise related to the requisite rapid gradient switching. Here, we demonstrate the use of a novel volumetric imaging method, called Inverse Imaging (InI), which uses minimal phase encoding to achieve an order-of-magnitude improvement in BOLD-contrast temporal resolution. Its minimal dependence on encoding gradients allows extremely short image acquisition times, with an associated trade-off involving somewhat reduced and spatially-varying spatial resolution.

The temporal resolution of MRI is limited by the time required to traverse k-space during signal acquisition. The collection of volumetric MRI data continues until the completion of k-space traversal in multiple 2D k-spaces or in a single 3D k-space. Classical gradient-echo or spin-echo image acquisition methods collect data from one k-space line during each excitation. Thus the total acquisition time for traditional 3D MRI data acquisition is the product of the number of slices and the number of phase encoding steps. In contrast to gradient-echo or spin-echo imaging, both echo-planar imaging (EPI) (Mansfield, 1977) and spiral imaging (Blum et al., 1987) utilize fast gradient switching to achieve 2D k-space traversal in a single RF excitation. With current state-of-the-art EPI or spiral imaging techniques, one 2D single slice image can be collected in approximately 80 ms, allowing whole head coverage with 3 mm isotropic resolution in two to four seconds. Small improvements in temporal resolution can be achieved by optimizing k-space sampling schemes and reconstruction methods: e.g., instead of completing the k-space traversal for every measurement, MRI data acquisition can be accelerated by coordinated alterations of in k-space trajectories and their associated image reconstruction algorithms, as in partial-k space sampling (McGibney et al., 1993). Alternatively, *a priori* information-based methods can improve the temporal resolution of MR dynamic measurements (Tsao et al., 2001).

Recently, parallel imaging methods have been introduced to reconstruct images using spatial information derived simultaneously from multiple coil array channels The techniques employed include k-space SMASH (Sodickson and Manning, 1997), k-space GRAPPA (Griswold et al., 2002) and image domain SENSE (Pruessmann et al., 1999), all of which share a similar theoretical background (Sodickson and McKenzie, 2001). While parallel MRI can accelerate data acquisition rates by reducing total k-space traversal at the cost of reduced signal-to-noise ratio (SNR), the resulting net acceleration rate is limited both by the number of array coils and the specific phase-encoding scheme employed.

Prior information can be incorporated by combining EPI with parallel MR imaging (Golay et al., 2000; Preibisch et al., 2003; Schmidt et al., 2005; Weiger et al., 2002), resulting in fMRI detection sensitivity improvements with sensitivity encoded parallel MRI techniques (Lin et al., 2005b). Prior-informed parallel MRI has been explored using a fixed regularization parameter with empirical singular value decomposition truncation (King, 2001; Sodickson, 2000). Incorporation of prior information can suppress noise amplification in parallel MRI reconstruction (Lin et al., 2002; Lin et al., 2005b; Lin et al., 2004; Tsao et al., 2002) and traditional parallel MRI has been used to solve under-determined ill-posed problems (Katscher and Manke, 2002; Tsao et al., 2003). However, either only minor acceleration has been achieved (4-fold acceleration using a 2-channel array in cardiac imaging) (Katscher and Manke, 2002), or the reconstruction process has depended on incorporation of low-resolution prior image information (Tsao et al., 2003; Tsao et al., 2005).

More extreme accelerations in MRI acquisition rates have been achieved by reconstructing each image from a single echo. For example, single-echo-acquisition (SEA) was achieved using a dedicated 64-channel linear planar array that eliminated phase encoding, instead using the spatial information obtained from an array of long and parallel coils. This planar pair element design proved to be crucial for achieving well-localized field sensitivity patterns (McDougall and Wright, 2005). In other work, Hennig developed the one-voxel-one-coil (OVOC) MR-encephalography technique, obtaining a reconstructed image by computing the product of a full FOV reference scan and the accelerated acquisition scan where traditional phase and frequency encoding can be selectively omitted. This approach uses simultaneous multi-channel acquisition with multiple small receiver coils sampled such that the signal received by each coil is read out separately. The effective voxel size observed by each receiver channel is determined by the sensitive volume of the corresponding coil element and the source spatial distribution is estimated by constrained reconstruction using images from each separate coil as references (Hennig et al., 2007). A similar reconstruction algorithm termed HYPR was also proposed in the context of MR angiography (Mistretta et al., 2006). Nevertheless, none of these approaches explicitly formulate the relationship between the spatial information contained in the different channels of a RF coil array with full gradient encoding or with minimal gradient encoding. Nor do they provide algorithms to estimate the significance of task-related signal changes that would allow dynamic statistical inferences to be made from a highly temporally resolved data set.

Mathematically, the maximum acceleration possible with parallel MRI acquisition is limited by the available independent spatial information encoded by the coil array elements. This limit manifests itself as a problem in solving an over-determined linear system. Increasing the number of channels can thus increase MRI sampling rates. To this end, dense head coil arrays consisting of 16 (Bodurka et al., 2004; de Zwart et al., 2002; de Zwart et al., 2004), 23, 32 and 90 elements (Wiggins et al., 2005a; Wiggins et al., 2005b) have been constructed in support of a range of parallel acquisition applications. In addition, a dedicated 64-channel linear planar array has been developed to achieve 64-fold acceleration (McDougall and Wright, 2005). Notably, the geometric configuration of our 32-channel head array is remarkably similar to that used for electrode and super-conducting quantum interference (SQUID) sensor arrays in modern EEG and MEG systems (Hamalainen et al., 1993). While the MEG sensors detect magnetic fields generated by neural currents (Hamalainen et al., 1993), MRI detects oscillating electromagnetic fields from magnetization precession (Haacke, 1999). In addition, while MEG derives all of its spatial information from the geometry of the detectors, current accelerated MRI methods still rely heavily on gradient encoding.

We have generalized parallel MRI reconstruction techniques to exceed the limitations encountered when utilizing an under-determined linear system by introducing single-shot volumetric MR Inverse Imaging (InI), an approach that employs an over-determined linear system in order to achieve dramatically reduced acquisition times. We demonstrate the use of single-shot volumetric InI in supporting dynamic spatially-resolved statistical inference in a functional neuroimaging experiment. Inspired by MEG and EEG source localization techniques, we use a generalization of prior-informed parallel MRI (Lin et al., 2005b; Lin et al., 2004) and an adaptation of MEG reconstruction methods to MRI, to reduce the whole-brain sampling time by minimizing the k-space traversal time. Rather than relying on gradient encoding, InI derives spatial information by solving the inverse problem utilizing information from all array channels. Thus, given the constraint imposed by the need to use echo times (TE) that are optimal for BOLD-contrast, InI can complete k-space traversal and acquire sufficient data for whole-brain image reconstruction in under 100 ms. Although we have previously shown the feasibility of a 2D InI implementation (Lin et al., 2006), we now demonstrate its application to functional imaging studies employing 3D whole-brain coverage and event-related designs (Rosen et al., 1998). Event-related fMRI is a widely utilized neuroimaging method to study not only spatial but also temporal activity of hemodynamic changes secondary to the neuronal events. Compared to the classical block-design fMRI, the timing information available in the event-related fMRI allows for the study of both transient and steady states of cerebrovascular responses. This experimental technique mitigates the difficulty of potential bias originated from the contexts or the history of previous stimuli events. Event-related fMRI also enables the analysis of data using post-hoc categorization (Wagner et al., 1998). In some experimental designs, such as “odd-ball” experiments, can only be implemented using event-related fMRI rather than block-design (Friston, 2007). All reasons described above encouraged us to study the feasibility of 3D InI acquisitions and reconstructions using event-related fMRI design.

The principal novelty of our method is its combination of dense coil arrays with a linear estimation approach, allowing the transition from a largely gradient encoded to a largely detector encoded image, thereby achieving an order-of-magnitude speedup in the frame rate of dynamic whole-brain MRI. In the following sections, we present the data acquisition strategy and mathematical algorithms underlying InI, which allow extension of the technique to include event-related, 3D functional imaging designs with whole-brain coverage. We next demonstrate the technique’s capabilities in measuring the spatiotemporal properties of the hemodynamic response to brief visual stimuli using a 100 ms temporal sampling rate on a 3T scanner with a 32-channel coil array. In comparison with conventional EPI, 3D InI exhibits comparable sensitivity and adequate spatial resolving power in detecting visual task-related activity when performance is examined at both the single subject and group levels.

Five healthy participants were recruited for the study. Informed consent approved by the Institutional Review Board was obtained from each participant.

Our participants were asked to maintain fixation at the center of a tangent screen while viewing a high-contrast visual checkerboard reversing at 8 Hz. The checkerboard subtended 20° of visual angle and was generated from 24 evenly distributed wedges (15° each) and eight concentric rings of equal width. The stimuli were generated using the Psychtoolbox (Brainard, 1997; Pelli, 1997) and presented using Matlab (Mathworks, Natick, MA). The checkerboard stimuli were presented for 500 ms duration and the onset of each checkerboard presentation epoch was randomized with a uniform distribution of inter-stimulus intervals varying from 3 to 16 s. Thirty-two visual stimulation epochs were presented during each of four 240 s runs, resulting in a total of 128 visual stimulation epochs per participant.

MRI data were collected on a 3T MRI scanner (Tim Trio, Siemens Medical Solutions, Erlangen, Germany). We used a custom-built 32-channel head array receive coil (Wiggins et al., 2006) and a body transmit coil. The array consisted of 32 circular surface coils tessellated evenly to cover the whole brain. Functional imaging included standard BOLD-contrast imaging using conventional EPI techniques as well as volumetric InI.

Using 3D InI, each acquired volume was obtained at a particular time instant by the combination of EPI frequency encoding along the inferior-superior direction and phase encoding along the anterior-posterior direction. Figure 1 shows a schematic diagram of the spatial encoding procedure used to combine EPI and InI. The spatial information in the left-right direction was recovered during the InI reconstruction computation.

The 3D InI spatial encoding scheme implemented using a 32-channel array of coils. One thick sagittal slice was spatially encoded using an EPI sequence, allowing resolution of the spatial information in both the anterior-posterior (A-P) and superior-inferior **...**

InI reconstruction requires collection of a reference scan that provides information about the entire 3D volume. With this reference scan data, also called the forward operator, accelerated acquisition is enabled by replacing time consuming spatial encoding, dependent upon gradient switching, with an alternative approach utilizing an image reconstruction algorithm.

The reference scan was implemented using a single-slice echo-planar imaging (EPI) readout approach. Specifically, we excited one thick sagittal slab covering the entire brain (FOV 256 mm × 256 mm × 256 mm; 64 × 64 × 64 image matrix), setting the flip angle to the Ernst angle of 30°. Partition phase encoding was used to obtain the spatial information along the left-right axis (inter-aural line). The EPI readout had frequency and phase encoding along the superior-inferior and anterior-posterior axes respectively. We used 100 ms TR, 30 ms TE, 2604 Hz bandwidth and a 10 s total acquisition time for the reference scan.

InI functional data acquisition used the same volume prescription, TR, TE, flip angle, and bandwidth as the reference scan. The principal difference was that the partition phase encoding was removed. The full volume was excited like in the reference scan, and the spins were spatially encoded by a single-slice EPI acquisition. This resulted in a projection image along the left-right direction. The InI reconstruction algorithm, described in the next section, was then used to estimate the spatial information along the *x–y* axis. In each run, we collected 2,400 measurements after collecting 32 measurements in order to reach the longitudinal magnetization steady state. A total of four runs of data were acquired from each participant.

To validate the InI functional results, conventional EPI data, using conventional BOLD-contrast detection methods, was also collected using identical stimulus and presentation paradigm timing. EPI functional data acquisition used TE = 30 ms, TR = 2 s, flip angle = 90°, FOV = 220 mm, 24 slices, slice thickness = 4 mm with 20% gap).

In addition to the EPI and InI functional data, anatomical MRI data for each participant were obtained in the same session using a high-resolution T1-weighted 3D sequence (MPRAGE, TR/TE/flip = 2530 ms/3.49 ms/7°, partition thickness = 1.33 mm, matrix = 256 × 256, 128 partitions, FOV = 21 cm × 21 cm). Using these data, the location of the gray-white matter boundary was estimated with an automatic segmentation algorithm to yield a triangulated mesh model with approximate 340,000 vertices (Dale et al., 1999; Fischl et al., 2001; Fischl et al., 1999). The spatial registration between the EPI data or InI reconstructed data and the anatomical data was done by using SPM5 (http://www.fil.ion.ucl.ac.uk/spm/), estimating a 12-parameter rigid body transformation between the EPI mean image or volumetric InI reference scan and MPRAGE anatomical study. Note that the actually accelerated InI acquisitions were not spatially registered using this tool. The registration matrix was subsequently applied to each time instant of the EPI or InI hemodynamic estimates, to transform the neural activity estimates for each functional run to an inflated cortical surface space (Dale et al., 1999; Fischl et al., 1999). The transformed results were also spatially smoothed with a 3D Gaussian kernel with 10 mm full-width-half-maximum (FWHM).

Both reference scans and InI scans use EPI read-out, which demonstrate classical Nyquist ghost artifact in the phase-encoding direction (anterior-posterior direction in our implementation). To correct this artifact and to make consistent measurements between reference and InI scans, we performed the same phase correction algorithm in each reference (with partition encoding) and InI (without partition encoding) image at all channels of the coil array respectively. Specifically, 64 interleave EPI read-outs were acquired in the frequency encoding direction (superior-inferior direction) without phase encoding blips to estimate the k-space shift between even and odd echoes along the frequency encoding direction. Such an echo-shifting was then corrected by appending a compensation phase to the even echoes to avoid N/2 Nyquist ghost. The InI acquisition and reference scans were processed using 2D and 3D fast Fourier transformations from the k-space domain to the image domain, respectively. The reference scan in each channel of the coil array was synthetically averaged across partitions to simulate the InI acquisitions by making projection images along the partition encoding direction in each channel of the coil array, such that ${d}_{i}^{\mathit{\text{SIM}}}(\stackrel{\rightharpoonup}{r}),$ the simulated InI acquisition at location and channel *i* , is calculated as

$${d}_{i}^{\mathit{\text{SIM}}}(\stackrel{\rightharpoonup}{r})={\displaystyle \sum _{\stackrel{\rightharpoonup}{\rho}}{A}_{i}(\stackrel{\rightharpoonup}{\rho})}.$$

[Eq. 1]

. Here represents the spatial location indices across different partition phase encoding steps with the same frequency and phase encoding numbers, indicated by the spatial index , and *A _{i}* () represents the reference scan image from location and channel

$${\theta}_{i}(t)=\angle ({\displaystyle \sum _{\stackrel{\rightharpoonup}{r}}({d}_{i}(\stackrel{\rightharpoonup}{r},t)/{d}_{i}^{\mathit{\text{SIM}}}(\stackrel{\rightharpoonup}{r}))}),$$

[Eq. 2]

where *d _{i}* (,

After pre-processing of the whole time series in the InI projection acquisitions, now we tried to estimate the hemodynamic response function (HRF) in each projection images across all channels of the coil array. This effort reduces the size of data in time-domain dramatically without affecting the subsequent image-domain image reconstruction. From a list of the stimulus onset times, we constructed a vector with entries containing a one indicate the occurrence of stimuli during each 500 ms stimulation period, and all other entries contain zeros. A contrast matrix *D* was constructed from the convolution between the vector and the hemodynamic response function *H* ,

$$D=\stackrel{\rightharpoonup}{p}\otimes H\phantom{\rule{thinmathspace}{0ex}},$$

[Eq. 3]

, where denotes the linear convolution.

We used a finite-impulse response (FIR) basis for the HRF and thus *H* is a discrete delta function with different delay indicating the HRF at a specific time instant. Specifically, *H* modeled a HRF of 30 second duration with six second pre-stimulus baseline of 100 ms TR. Thus *H* is an identity matrix of dimension 300 (i.e., *n _{h}* = 300). In addition to the phase drift

$${\beta}_{i}(\stackrel{\rightharpoonup}{r},\tau )={({D}^{H}D)}^{-1}{D}^{H}{d}_{i}(\stackrel{\rightharpoonup}{r},t),\tau =1\dots {n}_{h},$$

[Eq. 4]

where the superscript *H* indicates the transpose and complex conjugate. HRF estimation, *h _{i}*(,

The minimum-norm estimate (MNE) provided the reconstructed image *x*(,*τ*) of 3D volumetric index and time index *τ* , which can be expressed as

$$\begin{array}{c}\begin{array}{cc}x(\stackrel{\rightharpoonup}{\rho},\tau )\hfill & =A{(\stackrel{\rightharpoonup}{r})}^{H}{(A(\stackrel{\rightharpoonup}{r})A{(\stackrel{\rightharpoonup}{r})}^{H}+\lambda (\stackrel{\rightharpoonup}{r})C)}^{-1}h(\stackrel{\rightharpoonup}{r},\tau )\hfill \\ \hfill & ={W}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{r})h(\stackrel{\rightharpoonup}{r},\tau )\hfill \end{array}\hfill \end{array}$$

[Eq. 5]

where *C* is the noise covariance matrix of the array, λ() is the regularization parameter, and *W _{MNE}*() is the inverse operator for the InI location index .

$$\lambda (\stackrel{\rightharpoonup}{r})=\mathit{\text{Tr}}(C)/\mathit{\text{Tr}}(A{(\stackrel{\rightharpoonup}{r})}^{H}A{(\stackrel{\rightharpoonup}{r})}^{H})/{\mathit{\text{SNR}}}^{2}.$$

[Eq. 6]

Here *Tr*(•) indicates the trace of the matrix. This choice of regularization parameter followed our original 2D InI reconstruction algorithm and was inspired by the analogy between InI reconstruction (Lin et al., 2006) and distributed source modeling in MEG (Lin et al., 2005a). We used *SNR* = 5 in this study. Such a choice of SNR was derived from our previous InI study (Lin et al., 2006). Sensitivity of MNE reconstruction in InI was investigated in our previous study, which showed that the regularization parameter affects InI reconstructions moderately if it changes within ±10 fold. Our choice of regularization parameter should be within that range.

To allow statistical inference from the results of the InI time series reconstruction, we estimated the noise in the baseline by applying the MNE inverse operator to the baseline InI data. Then, dynamic statistical parametric maps were derived as the time-point by time-point ratio between InI reconstruction and this baseline noise estimates, given by

$$\begin{array}{cc}t(\stackrel{\rightharpoonup}{\rho},\tau )\hfill & =x(\stackrel{\rightharpoonup}{\rho},\tau )/\mathit{\text{diag}}({W}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{r}){\mathit{\text{CW}}}_{\mathit{\text{MNE}}}{(\stackrel{\rightharpoonup}{r})}^{H})\hfill \\ \hfill & ={W}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho})d(\stackrel{\rightharpoonup}{\rho},\tau )/\mathit{\text{diag}}(.{W}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho})C.{W}_{\mathit{\text{MNE}}}{(\stackrel{\rightharpoonup}{\rho})}^{H}),\hfill \\ \hfill & ={W}_{\mathit{\text{MNE}}-\mathit{\text{dSPM}}}(\stackrel{\rightharpoonup}{\rho})d(\stackrel{\rightharpoonup}{\rho},\tau )\hfill \end{array}$$

[Eq. 7]

where *diag*(•) is the operator to construct a diagonal matrix from the input argument vector. Here *x*(,*τ*) represents estimated signal vector and *diag*(*W _{MNE}*()

The statistical modeling of the conventional EPI BOLD-contrast data was implemented using a general linear model (GLM) with a design matrix based on the same finite impulse response (FIR) model (c.f. Eq. 3 and Eq. 4). Using least square fitting, we estimated the strength of the FIR HDR at each voxel at each voxel and each time instant, and the residuals in GLM led to an estimate of baseline variability. Taken together, the *t*-statistics of the evoked hemodynamic response were calculated for each image voxel at each time instant.

Here we want to emphasize that InI image reconstruction had separate processing for HRF estimation in the projection images in all channels of the coil array (time domain) and restoration of 3D spatial information using MNE (spatial domain). However, EPI data had only time domain HRF estimation. In the time domain estimation of HRF, both InI and EPI data used the same finite-impulse-response (FIR) basis in the GLM estimation of HDR and thus the results are comparable. The use of FIR basis function has no prior assumption on the shape of HDR. This is different from the traditional use of parameterized models to estimate canonical HDR using, for example, Gamma functions. Compared to traditional EPI reconstruction, the only additional processing in InI reconstruction is the MNE in the spatial domain to convert the 2D projected images into 3D volumetric images.

To compare the similarity of InI and EPI reconstructions, InI dSPMs of *t*-statistic maps were averaged between 3–7 s after stimulus onset. The center of mass of this temporally averaged *t*-statistic map from the InI *T _{INI}*() and EPI

$$\begin{array}{c}{\stackrel{\rightharpoonup}{\rho}}_{\mathit{\text{InI}}}={\displaystyle \sum _{\stackrel{\rightharpoonup}{\rho}}{T}_{\mathit{\text{InI}}}(\stackrel{\rightharpoonup}{\rho})\stackrel{\rightharpoonup}{\rho}/{\displaystyle \sum _{\stackrel{\rightharpoonup}{\rho}}{T}_{\mathit{\text{InI}}}(\stackrel{\rightharpoonup}{\rho})}}\hfill \\ {\stackrel{\rightharpoonup}{\rho}}_{\mathit{\text{EPI}}}={\displaystyle \sum _{\stackrel{\rightharpoonup}{\rho}}{T}_{\mathit{\text{EPI}}}(\stackrel{\rightharpoonup}{\rho})\stackrel{\rightharpoonup}{\rho}/{\displaystyle \sum _{\stackrel{\rightharpoonup}{\rho}}{T}_{\mathit{\text{EPI}}}(\stackrel{\rightharpoonup}{\rho}).}}\hfill \\ \eta =\Vert {\stackrel{\rightharpoonup}{\rho}}_{\mathit{\text{InI}}}-{\stackrel{\rightharpoonup}{\rho}}_{\mathit{\text{EPI}}}\Vert \hfill \end{array}$$

[Eq. 8]

Here ‖•‖ denotes the Euclidean norm.

We performed numerical simulations to evaluate the spatial resolution and localization accuracy of our InI reconstructions. The reference data for the forward operator *A*() and noise covariance matrix *C* . The simulation procedure started from creating a source vector *x*(), with asset to unit activity and other locations set to zero. We then estimated the idealized measurements from all coil array channels by computing the product of the forward operator *A*() and *x*().

$$s(\stackrel{\rightharpoonup}{r})=A(\stackrel{\rightharpoonup}{\rho})x(\stackrel{\rightharpoonup}{\rho})$$

[Eq. 9]

We created 100 realizations of synthetic noise with spatial coloring according to the noise covariance matrix:

$$n(\stackrel{\rightharpoonup}{r})={U}_{c}{S}_{c}^{1/2}\epsilon (\stackrel{\rightharpoonup}{r})$$

[Eq. 10]

where *ε*() is the noise vector with complex values following a Gaussian distribution of zero mean and unit variance. *U _{c}* and

$$d(\stackrel{\rightharpoonup}{r})=s(\stackrel{\rightharpoonup}{r})+\frac{1}{\mathit{\text{SNR}}}\sqrt{\frac{\text{max}({\left|s(\stackrel{\rightharpoonup}{r})\right|}^{2})}{\mathit{\text{Tr}}({S}_{c})}}n(\stackrel{\rightharpoonup}{r}).$$

[Eq. 11]

Then we followed Eq. 5 to obtain the data covariance matrix, the inverse operator, and the noise normalized inverse operator:

$$\begin{array}{c}{\widehat{\stackrel{\rightharpoonup}{x}}}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho})={W}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho})d(\stackrel{\rightharpoonup}{r})\hfill \\ {\widehat{\stackrel{\rightharpoonup}{x}}}_{\mathit{\text{MNE-dSPM}}}(\stackrel{\rightharpoonup}{\rho})={W}_{\mathit{\text{MNE-dSPM}}}(\stackrel{\rightharpoonup}{\rho})d(\stackrel{\rightharpoonup}{r})\hfill \end{array}.$$

[Eq. 12]

. The InI reconstruction obtained with this procedure is equivalent to the point spread function of the simulated source *x*() . Both
${\widehat{\stackrel{\rightharpoonup}{x}}}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho})\phantom{\rule{thinmathspace}{0ex}}and\phantom{\rule{thinmathspace}{0ex}}{\widehat{\stackrel{\rightharpoonup}{x}}}_{\mathit{\text{MNE-dSPM}}}(\stackrel{\rightharpoonup}{\rho})$ were scaled to a maximum of 1.

Similar to procedures used in MEG/EEG source analysis (Dale et al., 2000; Liu et al., 1998; Liu et al., 2002), we estimated the average point spread function (aPSF) at each location to quantify the spatial spread of the reconstruction:

$$\begin{array}{c}\hfill {\mathit{\text{aPSF}}}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho})=\frac{{\displaystyle \sum _{i,i\ne \rho}|{\mathbf{d}}_{i}(\stackrel{\rightharpoonup}{\rho})|{x}_{\mathit{\text{MNE}}}^{i}}}{l}\hfill \\ \hfill {\mathit{\text{aPSF}}}_{\mathit{\text{MNE-dSPM}}}(\stackrel{\rightharpoonup}{\rho})=\frac{{\displaystyle \sum _{i,i\ne \rho}|{\mathbf{d}}_{i}(\stackrel{\rightharpoonup}{\rho})|{x}_{\text{MNE-dSPM}}^{i}}}{l}\hfill \end{array},$$

[Eq. 13]

where |**d**_{i}()| indicates the distance between source location *i* and source location . ${x}_{\mathit{\text{MNE}}}^{i}(\text{or}\phantom{\rule{thinmathspace}{0ex}}{x}_{\mathit{\text{MNE}}-\mathit{\text{dSPM}}}^{i})$ represents vector entries in the InI reconstruction ${\widehat{\stackrel{\rightharpoonup}{x}}}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho})\phantom{\rule{thinmathspace}{0ex}}(\text{or}\phantom{\rule{thinmathspace}{0ex}}{\widehat{\stackrel{\rightharpoonup}{x}}}_{\mathit{\text{MNE-dSPM}}}(\stackrel{\rightharpoonup}{\rho}))$ exceeding 0.5 and *l* is the number of voxels to be spatially resolved by the InI reconstructions. This procedure allows estimation of the full-width-half-maximum (FWHM) of the point spread function. A 3D map of the spatial distribution of the average point spread function, aPSF, for either MNE or MNE dSPM estimates can be obtained by repeating the calculation across the whole source space, the 256 × 256 × 256 mm FOV and 4 mm isotropic spatial resolution.

Since InI is an intrinsically ill-posed inverse problem, the reconstructed image may not be the original spatial distribution of spins contributing to the actual measurements. Thus the analysis of localization accuracy on the discrepancy between the reconstructed sources and the original sources is desired. Quantification of localization accuracy was done by calculating the shift between the center of mass of InI reconstruction and the simulated source:

$$\begin{array}{c}\hfill {\mathit{\text{SHIFT}}}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho})=\left|\left({\displaystyle \sum _{i,i\ne \rho}{\widehat{\stackrel{\rightharpoonup}{x}}}_{\mathit{\text{MNE}}}(\stackrel{\rightharpoonup}{\rho}){x}_{\mathit{\text{MNE}}}^{i}}\right)-\stackrel{\rightharpoonup}{\rho}\right|\hfill \\ \hfill {\mathit{\text{SHIFT}}}_{\mathit{\text{MNE-dSPM}}}(\stackrel{\rightharpoonup}{\rho})=\left|\left({\displaystyle \sum _{i,i\ne \rho}{\widehat{\stackrel{\rightharpoonup}{x}}}_{\mathit{\text{MNE-dSPM}}}(\stackrel{\rightharpoonup}{\rho}){x}_{\mathit{\text{MNE-dSPM}}}^{i}}\right)-\stackrel{\rightharpoonup}{\rho}\right|\hfill \end{array}.$$

[Eq. 14]

A 3D SHIFT metric map for MNE and MNE-dSPM was generated in each simulation. Since the inverse operators depend on the SNR (and the measurement data), the SNRs were parametrically varied from 0.1 to 100. Figure 2 illustrates the procedure of spatial resolution analysis using simulations. The implementation of the image reconstruction and statistical analysis procedures were done with Matlab.

Figure 3 shows the spatial distribution of aPSF at different SNRs (SNR = 0.5, 1, 5, and 10). We observed global reductions in the aPSF for MNE reconstructions at higher SNRs. In particular, deep brain regions show a larger aPSF metric. This matches the physical intuition that at center of the head coil the B1 fields from all channels are less spatially disparate as well as the SNR is lower than that at periphery and thus the spatial resolution becomes worse. Quantitatively, the MNE-dSPM inverse still has an average aPSF of 8.64 mm when SNR is higher than 1. At extremely high SNR (SNR >50), MNE can provide an excellent spatial resolution (average aPSF of 0.15 mm). To translate the aPSF into a measure of spatial resolution, the point spread function reported here should be spatially convolved with the nominal spatial resolution of the fully gradient encoded scan. Thus, for example, the average spatial resolution at SNR=5 is approximately 8.7 mm with a standard deviation of 6.5 mm (see Table 1).

The distribution of spatial resolution (quantified by average point spread function) of reconstructions at different SNRs (SNR = 0.5, 1, 5 and 10). Only central 24 axial slices of total 64 slices in a brain volume are shown here.

The SHIFT metrics derived from MNE were shown in Figure 4. MNE is characterized by sporadic high SHIFT metrics at low SNRs (SNR =0.5 and 1). Larger errors occurred at source locations in deep brain areas. On average, the localization accuracy is higher than 5 mm when SNR is higher than 1. Details of the aPSF and SHIFT metrics are listed in Table 1.

Raw traces of InI acquisitions from three channels close to the occipital lobe, parietal lobe, and frontal lobe were shown in Figure 5. Clear cardiac and respiratory fluctuations were observed in the phase and magnitude plots of the time series from each RF coil channel. Figure 5 also shows that acquired InI images indeed had significant phase drifting at respiratory fluctuation. This supports the requirement of phase correction in Eq. 2. To illustrate the quality and consistency of the InI reconstruction, we show the difference magnitude images between simulated InI acquisitions ${d}_{i}^{\mathit{\text{SIM}}}(\stackrel{\rightharpoonup}{r})$ from reference scan and the those images of the mean of actually accelerated InI acquisitions in Figure 6. The discrepancy between the simulated and actual InI acquisitions was quantified as the percentage error with respect to the simulated InI acquisitions in each channel of the coil array separately. We observed that all 32 channels show a discrepancy of less than 10%. Note that the number of channel is arbitrary and it does not imply spatial location of each RF coil. This implies that the reference scan suffices an accurate forward operator for the subsequent inverse operator derivation and volumetric image reconstruction.

Difference magnitude images between the simulated and actually accelerated InI acquisitions from fully gradient-encoded reference scan at 32 channels. The difference between two was quantified by the percentage error labeled in the individual simulated **...**

To illustrate the spatial resolution of InI reconstruction, we performed InI reconstruction using the average InI acquisitions across time. The results are shown in Figure 7, which includes magnitude images of InI reconstructions and sum-of-squares reference images across 32-channels of the coil array at central 12 contiguous axial slices. Compared to the reference images, InI reconstructions preserved some features, including boundary, cerebral hemisphere boundary, and the contours of the whole brain. However, localized image features and contrasts are different. In the context of fMRI, we are interested in the time-domain contrast-to-noise ratio, rather than the spatial-domain contrast-to-noise ratio. Also, using univariate general linear model, the analysis treats different image voxel separately. Thus we tolerated such spatial heterogeneity and proceeded with time-series analysis.

(Top) InI reconstructed images from actually accelerated scan using MNE inverse operator. (Bottom) InI reconstructed images from fully gradient encoded reference scan using sum-of-squares algorithm.

Figure 8 shows two series (2.0 to 2.9 s and 6.0 – 6.9 s post stimulus) of functional activity estimated by InI dSPM overlaid on the sum-of-squares images from the reference scans at eight axial slices from a representative participant. The InI dSPM *t*-statistic maps were rendered on an inflated brain surface where light gray indicates gyri and dark gray indicates sulci. We used a critical threshold of *t* = 2 (uncorrected *p*-value = 0.028). The maps show a progressive activation in response to the reversing checkerboard around the calcarine sulcus. Starting from the most posterior part of the occipital lobe, the activity spreads anteriorly, superiorly, and inferiorly. Peak activity was observed between three and five s after stimulus onset. This response started to decrease approximately six seconds after stimulus onset

Single subject results showing successive frames of InI MNE-dSPM *t*-values superimposed on the sum-of-squares images from the reference scan at eight axial slices. The critical threshold used is *t* = 2.0 (uncorrected *p*-value = 0.028)

The area showing a positive visual response in the first three seconds after stimulus onset was used to define a region of interest (ROI) in the primary visual cortex (V1). The average values and standard deviations of this ROI time course are shown in Figure 9. Without utilizing any specific model of the hemodynamic response, InI measures revealed a sharp BOLD-contrast signal peak at four s after stimulus onset. We also observed a post-stimulus undershoot between 10 and 24 s after stimulus onset. The pre-stimulus interval shows fluctuations around the baseline.

Single subject InI dSPM and the location of the primary visual cortex (V1) ROI (inset). The time course shows the average (dark blue) and the standard deviation (light blue vertical error bars) of the InI dSPM *t*-values within the V1 ROI. The average EPI **...**

To compare InI with conventional EPI reconstruction, we compared data from the same participants studied with both techniques. Figure 10 shows the *t*-statistic maps from both InI and EPI acquisitions using an inflated cortical surface model. We averaged the InI dSPM *t*-statistics between three and five second after stimulus onset in order to match the two-second TR of the EPI acquisition from the same V1 ROI. To compensate for the differing sensitivity of the two methods, critical thresholds were chosen as *t* < 2 (uncorrected p-value <0.028) and *t* < 4 (uncorrected p-value <10^{−4}) for InI and EPI respectively in order to show similar size and pattern of cortical activity.

Single subject medial and ventral views of thresholded InI and EPI activity estimates seen in response to a 500 ms duration reversing checkerboard. The critical thresholds were chosen as *t* < 2 (uncorrected *p*-value <0.028) and *t* < **...**

Differing detection sensitivity between InI and EPI were quantified by examining the *t*-value maxima, the *t*-value baseline standard deviation, and the *z*-score transform of the *t*-values in primary visual cortex of each participant (see Table 2). Generally, EPI was associated with higher peak *t*-values compared to InI. However, the baseline variability of the EPI *t*-values was also higher than InI. The overall InI *z*-score transform of the *t* -values is thus higher than EPI. We noticed that the standard deviation of baseline t-statistics does not equal to 1 in V1 ROI. The reason for such a result is because we calculated the averaged t-statistics within V1 ROI first and then we calculated the standard deviation. The averaging within the V1 ROI certainly improves the noise fluctuation because of averaging across voxels within the ROI. This explained why the reported standard deviation of the t-statistics in the baseline period is less than 1.

The *t*-value maxima, *t*-value baseline standard deviation, and the *z*-score transform of the *t*-values in primary visual cortex for EPI and InI reconstructions of the five participants. On average, InI is associated with higher z-scores due its corresponding **...**

The distances η between the centers of mass of the EPI and InI reconstructions from five participants are listed in Table 3. On average, the spatial localization of InI and EPI peak activities differs by 5.2 mm. This localization discrepancy was quite stable among all five participants. The maximum distance is 6.8 mm and minimal distance is 3.4 mm, both of which are within 2 mm difference compared to the average.

Figure 11 shows single frames of InI dSPM *t*-values averaged over five participants at 100 ms temporal resolution. The shape of this group average time series is very similar in character to the individual participant time series shown in Figure 9, demonstrating the stable performance of InI acquisition and reconstruction techniques. Notably, the negative *t*-values observed in the single participant, shown in Figure 9, disappeared in this group average, implying that the signal decrease, as evidenced by the negative *t*-values, may be an individual or temporally variable effect. The individual frames of this group average show progressive increasing activity starting at 2.7 s after the stimulus onset (critical threshold *t* < 2; uncorrected *p*-value <0.028). The signal returned to baseline approximately 6.0 s after stimulus onset.

Single frames of the InI dSPM *t*-values in visual cortex averaged across five participants shown from the medial aspect of the left hemisphere using an inflated brain surface model. The critical threshold was t < 2.0 (uncorrected p-value < **...**

The time course of the InI dSPM *t*-values from the group average are shown in Figure 12. The shape of the hemodynamic response in the group average is similar to that observed in individual participants (see Figure 3). We observed a reduced variability along the time course, potentially resulting from averaging multiple participants. The peak activity was found at 4.5 s after stimulus onset and the post-stimulus undershoot was observed between 10 s and 24 s after stimulus onset in the group average.

We have shown that single-shot volumetric InI methods can achieve an order-of-magnitude acceleration in hemodynamic response sampling by combining dense head coil array parallel data acquisition with distributed source modeling. Applying volumetric InI to study the visual system using an event-related design, we found that the method is both sensitive and reliable, as demonstrated in the individual participant and group average results. Localization accuracy was examined by comparing the InI and conventional EPI task-related activity in primary visual cortex, with the comparison revealing a good spatial match (approximately 5 mm) between the two methods. Our results in the visual cortex activation show faster BOLD response compared to conventional EPI. This is comparable to results from literature (Janz et al., 1997; Pfeuffer et al., 2002). In comparison to previous studies, InI methods exhibit two principal advantages. First, they allow collection of whole-brain volumetric data using single-shot EPI acquisition, while previously it was only possible to obtain one single 2D image at high temporal resolution (20 ms). Volumetric acquisition is more attractive to neuroscientists primarily interested in investigating whole-brain spatiotemporal activity patterns. By sampling the entire brain, volumetric InI can avoid the need for tedious manual slice prescription based on the prior anatomical scans required to identify target brain areas. Second, our approach combines volumetric InI and event-related fMRI design (Rosen et al., 1998), a commonly employed experimental approach in functional neuroimaging studies. We found that volumetric InI acquisition and reconstruction had sufficient sensitivity to reliably detect task-related activity in primary visual areas with high temporal resolution. We also found that this study we had reduced false positive activation compared to our first attempt of InI (Lin et al., 2006) (see Figure 8). This is because the current implementation has better matching between the reference scan and the InI scan: both scan used the same 3D RF excitation, frequency encoding, and phase encoding. The only difference is that the reference scan had partition encoding steps, but InI scans did not. The improved matching between the reference scans and the InI acquisitions, as quantified in Figure 6, helps to generate a better reconstruction with higher sensitivity and reduced false positive activation.

The difference of SNR and sensitivity of EPI and InI have been analytically and empirically reported in our results. Both methods used BOLD contrast. However, due to different image reconstruction algorithm and sampling rate, EPI and InI have different sensitivity and degree of freedom. Compared to EPI, InI has an additional step of solving an inverse problem. Both EPI and InI time series modeling and statistical inference used the same general linear model. Using an MNE inverse operator in InI reconstruction resulted in a “decreased” sensitivity (as reported in the reduced peak *t*-statistics) as well as “decreased” detection variability (as reported in the reduced baseline *t*-statistics variation). This is due to the spatial smoothing intrinsic to the minimum L-2 norm inverse operator. Such a smoothing effect was analyzed quantitatively in the average point spread function simulations Note that InI and EPI also have different degrees of freedom, since they are of different sampling rates. We are aware of the potential temporal correlation inside the acquisitions. However, a further detail analysis on such correlation is beyond the scope of this study and we will pursue it in the near future.

One limitation in BOLD-contrast fMRI detection sensitivity is physiological noise originating from cardiac and respiratory fluctuations (Kruger and Glover, 2001). It has been shown that physiological noise sources are the dominant limiting factor in high-field fMRI (Giove et al.). Traditional EPI takes approximately two to four seconds to acquire a full 3D brain volume and therefore lacks sufficient temporal resolution to resolve the physiological noise which then becomes aliased into the recorded time series. It will be possible to utilize the high temporal resolution of InI to improve the sensitivity of BOLD-contrast fMRI by acquiring single-shot multiple-echo datasets at 10 Hz sampling rates. Due to this rapid sampling rate, volumetric InI data can satisfy the Nyquist sampling criterion and thereby allowing more effective digital filtering strategies to isolate and reduce cardiac and respiratory sources of noise. These reductions in temporal noise will result in improved detection of the relevant effects of interest in the experimental design.

InI solves an ill-posed inverse problem in image reconstruction. The limited spatial resolution from all channels of the coil array may not be able to provide a stable and unique solution for image reconstruction. This deficiency in independent coil information leads to limits in spatial resolution. Among all inverse solution alternatives, we want to bring a caution on the nominal spatial resolution and the actual image resolution. For example, using a classical equivalent current dipole (ECD), the spatial resolution can be regarded as infinitely high. However, this does not adequately reflect the physical resolution set by the measurement. The reconstructed InI images using MNE algorithm preserved some low-resolution image figures (Figure 7). This matches to the well-known spatial smoothing effect in a minimum- norm solution. Previously, we estimated that 2D InI usingin MNE reconstruction has a spatial resolution ranging between 5 mm and 16 mm, depending on the measurement SNR and the coil geometry (Lin et al., 2006). From our results, we can see that volumetric InI provides reasonable spatial accuracy relative to conventional EPI acquisition. However, a more systematic and extensive study of all factors modulating the spatial resolution, including SNR, coil geometry, field strength, and regularization parameters, will be required in order to quantify the factors influencing spatial variations in resolving power.

The spatial resolution is spatially-varying because the amount of independent spatial information from all channels in a coil array is varies among the image voxels within the field-of-view. Notably, in peripheral cortical areas the spatial resolution is higher and at the center of the brain the spatial resolution is lower (Lin et al., 2006). In volumetric InI, anisotropic spatial resolution effects only appear in the InI dimension (L-R direction in this study), while the other two spatial dimensions (A-P and S-I directions) still retain isotropic spatial resolution because gradient encoding is used. To improve the spatial resolution, there are two alternatives. First, increasing the number of channels in a coil array can provide more independent spatial information. However, the benefit of increasing channels will reach a plateau as the consequence of electromagnetic theoretical limitations (Ohliger et al., 2003; Wiesinger et al., 2004). In addition, at high field more independent spatial information can be obtained from the same geometry of a coil array as the consequence of a shorter wavelength. This implies that higher spatial resolution can be obtained at 7T, using the existing 32-channel coil array geometry and MNE reconstruction algorithm. Second, we can ameliorate the spatial blurring using other inverse reconstruction kernels, as discussed in the following section.

The volumetric InI technique solves an ill-posed inverse problem to obtain reconstructed images, using the MNE solution to estimate the spatial distribution of task-related activity. In MEG and EEG distributed source modeling, spatial filters using linear constraint minimal variance (LCMV) beamformers have been also investigated extensively (Hillebrand and Barnes, 2005; Robinson, 2004; Sekihara et al., 2002; Van Veen et al., 1997). In contrast to MNE, the LCMV approach minimizes the point spread function of the inverse operator, resulting in LCMV results that are more focal compared to MNE. However, one disadvantage of LCMV is its inefficient detection and separation of coherent sources (Van Veen et al., 1997). Interestingly, it has been shown that MNE and LCMV are mathematically related to each other; the difference residing in the fact that MNE is model-driven while LCMV is data-driven during the estimation of the data covariance matrix (Mosher et al., 2003).

Under the interpretation of Bayesian estimation theory, both LCMV and MNE use L-2 norm prior models to solve an inverse problem. Thus it may be possible to replace the L-2 norm model by an L-1 norm model to achieve a narrower point spread function in InI image reconstruction. The reason that the L-1 norm model has a higher spatial resolution (a narrower point spread function) is that the probabilistic distribution of L-1 norm follows a bi-exponential distribution and the L-2 norm model follows a normal distribution. Moreover, the bi-exponential distribution has a more concentrated probability around the mean compared to the normal distribution (Uutela et al., 1999). L-1 norm models have been previously studied in the context of MEG and EEG source localization under the term of “minimum-current estimate” MCE (Lin et al., 2005a; Matsuura and Okabe, 1995; Uutela et al., 1999). In future work, we will investigate the comparative spatial localization accuracy of the LCMV inverse and MCE approaches using InI data.

In this study we employed a single-shot volumetric InI acquisition with an EPI readout. Thus, the reconstructed projection image contains the expected EPI artifacts, including intra-voxel signal loss due to spatially inhomogeneous susceptibility distribution and geometrical distortion along the phase encoding direction due to its intrinsically slower bandwidth. Correction of these artifacts has been studied extensively. For example, to mitigate these artifacts, we can use field mapping to investigate the spatial distribution of the off-resonance effects and then use this information to reduce the susceptibility artifacts (Chen et al., 2006; Chen and Wyrwicz, 1999; Zeng and Constable, 2002). It is also possible to use parallel imaging techniques with EPI acquisitions to limit geometric distortion (Weiger et al., 2002), by systematically skipping multiple integer lines in the continuous sampling of different phase-encoding lines and then reconstructing the skipped phase encoding lines using spatial information embedded inside different array channels. Thus the effective bandwidth in the reconstructed image will be wider and will thereby reduce the distortion. Note that here the spatial information from the coil array channels is in an orthogonal direction between EPI phase encoding (anterior-posterior axis) and InI encoding (left-right axis). Thus the implementation of this approach will not reduce the available spatial information in the InI reconstructions, while SNR loss is a price that must be paid for reduced data sample numbers and changes in the parallel MRI reconstruction geometrical factor (g-factor).

The brain is a highly dynamic system. Consequently, detection and localization of static patterns of regional functional specialization will not be sufficient to fully understand the neural mechanisms underlying complex behavior. Nevertheless, current fMRI research mainly relies on methods that are most sensitive to these stationary patterns of task-related activity. Volumetric InI allows functional imaging measurements at high sampling rates over extended brain volumes using minimal gradient encoding. To our knowledge, this is the most rapid whole-brain fMRI achieved to date and this relatively high sampling rate may, in future work, enable the measurement of relative brain activity onset times and thereby provide a better understanding of the dynamic driving relationships among interacting neural subsystems across the whole brain. This goal cannot be achieved with conventional EPI methods due to their limited spatial coverage and temporal resolution.

Because of the relative sluggishness of the fMRI BOLD-contrast time course, it has been argued that improving the temporal resolution of fMRI might not reveal more information about brain dynamics. However, it is becoming increasingly clear that this popularly held view is incorrect, possibly resulting from a paucity of practical methods allowing ultrafast fMRI recordings in large brain volumes. For example, using limited spatial coverage, high temporal resolution fMRI (TR = 100 ms per slice or 400 ms per volume) can detect spatially distinct differences in neural activity onset times that can then be correlated with reaction time in order to localize the cortical areas active in the different phases of sensorimotor integration (Menon et al., 1998). In addition, fMRI experiments using jittered inter-stimulus-interval manipulations have detected millisecond-scale interactions among neural systems (Ogawa et al., 2000). However, no direct measurements of BOLD-contrast signals at high temporal resolution with spatial coverage of the entire brain have been achieved previously. Using volumetric InI, it may be possible to adapt time-resolved imaging techniques to experimental designs probing mutual interactions among a spatially distributed set of participating functionally specialized regions.

Motion correction is a critical process to improve detection power of fMRI. However, Motion correction on InI data is difficult because the acquisitions have only projection images. Such reduction of one dimensional spatial information poses a technical difficulty of utilizing post-processing method to do, for example, rigid body transformation, an intrinsic 3D processing. Thus currently InI may not have the optimal sensitivity due to motion artifacts. A potential method to mitigate this motion artifact issue is using navigator echoes during acquisition to correct motion effects on the flight at a moderate cost of temporal resolution (van der Kouwe et al., 2006). Even though volumetric InI allows dramatic improvement in sampling rates, it is still constrained by the need to use the optimal TE for detection of BOLD-contrast effects (approximately 30 ms at 3T). At higher field strengths, such as 7T, the optimal TE for BOLD-contrast would be 20 ms or less, allowing further acceleration of volumetric InI, possibly to 50 ms whole-brain sampling times. In addition, there are two potential approaches to mitigate the temporal resolution limitations: First, we may use different contrast mechanisms, such as steady-state free precession (SSFP), where TE is usually less than 5 ms (Miller et al., 2003; Miller et al., 2006). However, SSFP contrast is challenging at high field because of its high SAR from RF excitation. Second, we may use an echo-shifting technique to reduce the TR (Golay et al., 2000; Liu et al., 1993), as demonstrated in our previous 2D InI study (Lin et al., 2006). Both alternatives are capable of reducing the sampling time to around 20 ms. The resulting higher temporal resolution may allow study of relative cortical activity timing on an extraordinarily fine time scale, thereby facilitating study of the complex interactions among regionally specialized neural subsystems responsible for the mediation of complex behavior.

This work was supported by National Institutes of Health Grants R01 HD040712, R01 NS037462, P41 RR14075, the Mental Illness and Neuroscience Discovery Institute (MIND).

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

- Belliveau JW, Kennedy DN, Jr, McKinstry RC, Buchbinder BR, Weisskoff RM, Cohen MS, Vevea JM, Brady TJ, Rosen BR. Functional mapping of the human visual cortex by magnetic resonance imaging. Science. 1991;254:716–719. [PubMed]
- Belliveau JW, Rosen BR, Kantor HL, Rzedzian RR, Kennedy DN, McKinstry RC, Vevea JM, Cohen MS, Pykett IL, Brady TJ. Functional cerebral imaging by susceptibility-contrast NMR. Magn Reson Med. 1990;14:538–546. [PubMed]
- Blum MJ, Braun M, Rosenfeld D. Fast magnetic resonance imaging using spiral trajectories. Australas Phys Eng Sci Med. 1987;10:79–87. [PubMed]
- Bodurka J, Ledden PJ, van Gelderen P, Chu R, de Zwart JA, Morris D, Duyn JH. Scalable multichannel MRI data acquisition system. Magn Reson Med. 2004;51:165–171. [PubMed]
- Brainard DH. The Psychophysics Toolbox. Spat Vis. 1997;10:433–436. [PubMed]
- Chen NK, Oshio K, Panych LP. Application of k-space energy spectrum analysis to susceptibility field mapping and distortion correction in gradient-echo EPI. Neuroimage. 2006;31:609–622. [PubMed]
- Chen NK, Wyrwicz AM. Correction for EPI distortions using multi-echo gradient-echo imaging. Magn Reson Med. 1999;41:1206–1213. [PubMed]
- Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis. I. Segmentation and surface reconstruction. Neuroimage. 1999;9:179–194. [PubMed]
- Dale AM, Liu AK, Fischl BR, Buckner RL, Belliveau JW, Lewine JD, Halgren E. Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity. Neuron. 2000;26:55–67. [PubMed]
- de Zwart JA, Ledden PJ, Kellman P, van Gelderen P, Duyn JH. Design of a SENSE-optimized high-sensitivity MRI receive coil for brain imaging. Magn Reson Med. 2002;47:1218–1227. [PubMed]
- de Zwart JA, Ledden PJ, van Gelderen P, Bodurka J, Chu R, Duyn JH. Signal-to-noise ratio and parallel imaging performance of a 16-channel receive-only brain coil array at 3.0 Tesla. Magn Reson Med. 2004;51:22–26. [PubMed]
- Fischl B, Liu A, Dale AM. Automated manifold surgery: constructing geometrically accurate and topologically correct models of the human cerebral cortex. IEEE Trans Med Imaging. 2001;20:70–80. [PubMed]
- Fischl B, Sereno MI, Dale AM. Cortical surface-based analysis. II: Inflation, flattening, and a surface-based coordinate system. Neuroimage. 1999;9:195–207. [PubMed]
- Friston KJ. Statistical parametric mapping : the analysis of funtional brain images. 1st ed. Amsterdam ; Boston: Elsevier/Academic Press; 2007.
- Giove F, Mangia S, Bianciardi M, Garreffa G, Di Salle F, Morrone R, Maraviglia B. The physiology and metabolism of neuronal activation: in vivo studies by NMR and other methods. Magn Reson Imaging. 2003;21:1283–1293. [PubMed]
- Golay X, Pruessmann KP, Weiger M, Crelier GR, Folkers PJ, Kollias SS, Boesiger P. PRESTO-SENSE: an ultrafast whole-brain fMRI technique. Magn Reson Med. 2000;43:779–786. [PubMed]
- Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J, Kiefer B, Haase A. Generalized autocalibrating partially parallel acquisitions (GRAPPA) Magn Reson Med. 2002;47:1202–1210. [PubMed]
- Haacke EM. Magnetic resonance imaging : physical principles and sequence design. New York: J. Wiley & Sons; 1999.
- Hamalainen M, Hari R, Ilmoniemi R, Knuutila J, Lounasmaa O. Magnetoencephalography-theory, instrumentation, and application to noninvasive studies of the working human brain. Review of Modern Physics. 1993;65:413–497.
- Hennig J, Zhong K, Speck O. MR-Encephalography: Fast multi-channel monitoring of brain physiology with magnetic resonance. Neuroimage. 2007;34:212–219. [PubMed]
- Hillebrand A, Barnes GR. Beamformer analysis of MEG data. Int Rev Neurobiol. 2005;68:149–171. [PubMed]
- Janz C, Speck O, Hennig J. Time-resolved measurements of brain activation after a short visual stimulus: new results on the physiological mechanisms of the cortical response. NMR Biomed. 1997;10:222–229. [PubMed]
- Katscher U, Manke D. International Society for Magnetic Resonance in Medicine Tenth Scientific Meeting and Exhibition. Honolulu, Hawaii, USA: International Society for Magnetic Resonance in Medicine; 2002. Underdetermined SENSE using a-priori knowledge; p. 2396.
- King K. Proceedings of the 9th Annual Meeting of ISMRM. Glasgow, Scotland: International Society of Magnetic Resonance in Medicine; 2001. SENSE image quality improvement using matrix regularization; p. 1771.
- Kruger G, Glover GH. Physiological noise in oxygenation-sensitive magnetic resonance imaging. Magn Reson Med. 2001;46:631–637. [PubMed]
- Kwong KK, Belliveau JW, Chesler DA, Goldberg IE, Weisskoff RM, Poncelet BP, Kennedy DN, Hoppel BE, Cohen MS, Turner R, Cheng H, Brady TJ, Rosen BR. Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proc Natl Acad Sci U S A. 1992;89:5675–5679. [PubMed]
- Lin F-H, Kwong KK, Chen Y-J, Belliveau JW, Wald LL. International Society for Magnetic Resonance in Medicine Tenth Sceintific Meeting and Exhibition. Honolulu, Hawaii, USA: International Society for Magnetic Resonance in Medicine; 2002. Reconstruction of sensitivity encoded images using regularization and discrete time wavelet transform estimates of the coil maps; p. 2389.
- Lin FH, Belliveau JW, Dale AM, Hamalainen MS. Distributed current estimates using cortical orientation constraints. Hum Brain Mapp. 2005a [PubMed]
- Lin FH, Huang TY, Chen NK, Wang FN, Stufflebeam SM, Belliveau JW, Wald LL, Kwong KK. Functional MRI using regularized parallel imaging acquisition. Magn Reson Med. 2005b;54:343–353. [PMC free article] [PubMed]
- Lin FH, Kwong KK, Belliveau JW, Wald LL. Parallel imaging reconstruction using automatic regularization. Magn Reson Med. 2004;51:559–567. [PubMed]
- Lin FH, Wald LL, Ahlfors SP, Hamalainen MS, Kwong KK, Belliveau JW. Dynamic magnetic resonance inverse imaging of human brain function. Magn Reson Med. 2006;56:787–802. [PubMed]
- Liu AK, Belliveau JW, Dale AM. Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations. Proc Natl Acad Sci U S A. 1998;95:8945–8950. [PubMed]
- Liu AK, Dale AM, Belliveau JW. Monte Carlo simulation studies of EEG and MEG localization accuracy. Hum Brain Mapp. 2002;16:47–62. [PubMed]
- Liu G, Sobering G, Duyn J, Moonen CT. A functional MRI technique combining principles of echo-shifting with a train of observations (PRESTO) Magn Reson Med. 1993;30:764–768. [PubMed]
- Mansfield P. Multi-planar image formation using NMR spin echos. J. Physics. 1977;C10:L55–L58.
- Matsuura K, Okabe Y. Selective minimum-norm solution of the biomagnetic inverse problem. IEEE Trans Biomed Eng. 1995;42:608–615. [PubMed]
- McDougall MP, Wright SM. 64-channel array coil for single echo acquisition magnetic resonance imaging. Magn Reson Med. 2005;54:386–392. [PubMed]
- McGibney G, Smith MR, Nichols ST, Crawley A. Quantitative evaluation of several partial Fourier reconstruction algorithms used in MRI. Magn Reson Med. 1993;30:51–59. [PubMed]
- Menon RS, Luknowsky DC, Gati JS. Mental chronometry using latency-resolved functional MRI. Proc Natl Acad Sci U S A. 1998;95:10902–10907. [PubMed]
- Miller KL, Hargreaves BA, Lee J, Ress D, deCharms RC, Pauly JM. Functional brain imaging using a blood oxygenation sensitive steady state. Magn Reson Med. 2003;50:675–683. [PubMed]
- Miller KL, Smith SM, Jezzard P, Pauly JM. High-resolution FMRI at 1.5T using balanced SSFP. Magn Reson Med. 2006;55:161–170. [PubMed]
- Mistretta CA, Wieben O, Velikina J, Block W, Perry J, Wu Y, Johnson K, Wu Y. Highly constrained backprojection for time-resolved MRI. Magn Reson Med. 2006;55:30–40. [PMC free article] [PubMed]
- Mosher JC, Baillet S, Leahy RM. Equivalence of linear approaches in biomagnetic inverse solutions. 2003 IEEE workshop on statistical signal processing. St. Louis, Missouri: IEEE; 2003. pp. 294–297.
- Ogawa S, Lee TM, Kay AR, Tank DW. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc Natl Acad Sci U S A. 1990;87:9868–9872. [PubMed]
- Ogawa S, Lee TM, Stepnoski R, Chen W, Zhu XH, Ugurbil K. An approach to probe some neural systems interaction by functional MRI at neural time scale down to milliseconds. Proc Natl Acad Sci U S A. 2000;97:11026–11031. [PubMed]
- Ohliger MA, Grant AK, Sodickson DK. Ultimate intrinsic signal-to-noise ratio for parallel MRI: electromagnetic field considerations. Magn Reson Med. 2003;50:1018–1030. [PubMed]
- Pelli DG. The VideoToolbox software for visual psychophysics: transforming numbers into movies. Spat Vis. 1997;10:437–442. [PubMed]
- Pfeuffer J, Adriany G, Shmuel A, Yacoub E, Van De Moortele PF, Hu X, Ugurbil K. Perfusion-based high-resolution functional imaging in the human brain at 7 Tesla. Magn Reson Med. 2002;47:903–911. [PubMed]
- Preibisch C, Pilatus U, Bunke J, Hoogenraad F, Zanella F, Lanfermann H. Functional MRI using sensitivity-encoded echo planar imaging (SENSE-EPI) Neuroimage. 2003;19:412–421. [PubMed]
- Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med. 1999;42:952–962. [PubMed]
- Robinson SE. Localization of event-related activity by SAM(erf) Neurol Clin Neurophysiol 2004. 2004;2004:109. [PubMed]
- Rosen BR, Buckner RL, Dale AM. Event-related functional MRI: past, present, and future. Proc Natl Acad Sci U S A. 1998;95:773–780. [PubMed]
- Schmidt CF, Degonda N, Luechinger R, Henke K, Boesiger P. Sensitivity-encoded (SENSE) echo planar fMRI at 3T in the medial temporal lobe. Neuroimage. 2005;25:625–641. [PubMed]
- Sekihara K, Nagarajan SS, Poeppel D, Marantz A, Miyashita Y. Application of an MEG eigenspace beamformer to reconstructing spatio-temporal activities of neural sources. Hum Brain Mapp. 2002;15:199–215. [PubMed]
- Sodickson DK. Tailored SMASH image reconstructions for robust in vivo parallel MR imaging. Magn Reson Med. 2000;44:243–251. [PubMed]
- Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magn Reson Med. 1997;38:591–603. [PubMed]
- Sodickson DK, McKenzie CA. A generalized approach to parallel magnetic resonance imaging. Med Phys. 2001;28:1629–1643. [PubMed]
- Tsao J, Behnia B, Webb AG. Unifying linear prior-information-driven methods for accelerated image acquisition. Magn Reson Med. 2001;46:652–660. [PubMed]
- Tsao J, Boesiger P, Pruessmann KP. k-t BLAST and k-t SENSE: dynamic MRI with high frame rate exploiting spatiotemporal correlations. Magn Reson Med. 2003;50:1031–1042. [PubMed]
- Tsao J, Kozerke S, Boesiger P, Pruessmann KP. Optimizing spatiotemporal sampling for k-t BLAST and k-t SENSE: application to high-resolution real-time cardiac steady-state free precession. Magn Reson Med. 2005;53:1372–1382. [PubMed]
- Tsao J, Pruessmann K, Boesiger P. International Society for Magnetic Resonance in Medicine Tenth Scientific Meeting and Exhibition. Honolulu, Hawaii, USA: International Society of Magnetic Resonance in Medicine; 2002. Prior-information-enhanced dynamic imaging using single or multiple coils with k-t BLAST and k-t SENSE; p. 2369.
- Uutela K, Hamalainen M, Somersalo E. Visualization of magnetoencephalographic data using minimum current estimates. Neuroimage. 1999;10:173–180. [PubMed]
- van der Kouwe AJ, Benner T, Dale AM. Real-time rigid body motion correction and shimming using cloverleaf navigators. Magn Reson Med. 2006;56:1019–1032. [PubMed]
- Van Veen BD, van Drongelen W, Yuchtman M, Suzuki A. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans Biomed Eng. 1997;44:867–880. [PubMed]
- Wagner AD, Schacter DL, Rotte M, Koutstaal W, Maril A, Dale AM, Rosen BR, Buckner RL. Building memories: remembering and forgetting of verbal experiences as predicted by brain activity. Science. 1998;281:1188–1191. [PubMed]
- Weiger M, Pruessmann KP, Osterbauer R, Bornert P, Boesiger P, Jezzard P. Sensitivity-encoded single-shot spiral imaging for reduced susceptibility artifacts in BOLD fMRI. Magn Reson Med. 2002;48:860–866. [PubMed]
- Wiesinger F, Boesiger P, Pruessmann KP. Electrodynamics and ultimate SNR in parallel MR imaging. Magn Reson Med. 2004;52:376–390. [PubMed]
- Wiggins GC, Potthast A, Triantafyllou C, Lin F-H, Benner T, Wiggins CJ, Wald LL. International Society for Magnetic Resonance in Medicine Thirteenth Scientific Meeting and Exhibition. Miami, Florida, USA: International Society for Magnetic Resonance in Medicine; 2005a. A 96-channel MRI system with 23- and 90-channel phase array head coils at 1.5 Tesla; p. 671.
- Wiggins GC, Triantafyllou C, Potthast A, Reykowski A, Nittka M, Wald LL. International Society for Magnetic Resonance in Medicine Thirteenth Scientific Meeting and Exhibition. Miami, Florida, USA: International Society for Magnetic Resonance in Medicine; 2005b. A 32 Channel Receive-only Phased Array Head Coil for 3T with Novel Geodesic Tiling Geometry; p. 671.
- Wiggins GC, Triantafyllou C, Potthast A, Reykowski A, Nittka M, Wald LL. 32-channel 3 Tesla receive-only phased-array head coil with soccer-ball element geometry. Magn Reson Med. 2006;56:216–223. [PubMed]
- Zeng H, Constable RT. Image distortion correction in EPI: comparison of field mapping with point spread function mapping. Magn Reson Med. 2002;48:137–146. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |