NeuroImage 102 (2014) 861–874

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

Differentiating BOLD and non-BOLD signals in fMRI time series from anesthetized rats using multi-echo EPI at 11.7T Prantik Kundu a,b,⁎,1, Mathieu D. Santin c,1, Peter A. Bandettini a,d, Edward T. Bullmore b,e,f, Alexandra Petiet c a

Section on Functional Imaging Methods, National Institute of Mental Health, Bethesda, MD, USA Behavioural Clinical Neuroscience Institute, University of Cambridge, Cambridge, UK Center for Neuroimaging Research, Brain and Spine Institute, Paris, France d Functional MRI Core Facility, National Institute of Mental Health, Bethesda, MD, USA e NIHR Cambridge Biomedical Research Centre, Cambridgeshire Peterborough NHS Foundation Trust, UK f ImmunoPsychiatry, Alternative Discovery & Development, GlaxoSmithKline, Stevenage, UK b c

a r t i c l e

i n f o

Article history: Accepted 15 July 2014 Available online 24 July 2014

a b s t r a c t The study of spontaneous brain activity using fMRI is central to mapping brain networks. However, current fMRI methodology has limitations in the study of small animal brain organization using ultra-high field fMRI experiments, as imaging artifacts are difficult to control and the relationship between classical neuroanatomy and spontaneous functional BOLD activity is not fully established. Challenges are especially prevalent during the fMRI study of individual rodent brains, which could be instrumental to studies of disease progression and pharmacology. A recent advance in fMRI methodology enables unbiased, accurate, and comprehensive identification of functional BOLD signals by interfacing multi-echo (ME) fMRI acquisition, NMR signal decay analysis, and independent components analysis (ICA), in a procedure called ME-ICA. Here we present a pilot study on the suitability of ME-ICA for ultra high field animal fMRI studies of spontaneous brain activity under anesthesia. ME-ICA applied to 11.7 T fMRI data of rats first showed robust performance in automatic high dimensionality estimation and ICA decomposition, similar to that previously reported for 3.0 T human data. ME sequence optimization for 11.7 T indicated that 3 echoes, 0.5 mm isotropic voxel size and TR = 3 s was adequate for sensitive and specific BOLD signal acquisition. Next, in seeking optimal inhaled isoflurane anesthesia dosage, we report that progressive increase in anesthesia goes with concomitant decrease in statistical complexity of “global” functional activity, as measured by the number of BOLD components, or degrees of freedom (DOF). Finally, BOLD functional connectivity maps for individual rodents at the component level show that spontaneous BOLD activity follows classical neuroanatomy, and seed-based analysis shows plausible cortical–cortical and cortical–subcortical functional interactions. © 2014 Published by Elsevier Inc.

Introduction Resting state fMRI presents a remarkable opportunity to probe the functional organization of the brain based on the blood oxygenation level dependent (BOLD) response to spontaneous neuronal activity (Ogawa et al., 1990; Biswal et al., 2010). Many analyses of spontaneous BOLD signals are possible, including mapping the correlation of time series between regions of interest to determine inter-regional functional connectivity, and decomposition of data into spatiotemporal signal components using independent components analysis (ICA) (Biswal et al., 1995; Biswal and Ulmer, 1999). A fundamental limitation of fMRI that hinders these analyses is the presence of scanner and physiological artifact, which can account for greater than 50% of total fMRI time ⁎ Corresponding author. Fax: +1 301 402 1370. E-mail address: [email protected] (P. Kundu). 1 Authors contributed equally to this research.

http://dx.doi.org/10.1016/j.neuroimage.2014.07.025 1053-8119/© 2014 Published by Elsevier Inc.

series variance (Bianciardi et al., 2009). This makes functional BOLD signals a minority signal source, to an extent that varies widely across scanner platforms and experimental paradigms. Both spontaneous BOLD signals and fMRI artifacts have been evaluated in studies involving human fMRI, but small animal fMRI studies are comparatively few in number (Keilholz et al., 2004; Pawela et al., 2008; Bosshard et al., 2010; Hutchison et al., 2010; Jonckers et al., 2013). The rodent is a particularly interesting model system for studying functional connectivity, as it may provide new understanding of the effects of pharmacological agents or genetic manipulations (Borsook et al., 2006). However, the small size of the rat brain necessitates fMRI scanning using ultra high field MRI to provide adequate spatial and temporal resolution, which is problematic because complex artifacts are ordinarily encountered at these field strengths but poorly documented (Kalthoff et al., 2011). Regardless of these limitations, the rat model system is still very interesting and important for neurobiological studies based on BOLD MRI contrast.

862

P. Kundu et al. / NeuroImage 102 (2014) 861–874

A recent approach to analysis of spontaneous BOLD signals with a high potential for general application in non-human fMRI has been proposed by us (Kundu et al., 2011, 2013). This approach distinguishes functionally-related BOLD signals from non-BOLD artifacts using a physical approach to BOLD fMRI time series analysis, without depending on other information pertaining to neuroanatomy or specific artifact processes. Our approach firstly utilizes an fMRI pulse sequence that acquires T2* weighted images at multiple imaging contrasts corresponding to multiple echo times (TE), called multi-echo (ME) fMRI. Our method then exploits a property of BOLD signals across TEs: the percent signal change of a BOLD signal fluctuation scales linearly with TE. This is a consequence of the variation of signal magnitude change with TE due to the decay rate change (R2*), which results from oxygenation changes in response to functional activity (Menon et al., 1993). After a decomposition of multi-echo fMRI data into linear components by spatial independent component analysis (ICA; Hyvarinen, 1999), component signal changes are analyzed to assess linear TE-dependence (Kundu et al., 2011). ICA components are found to represent separated maps of distinct physical sources based on “independence” or sparsity of their latent statistical distributions (Daubechies et al., 2009). The TEdependence analysis of ICA component signal changes enables automated and physically principled identification of BOLD components based on their TE-dependence. Evaluation of these BOLD components based on their spatial localization and time course properties indicates functional networks. Many other components of ME-fMRI signals such as those from scanner and physiological artifacts do not show linear TE-dependent scaling of signal changes, and can in fact be described with a model for TE-independence. The TE-dependent BOLD signal and TE-independent non-BOLD artifact signal classes encompass the majority of fMRI signal variance (N95%), with the residual being constituted by thermal (normally distributed) noise. Statistics are then computed to express the overall component-level goodness of fit to BOLD and non-BOLD signal models: TE-dependence is expressed by the pseudo F-statistic κ, and TE-independence is expressed by the pseudo F-statistic ρ. κ and ρ are equivalently scaled. Functionally related BOLD signals are categorized as “high κ” and “low ρ”, and artifacts are the reverse. This approach, called multi-echo independent components analysis (ME-ICA), stably and automatically segregates BOLD and nonBOLD signal components in fMRI data, minimizing user-driven differentiation of signal and noise and avoiding arbitrary data filtering choices (such as bandpass filtering). Most any BOLD fMRI dataset imaged with a multi-echo (ME) sequence can be analyzed using ME-ICA. In this study we proposed that multi-echo fMRI of the rodent brain can allow the separation of BOLD and non-BOLD signals in a similar fashion to that previously demonstrated in human studies. However, unlike our prior human studies that were mostly conducted at high field (3 T), in this study we acquired ME-fMRI at 11.7 T to optimize sensitivity and resolution for rat imaging. Multi-echo fMRI pulse sequences acquire functional volumes using one excitation pulse and several readouts (separated by phase rewinders). After acquisition, a procedure involving essential preprocessing (e.g. slice-time correction and volume alignment), PCA dimensionality reduction and ICA decomposition, and BOLD vs. non-BOLD component classification was used to enable several analyses. BOLD ICA component maps and time courses were firstly evaluated to assess the spatial and temporal properties of rodent functional brain networks — the most common application of ICA in resting state fMRI studies to date. ME-ICA achieves ICA decompositions that explain a larger portion of total dataset variance than conventional ICA (over 90% vs. approximately 75%), based on an approach to dimensionality estimation involving principal components analysis (PCA) and computation of κ and ρ for PCA components, called ME-PCA. This approach sensitively and specifically identifies non-Gaussian (i.e. not thermal noise) components of the PCA space for subsequent ICA decomposition (i.e. unmixing), altogether using ICA to explain a greater proportion of dataset variance. Such a comprehensive ICA decomposition enables

further assessments of functional brain organization. Global functional brain variability can be assessed in terms of BOLD spatiotemporal degrees of freedom (DOF), which is approximated by the number of isolated BOLD independent components from high-dimensional decomposition as implemented in ME-ICA (Kundu et al., 2013). BOLD DOF can also be estimated after computing the eigenspectrum of a denoised BOLD signal time course dataset, which indicates the rank (in the linear algebra sense) of data. Based on the formulation of ICA as well as the interpretation of rank, BOLD DOF is indicative of the extent of unique information in a dataset, which is conceptually related to its compressibility or complexity (Comon, 1994). ME-ICA also enables seed-based functional connectivity analysis using a technique called independent component regression (ME-ICR). In this approach, a vector of component coefficients representing a given region of interest is used to compute correlation with coefficient vectors from other regions. This approach to computing functional connectivity is statistically wellconditioned because coefficients are independently distributed. BOLD DOF can then be used to estimate the standard error term for correlation, which in turn is used to standardize correlation values. Resulting ME-ICR correlation values are found to be more regularly distributed than ordinary time course correlation values, which are known to be ill-conditioned for conventional statistical inference. The purpose of this study was twofold. First it was to provide proofof-principle that ME-ICA and its related methodological developments can be applied to BOLD fMRI data from animal models. Moreover, this study focused primarily on the analysis of individual animals. This was done to establish a platform for studying disease and pharmacological models based on individual variation — disease severity or treatment dosage, for example. In doing so, we hope to lay the groundwork for analysis of BOLD signals in other model systems and varying neurobiological conditions. Second, this study sought to optimize the experimental parameters towards detection of spontaneous BOLD signals in individual rats using multi-echo fMRI. This involved testing multiple pulse sequence configurations to find one that optimized BOLD sensitivity versus imaging resolution. Our prior studies indicated that an effective measure of general BOLD sensitivity was BOLD DOF. Across cohorts of 3 T human data, we previously found that greater subject motion decreased BOLD DOF, because spin-history artifacts and spin-density fluctuations decrease sensitivity to subtler BOLD fluctuations (Kundu et al., 2013). Assuming other experimental parameters are fixed, the overall BOLD sensitivity of varying pulse sequences would be similarly captured in BOLD DOF. We explored BOLD DOF in an experiment where anesthesia dosage was varied parametrically. This measure was expected to be particularly helpful for the analysis of spontaneous BOLD, since ordinary model-specific sensitivity measures such as contrast-to-noise ratio are not applicable. Optimizing the imaging protocol also involved determining an anesthesia protocol that would not suppress amplitude of BOLD activity to levels below detectability (Kannurpatti et al., 2008). Thus, after determining an optimal pulse sequence, a separate experiment was conducted to determine an isoflurane anesthesia dosage that sufficiently suppressed animal movement without over-suppressing BOLD functional variability. In three rodents, inhaled isoflurane anesthesia was varied stepwise (1.0, 1.5, 2.0%) and BOLD DOF was then counted via the BOLD vs. non-BOLD independent component selection procedure. The variation of BOLD DOF versus anesthesia dosage helped determine if there was a simple relationship between the two, or if there was particularly high BOLD signal complexity in rodents when using low anesthesia dosage.

Methods Animals Experiments were conducted on 7 healthy male Sprague–Dawley rats of 250–350 g body weight. MRI images were acquired in vivo.

P. Kundu et al. / NeuroImage 102 (2014) 861–874

This study was approved by the local ethics board of the Brain and Spine Institute, Pitié Salpêtrière Hospital (Paris, France). Acquisition MRI data were acquired using a Bruker BioSpec 11.7 T small animal MRI scanner (Bruker BioSpin, Ettlingen, Germany) equipped with a 72 mm-diameter resonator used for emission and a 4-channel phased-array surface coil used for reception. The scanner was interfaced to a console running Paravision 5.1. Anatomical images were acquired using a 2D multi-slice spin-echo sequence with parameters: TR = 6 s, TE = 15 ms, spectral width 50 kHz, matrix size 128 × 128 × 100 yielding 100 coronal slices at 250 μm isotropic voxel size, requiring 12 min 48 s for acquisition. For multi-echo functional scanning, 3 fMRI pulse sequences using 2-D echo-planar imaging (EPI) acquisition were piloted, detailed in Table 1. The sequences varied by voxel size and the number of volumes acquired; all sequences were used to acquire whole-brain data, with flip angle (FA) equal to Ernst angle. Anesthesia During MRI scanning, animals were anesthetized and their heads were immobilized with a tooth bar and ear bars. Respiration was monitored for the duration of the experiment to ensure animal stability. Body temperature was monitored throughout the experiment based on rectal temperature, and was maintained at 37 °C with a thermoregulating circulating water-filled heating system integrated to the cradle. For experiments to determine a preferred pulse sequence based on BOLD acquisition sensitivity as reflected by BOLD DOF, medetomidine hydrochloride (Domitor®, Pfizer) was infused in the tail vein (200 μg/kg/h) using a subcutaneous bolus injection of 0.45 mg/kg as described by Pawela et al. (Pawela et al., 2009). Medetomidine was administered in parallel with isoflurane for induction at 5% mixed in oxygen and air at 1:5 (flow of 1.5 L/min) and maintenance (1.5%). After completing these parallel-anesthesia experiments, anesthesia was reversed with atipamezole hydrochloride (Antisedan®, Pfizer) at 0.3 mg/kg. These experiments were conducted on 4 rats; each of the three sequences was piloted on a different rat. The sequence that was determined to yield maximal BOLD DOF was used to scan the fourth rat. The preferred pulse sequence was then used to determine BOLD variability in response to varying isoflurane anesthesia dosage. In this experiment, only isoflurane anesthesia was used in order to allow maximal expression of BOLD variation (Williams et al., 2010). “Global” BOLD variability was expressed by BOLD DOF when BOLD acquisition sensitivity was held constant. After isoflurane anesthesia induction (5% isoflurane mixed in oxygen and air at 1:5, flow of 1.5 L/min), three different maintenance anesthesia levels were tested in order: 1.0%, 1.5%, and 2.0%. After each fMRI scan, the dosage of inhaled isoflurane anesthesia was increased to the next level, and the next fMRI acquisition followed. Image processing

863

used for further volumetric preprocessing (Cox, 2012) (Fig. 1). Each NIFTI 3D + t dataset was first corrected for slice timing offsets (3dTshift). Next, the 3D + t dataset of the first TE (i.e. with highest image signal-to-noise ratio) was used to determine volume-to-volume motion correction relative to a middle volume, computed as parameters of rigid-body realignment using rotations and translations (3dvolreg). The motion correction parameters computed from the first echo dataset were applied to the 3D + t datasets of the remaining echoes (3dAllineate). A brain mask was computed from the motion correction base volume using a spherical inflation model (3dSkullStrip). To do this effectively, image regions below the inferior edge of the brain were deleted, the brain was “shrunk” in the rostral–caudal direction using a voxel-dimension resize to create a more spherical image (3drefit), and the -pushout routine of the inflation procedure was enabled. After masking, the shrunk direction was re-expanded by restoring the original voxel dimension size. Separately, the high-resolution anatomical image was brain-extracted with a standard application of the inflation approach (3dSkullStrip). This masked anatomical was aligned to the functional using an affine (12-parameter) coregistration routine (3dAllineate). A Gaussian FWHM 0.5 mm smoothing kernel was applied to each 3D + t dataset within the brain mask to attenuate thermal noise by essentially averaging local time series, while also producing a masked 3D + t dataset (3dBlurInMask). Lastly, smoothed and masked 3D + t images of different TEs were concatenated in the zdirection as a convenient way to represent 3D + t + echo data without using 5-D arrays. This dataset was input to the ME-ICA decomposition and denoising procedures implemented in tedana.py, which unfolded the z-concatenated data into the 5-D form for further analysis. Component level TE-dependence analysis Component-level TE-dependence was utilized in this study. This analysis implemented ME-PCA dimensionality reduction, ME-ICA decomposition, and denoising. The procedure is summarized here in the list form, and a more detailed description is available in Kundu et al. (2013). The following steps were executed in order: 1 Computation of T2* and S0 parameter maps from fits of voxel-wise time course means from different TEs to the log-linear transformation of mono-exponential decay (Eq. (1))
  SðTEn Þ ¼ S0 exp −R2 TEn

ð1Þ

where n is the echo number, TEn is the nth echo time, S(TEn) is the time series mean at echo n, S0 is the initial signal intensity, and R2* is the inverse of relaxation time or 1/T2*. S0 and R2* are solved for, by fitting the signal means from 3 or more TEs. 2 Weighted averaging of time series of different TEs to produce “optimally combined” data, where the time series at each TE was weighted (Eq. (2)) by    ω T 2;v ¼ n

  TEn  exp −TEn =T 2;v N X

   TEn  exp −TEn =T 2;v

ð2Þ

n

Python routines (using Numpy, Nibabel) were used to format scanner images into NIFTI image time series datasets of volume by time (3D + t), one per TE (Taylor and Perez, 2011). AFNI tools were

where v is a voxel, N is the total number of echoes, n is a particular echo, and ω is the weight for each echo number n.

Table 1 Table of multi-echo fMRI pulse sequences piloted on the Bruker BioSpec 11.7 T small animal MRI scanner. Sequence

Geometry

Timing

1 2 3

0.3 × 0.3 × 1.2 mm voxel size; 16 coronal slices w/0.2 mm gap; FA = 77 (Ernst) 0.4 × 0.4 × 0.6 mm voxel size; 32 coronal slices w/0.2 mm gap; FA = 77 (Ernst) 0.5 × 0.5 × 0.5 mm voxel size; 40 coronal slices w/0.2 mm gap; FA = 77 (Ernst)

TR = 3 s; 3 TEs = 8.5, 27.0, 45.5 ms; 110 reps (330 s duration) TR = 3 s; 3 TEs = 7.6, 19.6, 31.6 ms; 500 reps (1500 s duration) TR = 3 s; 3 TEs = 8.25, 20.25, 32.25 ms; 400 reps (1200 s duration)

864

P. Kundu et al. / NeuroImage 102 (2014) 861–874

Despike/ Volume/ Anatomical Co-register

ICA Mixing Matrix

Component Maps & Time Series

MEICA

Separated noise and BOLD time series

T2* Weighted Combination

Raw

Least squares fit ME-ICA Automatic Processing

Results

Fig. 1. Diagram illustrating essential steps in ME-ICA procedure. Raw data is acquired using a multi-echo fMRI sequence, and signals are composed into separate volumetric time series datasets, one for each signal echo acquired. EPI volumes are temporally and spatial aligned, co-registered to anatomy, and smoothed to attenuate thermal noise (0.5 mm FWHM Gaussian kernel). (c) Multi-echo data are “optimally” averaged with T2* weights to further attenuate thermal noise and approximate acquisition at optimal contrast (TE = T2*) at each voxel. Multiecho ICA, including multi-echo PCA dimensionality reduction and FastICA unmixing computes spatial ICA. (e) A time domain mixing matrix is generated and is fit to the full-rank optimally combined functional data, producing component maps. Multi-echo component maps are generated for functional connectivity mapping and TE-dependence analysis. BOLD components are concatenated to produce a BOLD component coefficient dataset, and non-BOLD components are projected out of optimally combined time series to produce denoised time series.

3 Computing a data mask that excluded low intensity signals. 4 Variance normalization of the time series within the intensity mask, producing a voxel-by-time series data matrix from the mask, then mean-centering and variance normalizing the data matrix. 5 PCA of time series, implemented as singular value decomposition (SVD, Eq. (3)) with partial matrices X ¼ USV

T

ð3Þ

where X is the variance-normalized data, U and V are the left and right singular vectors, and S is the vector of singular values. 6 Principal component signal amplitudes were computed by multiple least squares fit of the PCA time courses (V, the right singular matrix) to the preprocessed signal-unit time series of each TE. 7 For each component at each voxel, the principal component signal amplitudes at different TEs were fit to linear TE-dependence and TE-independence models, and corresponding F-statistics for goodness of fit were computed. The TE-dependence and TE-independence models are (Eq. (4a) and Eq. (4b)) 

ΔSTE =STE ¼ −ΔR2 TE

ð4aÞ

ΔSTE =STE ¼ ΔS0 =S0

ð4bÞ

where ΔSTE is the signal change from mean for a fluctuation at a TE (i.e. its β weight from least squares fit), and STE is the signal mean at a TE, ΔR2* is the change in susceptibility-weighted transverse relaxation time that is solved for in the TE-dependence BOLD model, S0 is the initial signal intensity, and ΔS0 is the change in initial signal intensity that is solved for in the TE-independence non-BOLD model. The respective F-statistics are then determined (Eq. (5a) and Eq. (5b)): α0 −αR

2

F R ¼ 2

αR

2

d: f:0 −d:f:R

ð5aÞ

2

d:f:R

2

α0 −αS0 F S0 ¼

αS0 d:f:0 −d:f:S0 d:f:S0

ð5bÞ

where α0 is the null variance (∑β2c,v,TE), αR is the variance explained 2 by the fit to the TE-dependence model, αS0 is the variance explained by the fit to the TE-independence model, d.f.0 is the total number of degrees of freedom for TE-dependence models (equal to number of echoes, 3), and d:f:R and d:f:S0 are the degrees of freedom used in 2 the respective fits (1 each). κ and ρ statistics were computed for each principal component (Eq. (6a) and Eq. (6b)), V X

κc ¼

v

V X

ρc ¼

v

p

zc;v F c;v;R 2 X p zc;v

ð6aÞ

p

zc;v F c;v;S0 X p zc;v

ð6bÞ

where c is the component index, v is the voxel index, V is the total number of voxels, z is the normalized signal power, and p is a power factor (default 2). 9 Thresholds for κ, ρ, and eigenvalue were determined, after sorting values, as the inflection points of the respective Scree plots (i.e. using an elbow-finding routine). 10 Optimally combined time series were dimensionally reduced to the principal components with above-threshold κ, or ρ, or eigenvalue (i.e. isolating any correlated signal associated with MR contrast), followed by re-centering and re-normalizing the time series, then the full data matrix. 11 FastICA of dimensionally reduced data using the tanh contrast function produced time-domain independent component mixing matrix (variance normalized). 12 Least-squares fitting of the FastICA mixing matrix to separate-TE time series data, in order to compute per-voxel TE-dependence and TEindependence model fits, respective F-statistics, and componentlevel κ and ρ values that indicated BOLD and non-BOLD weighting, respectively. This involved the application of Eqs. (4a), (4b)–(6a), (6b) for ICA component β weights, in essence replicating steps 7–8 after the ICA step for the ICA time course mixing matrix. 13 Thresholds κthresh and ρthresh were computed as the inflection points of their respective Scree plots. BOLD components were identified as high-κ/low-ρ, i.e. κ N κ thresh and ρ b ρthresh, and non-BOLD components were identified as the reverse. 14 Individual BOLD ICA components were mapped by fitting the ICA mixing matrix to the preprocessed, optimally combined, and

P. Kundu et al. / NeuroImage 102 (2014) 861–874

variance normalized time series, in effect producing partial correlation coefficients, R, per-component, per-voxel. R values were then converted to Z values using the variance stabilizing Fisher R–Z transformation and performing spatial Z-normalization for each map, which allowed thresholding each component (Z N 3, uncorrected P b 0.01 used here). Thresholded component maps were then overlaid on co-registered anatomicals, upsampled with linear interpolation from 0.5 mm resolution to 0.2 mm resolution (only for display purposes), and colorized on a yellow-cyan color spectrum. Conventional ICA decomposition To compare ME-ICA decomposition to a conventional ICA approach that does not utilize multi-echo fMRI information, the component representation from ME-ICA for an individual animal was compared for the same dataset to the representation from a conventional ICA pipeline — probabilistic ICA as implemented in FSL MELODIC ICA (Beckmann and Smith, 2004). Probabilistic ICA can be applied to one of the echo time series datasets from a multi-echo run. A given echo time series fairly approximates the acquisition of a single-echo dataset of the same TE, since multi-echo EPI just involves repeated EPI readouts after a conventional gradient-echo excitation pulse. The preprocessed first echo dataset (i.e. after slice timing alignment, rigid body motion correction, and masking) was input to MELODIC with default parameters: tanh contrast function; automatic dimensionality estimation based on probabilistic PCA (Tipping and Bishop, 1999). The component maps were produced using the default 3-class mixture-model thresholding (Z N 3), and maps were ordered by variance explained (the default) and rendered similarly to ME-ICA component maps (above, step 14). Multi-echo independent coefficients regression (ME-ICR) For individual datasets, ME-ICR seed-based functional connectivity maps were computed using the high-κ/low-ρ (determined at step 13) BOLD component basis extracted from ME-ICA decomposition of a given dataset. After selecting a seed voxel a, the vector of BOLD independent coefficients from that voxel was extracted and used to compute Pearson correlation with BOLD independent coefficient vectors from all other voxels b (Eq. (7)). This produced a map of correlation values based on the computation (Eq. (7)): Ra;b ¼

Ia Ib  Ia Ib

ð7Þ

where Ra,b is the Pearson correlation of coefficient vectors Ia and Ib. Next, each correlation value was converted to a standard score using the standard Fisher R–Z transform, using the standard error term for normalization (Eq. (8)).   pffiffiffiffiffiffiffiffiffiffiffiffiffi Z a;b ¼ arctanh Ra;b  Nc −3:

ð8Þ

The DOF for standard error was the number of BOLD independent components (Nc). The significance (p) value for this standard score (Z) was computed from standard normal cumulative distribution function, and statistically significant Z-values (P b 0.05) were anatomically mapped. Results Multi-echo images and optimal combination Multi-echo images from fMRI acquisition from a representative animal are shown (Fig. 2, top). From the 3-echo sequence #3, a corresponding number of images are produced from each TR. With increasing TE, image intensity decreases. In this sample, images at TE = 7.6, 19.6, 31.6 ms show anatomical features such as veins (dark channels), cortex,

865

and subcortex. Signal dropout near the temporal bone (white arrow) and ventral regions increases with TE, indicating particularly short T2* in these regions (optimal contrast is at TE = T2*, dropout occurs when TE N N T2*). Maps of T2* and S0 were produced from the log-linear fit of signal means (Fig. 2, bottom). The S0 map reveals high initial signal intensity in cortex and low signal intensity in draining veins. The T2* map shows cortical T2* (at the present EPI acquisition parameters) on the order of 10 ms, and subcortical T 2 * on the order of 20 ms. The image resulting from T2*-weighted optimal combination of the three images produced an image that approximated acquisition with TE = T2* at each voxel. In this synthetic image, the brightness of the cortex in the image at TE = 7.5 ms is reduced while the contrast of the subcortex in the image at TE = 19.6 ms is preserved, altogether producing an image with optimal functional contrast over the whole brain. Linear TE-dependence analysis of component signal changes FastICA decomposition of optimally combined data produced approximately statistically independent spatial maps for each component, alongside a corresponding time-domain mixing matrix. The component maps directly from FastICA were in units of partial correlation. Fitting the mixing matrix to the separate signal-unit datasets of each TE yielded signal change values for each component at each TE. The variation of signal change magnitude with TE was assessed for each component to compute maps with the following metrics at each voxel: ΔT2* (change in susceptibility-weighted transverse relaxation due to BOLD contrast modulation), ΔS0% (percent change in spin-density due to artifact), FR2*, and FS0 (goodness of fit to signal models of ΔT2* and ΔS0). Counting the signal change maps at separate TEs and optimal combination, a 3echo acquisition produces 8 unique maps in addition to the 1 gained from FastICA alone, usable as aids in component assessment and classification into BOLD or non-BOLD categories without referencing spatial localization or time courses. For a given component (Fig. 3, representative of an insula network), the average of FR2* and FS0 values, weighted by Z-normalized partial correlation values, produced the metrics κ and ρ for component-level TEdependence and TE-independence, respectively. The percent signal change maps for the component varying over TE show intensification of signal with TE in areas of peak signal. Plotting the percent signal changes with TE at a representative voxel reveals a linear trend. Fitting the linear TE-dependence model to these signal changes yields FR2* = 198, which corresponds to R2 = 0.99. κ and ρ signify component-level BOLD and non-BOLD signal modulation with TE, respectively. Functional networks express high κ and low ρ, indicating high BOLD weighting and low non-BOLD (artifact) weighting. The weighted average of FR2* over all voxels is κ = 54, and that of FS0 is ρ = 15.3. Reference points for significance of κ for 3-echo data at P = 0.05 → F = 18.5 and P = 0.025 → F = 38.5 further indicate that the present component is BOLD-weighted. The 4-fold difference between κ and ρ along with failure of ρ to pass these reference values as well as the respective threshold values suggest that the present component is a functional network and not physiological artifact. Whole multi-echo dataset decomposition and component classification The decomposition of preprocessed multi-echo data into ICs was done in a two-step process involving PCA dimensionality reduction followed by ICA decomposition, which is the standard pipeline for most ICA applications. The κ, ρ, and variance (%) spectra produced by the PCA decomposition typically demonstrated an inflection point that could be used for component selection thresholds (Fig. 4, left). Any component with any metric exceeding its respective threshold was included for ICA. Components with no such metric were treated as thermal noise. Examining the spectra revealed that all PCA components explaining significant variance also had high values for κ and ρ. Additionally, most

866

P. Kundu et al. / NeuroImage 102 (2014) 861–874

Sagittal views of raw multi-echo data and parameter maps TE=7.5ms TE=19.6ms TE=31.6ms D

(a)

Rostrocaudal axis

Ctx

Cb

V

T2* Map

S0 Map D

(b) Cb

T2* Weighted Image

Rostrocaudal axis

Ctx

V 50ms Fig. 2. (a) Raw multi-echo images from 2-D EPI sequence #3 showing signal decay over TE = 7.5,19.6,31.6 ms. Sagittal images are shown, labeled to indicate (D)orsal and (V)entral directions. Middle and last TEs show substantially less signal amplitude than first TE, indicating a fast decay due to short T2* at 11.7 T (b) T2*and S0 maps are computed by fitting signal means across TEs (over time) to a monoexpoential decay, voxelwise. Cortical T2* is 10 ms, subcortical T2* is 20 ms. The T2* weighted image has optimal image contrast and reduces thermal noise in time series.

components with a high (supra-threshold) κ value also had a high ρ value. Not all components with high κ or ρ had high variance. This pattern concurred with the variance and TE-dependence analysis findings from human subject data. It suggested that PCA could not separate TEdependent and TE-independent signals on the sole basis of covariance, but it does pool robust MR signals into high variance components. Furthermore, this analysis showed that significant TE-dependent BOLD signals could be captured in principal components but could go undetected due to low variance. The multifaceted component selection from using κ, ρ, and variance thus increased the sensitivity of dimensionality reduction. Dimensionally reduced and optimally combined (T2* weighted) data were decomposed using FastICA using the tanh contrast function. The resulting set of components were analyzed for TE-dependence and TEindependence to produce κ, ρ, and variance statistics. The ranked series of κ values, the κ spectrum, alongside ρ and variance spectra were first used to evaluate decomposition performance (Fig. 4, right). The κ spectrum shows segregation into high and low κ regimes, signified by the characteristic set of high values followed by a long tail, separated by a clear inflection point (elbow). Referencing the respective ρ values shows that components with highest κ values have low ρ values. On the other hand, a number of low κ components have very high ρ values, which in some cases exceeded highest κ values. This shows that while PCA failed to separate TE-dependent and TE-independent signals based on their source distributions, FastICA was unequivocally able to do so. As for human data, this clear determination was possible without invoking network templates, time course spectral signatures, or other assumptions therein. Finally, examining variance spectra shows a majority of low variance components (b1%), with the exception of extreme variance components (10–40%). These extreme variance components coincide with the inflection point, suggesting low TE-dependence. This pattern is suggestive of a pervasive physiological artifact, scanner drift, or other low-rank phenomena.

BOLD component evaluation The maps and time courses of high κ components were assessed for spatial localization of signals and time course spectral properties. Three examples are shown (Fig. 5). The component with highest κ is shown to have specific cortical localization to the primary motor cortex (top, ventral view shown). The time course and its frequency spectrum (FFT) shows a single dominant ultra-low frequency band, with the remainder of the spectrum showing approximately equal amplitude frequency bands throughout. Referencing the time course, the ultra-low frequency band corresponds to a negatively trending signal attenuation during the first half of the acquisition duration. While the quality of spatial localization (i.e. image contrast) for this component agrees with that in human resting fMRI data, the time course spectral properties do not correspond to wakeful human data. The conspicuous “low frequency” 0.01–0.1 Hz range of activity is apparent but heavily attenuated compared to human data, in which it explains high amounts of variance (75% or more, normalizing for bandwidth). The component with next highest κ value spatially localizes to the primary sensory cortex, lateral and adjacent to the primary motor cortex component (middle). Its map shows localization quality equivalent to the primary motor cortex component. Its time course shows a negative attenuation similar to that of the motor cortex, which is reflected in its frequency spectrum. In contrast to the motor cortex component, the time course of the sensory cortex component firstly has a fully attenuated low frequency signature, but then it also has a relatively enhanced “high frequency”, i.e. N 0.10 Hz signature. The present signal cannot be thermal noise since it has strong spatial coherence as evidenced by the corresponding component map. However, this component may have greater inclusion of TE-independent (i.e. spin-density) signals due to the higher ρ value. Finally, the component with lowest κ of the set has a map showing spatial localization to the caudate putamen (bottom). The low κ indicates relatively weak TE-dependence, due possibly to

P. Kundu et al. / NeuroImage 102 (2014) 861–874

867

(a) Axial % Signal change maps with TE Rostrocaudal axis

L

R

TE=7.5ms

TE=19.6ms

(b) % Signal change with TE

TE=31.6ms

(c) TE-dependence (F) map

% Signal Change

Rostrocaudal axis

12

L

9

F=198

6

R = 0.9882

3 0

0

10

20

30

R

40

=54.0 =15.3

TE(ms)

Fig. 3. (a) ICA component with strong TE-dependent scaling of percent signal changes in area of high signal change (blue box, above), shown in axial slice labeled to indicate (L)eft and (R)ight directions. (b) Signal changes at indicated voxel versus TE are plotted and fit to linear TE-dependence model (lower left). (c) Based on component images of each TE, an F-statistic fitparameter map (lower right) shows areas of high linear scaling. The metrics κ and ρ are averages of F values over all voxels (weighted by normalized signal intensity), for TE-dependence and TE-independence, respectively.

(b) ME-ICA Decomposition

Variance (%)

Variance (%)

(a) ME-PCA Dimension Est.

Principal Components (by variance) ME-ICA Dimensionality Estimate + additional BOLD weighted comps Probabilistic Dimensionality Estimate

Independent Components (by ) BOLD components are low variance

Fig. 4. (a) ME-ICA uses κ and ρ for PCA dimensionality reduction, called ME-PCA. After κ, ρ, and variance values are sorted and assessed for a threshold (e.g. “elbow” of ME-PCA variance), all components with at least one parameter passing a respective threshold are included for ICA decomposition. ME-PCA allows low-variance PCA components with detectable BOLD or non-BOLD signal to be used in ICA decomposition. For rat rs-fMRI, variance based component selection (blue line) discards BOLD signals (orange arrows). (b) After ICA decomposition, κ and ρ are computed for each ICA component to detect BOLD ICA components. Notably, ICA components with substantial BOLD or non-BOLD weighting can have relatively low variance compared to artifacts.

868

P. Kundu et al. / NeuroImage 102 (2014) 861–874

Axial Slice

Signal Time Course

Fourier Transform

L =102.8 =12.6

Frequency Amplitude

=74.5 =22.5

A

Signal

R/P

=54.0 =15.3

Time (sec)

Frequency (Hz)

Fig. 5. Three BOLD weighted components from rat #2 identified by high κ and low ρ. Maps (left) show localization to caudate-putamen, primary motor, and primary sensorimotor networks. Component time courses (middle) show initial drift (8 min.) followed by stabilization. Frequency spectra (right) do not indicate predominantly low frequency (b0.10 Hz) signal. Primary somatosensory cortex shows predominantly higher frequency (N0.10 Hz) fluctuations.

less robust BOLD activity. Its localization quality is equivalent to the other two components, despite being a subcortical component. Also like the former two, its time course shows a low frequency attenuation, which is reflected in its frequency spectrum. The remainder of the spectrum shows no conspicuously high amplitude ranges of frequencies. Yet, the signal is sufficiently correlated to produce the quality of localization in the map, and the ρ value is too low to suggest a pronounced spin-density artifact. Non-BOLD component properties The maps and time courses of components with low κ and/or high ρ values were assessed for spatial and temporal properties consistent with artifact. The components with highest ρ values had spatial patterns resembling imaging artifacts (Fig. 6, top). Other such components had spatial localization along vasculature, with time courses showing periodic high frequency signal (bottom). Finally, a minority of low κ components with localization and time course properties that were by-andlarge artifactual showed anomalous patterns in a few slices which if viewed alone could make component classification ambiguous. These patterns indicated, for example, signal localized along ventricular walls or symmetric patterns localized near cortex (middle). However, TE-dependence analysis clearly indicated that these components were not functionally related BOLD components of the same regime as high-κ components. Determining an optimal pulse sequence The three pilot sequences varied in terms of voxel volume and scan length. For imaging of the anesthesia state induced with medetomidine, all sequences produced a comparable number of BOLD components. Sequence 1, a five-minute scan with voxel volume (0.108 mm3, 0.3 mm × 0.3 mm in-plane resolution), produced 5 high-κ components including ones representing dorsal and ventral

striatum. Sequence 2, which was acquired for a much longer period (25 min) and scanned at higher in-plane resolution (0.4 mm) but lowest overall volume (0.096 mm3) produced 6 high-κ components covering the cortex, but did not yield subcortical components. Sequence 3, which was acquired for 25% shorter period (20 min) but with isotropic voxel size (0.5 mm) giving largest voxel volume (0.125 mm3) also produced 6 high-κ components, including one in the subcortex. Using the same sequence in a fourth rat produced the same number of components, also with one in the basal ganglia. Assuming the small differences in voxel volume did not fundamentally change BOLD sensitivity, the lack of proportional increase of BOLD DOF with scan length suggested that the modes of BOLD variability are relatively limited at the present anesthesia dosage. Longer sequences were still favored, however, under the consideration that lower anesthesia doses could yield greater network variability over time. To have versatility in functional mapping over various imaging planes, sequences closer to isotropic voxel size were also favored. These considerations alongside relative reliability observed from repeated acquisition led to the election of the 0.5 mm isometric sequence as a preferred sequence. The networks found using this sequence and the medetomidine-isoflurane anesthesia protocol are shown (Fig. 7). Variation of BOLD DOF with anesthesia type and dosage BOLD DOF was assessed as a function of varying isoflurane anesthesia level. Keeping in mind the application of the present experimental protocol towards future pharmacological models, a high BOLD DOF is generally favorable because it indicates high sensitivity to more complex fluctuations. On the other hand, retaining a stable BOLD DOF with varying anesthesia is favorable for comparing functional connectivity and its response to pharmacological or surgical treatments across different anesthesia levels. Very low BOLD DOF is unfavorable since statistical inference for functional connectivity becomes untenable, and major brain regions will not be significantly or consistently represented in

P. Kundu et al. / NeuroImage 102 (2014) 861–874

Coronal Slice

Signal Time Course

869

Fourier Transform

Frequency Amplitude

=16.7 =107.4

Signal

D =24.8 =32.8 V/L

R

=24.7 =22.9

Time (sec)

Frequency (Hz)

Fig. 6. Artifact components identified by low κ and/or high ρ (thresholds based on κ-spectrum elbows) with ambiguous spatial localization and/or time course spectrum. Component maps (left) show coronal slices. Time series (middle) and frequency spectra (right) are plotted, showing no distinguishing time course of frequency characteristic. Component with the smoothest image (bottom), most similar to a functional component, has the strongest spin-density weighting, definitely indicating artifact.

connectivity maps. Previous studies have already indicated that BOLD signal amplitude in rodents is attenuated with increased anesthesia levels, with decreased neuronal activity being one model proposed of the isoflurane mechanism (Kannurpatti et al., 2008). We found that BOLD DOF decreases with increasing isoflurane anesthesia levels; in the three animals that were tested, BOLD DOF consistently decreased across 1.0%, 1.5%, and 2.0% dosage levels (Table 2b). The mean and standard deviation of BOLD DOF at the 1.0% dosage was 17.3 ± 6.8. Most interestingly, rat 1 at the 1% anesthesia level produced 25 BOLD components from a single 20 minute scan. These components are mapped and described in the following section. With the increase to the 1.5% dosage, mean BOLD DOF was reduced to 10.3 ± 5.8. By the 2.0% dosage, it was reduced to 6.0 ± 2.6, meaning the BOLD DOF is reduced by a factor 2–3 in response to doubling the dose of anesthesia. Notably, BOLD DOF at the 2.0% inhaled isoflurane dosage was consistent with that for the combined medetomidine-isoflurane treatment. These findings suggest careful consideration of anesthesia choices in spontaneous BOLD experiments, as well as careful interpretation of BOLD variability in response to anesthesia. BOLD component mapping in individual rodents Twenty-four BOLD independent components derived from an individual rat anesthetized with 1.0% inhaled isoflurane were mapped (Fig. 8). These components are representative of an individual animal, derived from a single ME dataset, from one ME-fMRI run. This number of BOLD DOF reflected a complexity of spontaneous BOLD signals on par with that observed in human data acquired at 3 T. Inspecting spatial maps of components revealed clear localization to brain regions approaching the anatomical scale described by the Paxinos atlas (Paxinos and Watson, 2006), with several components indicating anatomical delineations of functional areas consistent with that seen in canonical functional neuroanatomy. Both subcortical and cortical areas were identified, with some cortical areas showing segmentation into lateralized regions. Lateralized components can be observed in highdimensional sparse decomposition of fMRI data from individual human subjects, but to the best of our knowledge this is the first such

observation for rodent. Numerous bilateral cortical areas were identified, including the primary and secondary cingulate areas, primary visual areas, a number of sensory cortex areas in separate components, among others. A number of largely unilateral components were also observed, notably including left and right motor cortex, possibly in more than one pair of components. Subcortical components represented caudate-putamen, hippocampus, and nucleus accumbens. Comparison of ME-ICA decomposition to a conventional ICA decomposition Component maps obtained from ME-ICA were compared to the maps obtained from probabilistic ICA as implemented in FSL MELODIC. This comparison was made for the dataset that ME-ICA identified to have greatest BOLD complexity (25 components), as that dataset would likewise give MELODIC the best chance of extracting networklike maps. MELODIC was applied to the earliest TE preprocessed time series from the multi-echo dataset (approximating a single-echo dataset acquired at TE ~ 8.5 ms). A total of 60 components were identified by MELODIC (as determined by probabilistic PCA) compared to 186 identified by ME-ICA (determined by ME-PCA). The top components from MELODIC as ranked by variance explained are shown (Fig. 9). In contrast to highest κ components derived from ME-ICA, the highest variance-ranked MELODIC components show substantial artifactual contributions, indicated by spatial patterns such as alternating positive–negative bands and un-clustered voxels. In general, MELODIC components do not show the level of precision in anatomical localization that is provided by ME-ICA component maps, lacking features such as linear delineation of cytoarchitectonic boundaries. Only two networklike maps were identified from the MELODIC component set, and only one of these were found in the top-ranked group, most comparable to area M2 as shown in the ME-ICA component maps. Seed-based functional connectivity in an individual rodent using ME-ICR The stacking of ME-ICA component maps from an individual scan produces a coefficient vector dataset representing an individual rat's spontaneous BOLD activity that can be used to compute seed-based

870

P. Kundu et al. / NeuroImage 102 (2014) 861–874

Prelimbic Cortex

Secondary Somatosensory Cortex

Primary Somatosensory Cortex

Piriform Cortex

Retrosplenial Granular Cortex

Cingulate Cortex

Secondary Somatosensory Cortex

Primary Somatosensory Cortex

Caudate-Putamen

Rat #4

Rat #3

Retrosplenial Granular Cortex

Fig. 7. Consistent high κ (BOLD) components identified across different rats (#3 and #4) using sequence #3 at 1.2% inhaled isoflurane/oxygen anesthesia. Each component is shown in (top to bottom) coronal, axial, and sagittal views. Note that underlay images represent different anatomical image MRI acquisition parameters. Retrosplenial (cingulate), primary and secondary somatosensory cortex components align well across rats. Notably, retrosplenial cortex corresponds to the “default mode” configuration. Piriform cortex (amygdala/insula and caudate– putamen) manifests more variably across animals.

functional connectivity using ME-ICR. In human imaging, ME-ICR yields functional connectivity maps with high specificity that can be thresholded according to conventional techniques of statistical inference. While components map out the spatial basis of functional connectivity over the entire brain, ME-ICR seed connectivity maps probe the functional connectivity of a region of interest to other regions, which emerges from the overlap or intersection of basic connectivity processes involving the seed region. Seed-based connectivity is computed in standard Z-units, such that connectivity values surviving a threshold

Table 2 BOLD spatiotemporal degrees of freedom (DOF) for three dosage levels of isoflurane anesthesia repeated across three rats. Isoflurane anesthesia vs. BOLD sensitivity (as DOF) Sequence: 3 Rat

1 2 3 Mean Std

of Z ≥ 1.95 indicates significance of P ≤ 0.05. For the component set shown in Fig. 8, example ME-ICR maps are shown in Fig. 10 from mapping functional connectivity within an individual rat using four seeds at the P ≤ 0.05 significance level. Overall, maps reflect the high spatial specificity of seed-based functional connectivity that was originally observed from human ME-ICR analyses. The right secondary somatosensory motor cortex (S2) shows connectivity to the ipsilateral primary somatosensory motor cortex (S1) as well as the contralateral secondary somatosensory cortex. This pattern is reflected when placing a seed in the left S2. Seeding in left globus pallidus (GP) indicates functional connectivity to bilateral primary and secondary motor cortex (M1M2). Finally, seeding in the left ventral pallidum (VP) indicates connectivity to bilateral caudate-putamen (CPu). Discussion

Isoflurane anesthesia dosage 1.0%

1.5%

2.0%

25 15 12 17.3 6.8

17 6 8 10.3 5.8

9 4 5 6 2.6

In this study we applied ME-ICA to fMRI data from anesthetized rodents that were scanned in an 11.7 T small animal MRI. This study helped determine the efficacy of applying ME-ICA to data acquired from animals, as well as data acquired at ultra high field. Both conditions differ from the original development platform of ME-ICA, which was fMRI data from wakefully resting human subjects scanned at 3 T (Kundu et al., 2011). We hypothesized that analysis of multi-echo

P. Kundu et al. / NeuroImage 102 (2014) 861–874

871

M1M2

V1

RSGc

S1

S1

S1

DCPu

DCPu

S1

VCPu

M2

S1

S1

Au

S2

S1

S1

M1M2

S1BF

S1

S2

Ect

S1

Cg2

Fig. 8. Spatial localization maps for 24 independent components of spontaneous BOLD activity, ranked by BOLD weighting (κ), of an individual rat anesthetized with 1.0% inhaled isoflurane. Each component is shown in (top to bottom) coronal, axial, and sagittal views. Labeled maps include: left and right primary and secondary motor cortices (M1/M2), primary visual cortex (V1), retrosplenial granular cortex (RSGc), primary and secondary sensory cortices (S1/S2), dorsal and ventral caudate-putamen (DCPu/VCPu), primary and secondary auditory cortex (Aud), cingulate cortex (Cg2).

signals according to physical criteria of TE-dependence and TEindependence, combined with signal decomposition according to the statistical criterion of “independence” or sparsity was sufficiently general to first elucidate and then classify BOLD and non-BOLD component in these novel datasets. We succeeded in demonstrating the stability of high-dimensional ME-ICA decomposition in explaining 90–95% dataset variance. Notably, even for datasets with high DOF, BOLD signals comprised on the order of 15% of explained variance, with the remainder

being by non-BOLD artifact. ME-ICA κ and ρ spectra derived from rat data were then found to follow the patterns of spectra obtained from human data, thus allowing automatic component classification by similar means. BOLD functional network components elucidated by the ME-ICA selection procedure were found to be similar across rats at the same anesthesia dosage; and at lower anesthesia dosage, observed functional neuroanatomy increasingly resembled canonical cytoarchitectonic descriptions of rat neuroanatomy. The use of

872

P. Kundu et al. / NeuroImage 102 (2014) 861–874

M2

Fig. 9. Spatial localization maps for 12 independent components derived from a conventional ICA pipeline that does not take advantage of multi-echo information, ranked by variance explained (from MELODIC probabilistic ICA, which does not utilize multi-echo optimal combination or MR signal specific ME-PCA high-dimensionality estimation). Most components show ambiguous or artifact patterns, with only one component of the first 12 clearly showing a functionally-related pattern — the secondary motor cortex (M2).

component maps in the seed-based functional connectivity mapping procedure known as ME-ICR was furthermore demonstrated. Given these results, we determined that ME-ICA is a valid and effective means of studying spontaneous BOLD fluctuations in anesthetized rat, and may have further applications in other field strengths and other animal models. BOLD sensitivity was measured as the DOF of BOLD spatiotemporal fluctuations, counted as the number of high-κ/low-ρ ICA components. BOLD DOF estimates the level of structured information content in the BOLD subspace of an fMRI dataset, and is closely related to concepts of rank, compressibility, or complexity. Higher BOLD

S1

S1

DOF indicated that BOLD fluctuations were less compressible and had greater information content. Lower BOLD DOF meant that the BOLD subspace had a smaller proportion of unique information. BOLD DOF, as a measure of structured information content, is an effective indicator of BOLD sensitivity, expected to decrease as thermal (Gaussian) noise dominates a dataset, and increase as the sampled signal increases in complexity. The pulse sequence configurations used in our tests of BOLD sensitivity acquired signals in either elongated or isotropic voxels, but all with comparable voxel volume (approximately 0.1 mm3). The fMRI sequences also varied in acquisition time from 5 to over 20 min. All

M1M2

M1M2

S2 S2

S2

CPu GP

CPu VP

Fig. 10. Seed-based functional connectivity using ME-ICR based on component coefficient data of the ME-ICA decomposition represented in Fig. 8 for four seeds from left to right, marked with green cross-hairs. Seed regions are marked with green crosshairs. Regions with significant (P ≤ 0.05) surviving functional connectivity are mapped on a yellow-red color scale. Indicated regions are: primary and secondary somatosensory cortices (S1 & S2), primary and secondary motor cortex (M1M2), globus pallidus (GP), medial forebrain bundle (mfb), and caudate-putamen (CPu).

P. Kundu et al. / NeuroImage 102 (2014) 861–874

sequences piloted in this study had broadly comparable BOLD sensitivity in terms of DOF, so the sequence with isotropic voxel size was preferred. Further work is required to test sequences more parametrically, possibly using different imaging trajectories (3-D EPI, for example) or using phase images to produce susceptibility variation maps (fQSM). Varying acquisition times (i.e. increasing data volume) nominally increased the sensitivity of decomposition to additional BOLD covariance patterns, but not in proportion to the extra acquisition time required. This finding suggests that, depending on anesthesia protocol, shorter fMRI scans can efficiently capture spontaneous BOLD signal variance in rat, producing component maps that represent cortical and subcortical areas. Further study of sequence optimization in the context of low anesthesia dosage could be informative. Somatosensory and motor regions appeared most consistently, but frontal cortex and piriform regions were found as well. Importantly, inspecting these BOLD component time courses shows equal proportions of b0.1 Hz oscillations and N0.1 Hz oscillations, which differs from the predominantly N0.1 Hz BOLD oscillations observed in human data both after ME-ICA denoising and conventional analysis without filtering (Fransson, 2005). Without κ and ρ characterization, these time courses could be ascribed to noise-related processes. Once the BOLD origin of these components had been established, however, spatial component maps indicated that these higher frequency fluctuations were highly localized to distinct cortical and subcortical areas. The absence of low-frequency dominance suggests that these BOLD components may be less influenced by very low frequency phenomena such as vasomotion. Further studies with larger sample sizes will be able to provide more quantitative analysis of BOLD component spectral properties and underlying physiological mechanisms. The variation of BOLD DOF with progressively increased anesthesia dosage confirms prior suggestions that the use of anesthesia in fMRI experiments requires careful consideration. A number of studies have indicated to date that the type and dosage of anesthesia has a significant impact on spontaneous BOLD signals and the analysis of functional connectivity. Domitor (medetomidine), for example, has been shown to increase venous contrast at ultra high field potentially due to decreased cerebral blood flow (Ciobanu et al., 2012). Isoflurane anesthesia may decrease BOLD signal amplitude by decreasing metabolism and reducing venous deoxyhemoglobin concentrations (Kannurpatti et al., 2008). In experiments that combine medetomidine and isoflurane, complex vascular effects could also arise, such as the vasoconstrictor effect of medetomidine counteracting the vasodilator effect of isoflurane (Fukuda et al., 2013). In this study, we showed that increased isoflurane dosage reduces the effective DOF of spatiotemporal BOLD fluctuations. This effect could be explained by decreased metabolism leading to a bulk decrease in neuronal activation and thus inter-regional communication, altogether decreasing variability or complexity of activity as represented in spatiotemporal BOLD DOF. Furthermore, a strong dose effect was noted in the reduction of BOLD DOF with the increase in anesthesia. The BOLD DOF at the 1.0% dosage were substantially higher than expected—on a par with human data—leading to component maps with remarkable localization of TE-dependent fluctuations to well-known rat neuroanatomy. We note that for our non-randomized isoflurane anesthesia design, the accumulation of isoflurane in fat may have resulted in a greater effective anesthesia dose over time. This increased effective dose may be reflected in our analysis of BOLD DOF versus anesthesia level. Variations of BOLD DOF may also be related to physiological parameters that vary in response to anesthesia and affect hemodynamics such as such as blood pressure, respiration, temperature, and blood pCO2. Global statistical properties of the BOLD signal subspace, such as BOLD DOF, may ultimately respond to many such physiological parameters. This phenomenon draws attention to the need for careful assessment of global statistical properties in studies of functional connectivity in animal models. More quantitative assessment of BOLD DOF and connectivity across various anesthesia levels and in response to varying physiological parameters make for a wide variety of further studies on

873

hemodynamics and neurovascular coupling. Altogether, ME-fMRI in combination with ME-ICA is indicated to be an effective tool for studies of functional neuroanatomy in animal models and neurovascular responses to pharmacological manipulations such as anesthesia. Acknowledgments P.K. was supported by the NIH-Cambridge Scholars Program. Data analysis was conducted in part using the NIH Helix and Biowulf resources. We also thank IHU and IPRIAC for their support for the 11.7 T MRI system. ETB is employed half-time by the University of Cambridge and half-time by GlaxoSmithKline (GSK); he holds stock in GSK. References Beckmann, C.F.,Smith, S.M., 2004. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans. Med. Imaging 23 (2), 137–152. Bianciardi, M.,Fukunaga, M.,van Gelderen, P.,Horovitz, S.,de Zwart, J.,Shmueli, K.,Duyn, J., 2009. Sources of functional magnetic resonance imaging signal fluctuations in the human brain at rest: a 7 T study. Magn. Reson. Imaging 27, 1019–1029. Biswal, B.B., Ulmer, J.L., 1999. Blind source separation of multiple signal sources of fMRI data sets using independent component analysis. J. Comput. Assist. Tomogr. 23 (2), 265–271. Biswal, B.,Yetkin, F.,Haughton, V.,Hyde, J., 1995. Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537–541. Biswal, B., Mennes, M., Zuo, X., Gohel, S., Kelly, C., Smith, S., Beckmann, C., Adelstein, J., Buckner, R., Colcombe, S., 2010. Toward discovery science of human brain function. Proc. Natl. Acad. Sci. 107, 4734. Borsook, D.,Becerra, L.,Hargreaves, R., 2006. A role for fMRI in optimizing CNS drug development. Nat. Rev. Drug Discov. 5, 411–425. Bosshard, S.C.,Baltes, C.,Wyss, M.T.,Mueggler, T.,Weber, B.,Rudin, M., 2010. Assessment of brain responses to innocuous and noxious electrical forepaw stimulation in mice using BOLD fMRI. Pain 151 (3), 655–663. Ciobanu, L.,Reynaud, O.,Uhrig, L.,Jarraya, B.,Le Bihan, D., 2012. Effects of anesthetic agents on brain blood oxygenation level revealed with ultra-high field MRI. PLoS ONE 7 (3), e32645. Comon, P., 1994. Independent component analysis, a new concept? Signal Process. 36 (3), 287–314. Cox, R., 2012. AFNI: what a long strange trip it's been. NeuroImage 62, 743–747. Daubechies, I., Roussos, E., Takerkart, S., Benharrosh, M., Golden, C., D'ardenne, K., Richter, W., Cohen, J.D.,Haxby, J., 2009. Independent component analysis for brain fMRI does not select for independence. Proceedings of the National Academy of Sciences 106 (26), 10415–10422. Fransson, P., 2005. Spontaneous low-frequency BOLD signal fluctuations: an fMRI investigation of the resting-state default mode of brain function hypothesis. Hum. Brain Mapp. 26, 1529. Fukuda, M., Vazquez, A.L., Zong, X., Kim, S.G., 2013. Effects of the α2‐adrenergic receptor agonist dexmedetomidine on neural, vascular and BOLD fMRI responses in the somatosensory cortex. Eur. J. Neurosci. 37 (1), 80–95. Hutchison, R.M.,Mirsattari, S.M.,Jones, C.K.,Gati, J.S.,Leung, L.S., 2010. Functional networks in the anesthetized rat brain revealed by independent component analysis of restingstate fMRI. J. Neurophysiol. 103, 3398–3406. Hyvarinen, A., 1999. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural Netw. 10, 626–634. Jonckers, E.,Delgado y Palacios, R.,Shah, D.,Guglielmetti, C.,Verhoye, M.,Van der Linden, A., 2013. Different anesthesia regimes modulates the functional connectivity outcome in mice. Magn. Reson. Med. http://dx.doi.org/10.1002/mrm.24990. Kalthoff, D.,Seehafer, J.U.,Po, C.,Wiedermann, D.,Hoehn, M., 2011. Functional connectivity in the rat at 11.7 T: impact of physiological noise in resting state fMRI. NeuroImage 54, 28282839. Kannurpatti, S.S., Biswal, B.B., Kim, Y.R., Rosen, B.R., 2008. Spatio-temporal characteristics of low-frequency BOLD signal fluctuations in isoflurane-anesthetized rat brain. NeuroImage 40, 17381747. Keilholz, S.D., Silva, A.C., Raman, M., Melke, H., Koretsky, A.P., 2004. Functional MRI of the rodent somatosensory pathway using multislice echo planar imaging. Magn. Reson. Med. 52 (1), 89–99. Kundu, P.,Inati, S.J.,Evans, J.W.,Luh, W.M.,Bandettini, P.A., 2012. Differentiating BOLD and non-BOLD signals in fMRI time series using multi-echo EPI. Neuroimage 60 (3), 1759–1770. Kundu, P., Brenowitz, N.D., Voon, V., Worbe, Y., Vrtes, P.E., Inati, S.J., Saad, Z.S., Bandettini, P.A., Bullmore, E.T., 2013. Integrated strategy for improving functional connectivity mapping using multiecho fMRI. Proc. Natl. Acad. Sci. 201301725. Menon, R.,Ogawa, S.,Tank, D.,Ugurbil, K., 1993. Tesla gradient recalled echo characteristics of photic stimulation-induced signal changes in the human primary visual cortex. Magn. Reson. Med. 30, 380. Ogawa, S., Lee, T.M., Kay, A.R., Tank, D.W., 1990. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc. Natl. Acad. Sci. U. S. A. 87, 9868–9872. Pawela, C.P., Biswal, B.B., Cho, Y.R., Kao, D.S., Li, R., Jones, S.R., Schulte, M.L., Matloub, H.S., Hudetz, A.G., Hyde, J.S., 2008. Resting-state functional connectivity of the rat brain. Magn. Reson. Med. 59, 1021–1029.

874

P. Kundu et al. / NeuroImage 102 (2014) 861–874

Pawela, C.P.,Biswal, B.B.,Hudetz, A.G.,Schulte, M.L.,Li, R.,Jones, S.R.,Cho, Y.R.,Matloub, H.S., Hyde, J.S., 2009. A protocol for use of medetomidine anesthesia in rats for extended studies using task-induced BOLD contrast and resting-state functional connectivity. Neuroimage 46 (4), 1137–1147. Paxinos, G., Watson, C., 2006. The Rat Brain in Stereotaxic Coordinates: Hard Cover Edition. Academic Press.

Taylor, J., Perez, F., 2011. NIPY: Neuroimaging in Python. Tipping, M.E.,Bishop, C.M., 1999. Probabilistic principal component analysis. J. R. Stat. Soc. Ser. B Stat Methodol. 61, 611–622. Williams, K.A., Magnuson, M., Majeed, W., LaConte, S.M., Peltier, S.J., Hu, X., Keilholz, S.D., 2010. Comparison of α-chloralose, medetomidine and isoflurane anesthesia for functional connectivity mapping in the rat. J. Magn. Reson. Imaging 28 (7), 995–1003.