# **INVESTIGATING THE HUMAN BRAINSTEM WITH STRUCTURAL AND FUNCTIONAL MRI**

# **Topic Editors Florian Beissner and Simon Baudrexel**

## *FRONTIERS COPYRIGHT STATEMENT*

© Copyright 2007-2014 Frontiers Media SA. All rights reserved.

All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.

**ISSN** 1664-8714 **ISBN** 978-2-88919-247-2 **DOI** 10.3389/978-2-88919-247-2

## *ABOUT FRONTIERS*

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

## *FRONTIERS JOURNAL SERIES*

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing.

All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

## *DEDICATION TO QUALITY*

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view.

By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

## *WHAT ARE FRONTIERS RESEARCH TOPICS?*

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area!

Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# **INVESTIGATING THE HUMAN BRAINSTEM WITH STRUCTURAL AND FUNCTIONAL MRI**

## Topic Editors:

**Florian Beissner,** Hannover Medical School, Hannover, Germany **Simon Baudrexel,** University Hospital, Frankfurt am Main, Germany

A compilation of key results from the nine articles in this Research Topic (with permission of the authors).

The brainstem is one of the least understood parts of the human brain despite its prime importance for the maintenance of basic vital functions. Owing to its role as a relay station between spinal cord, cerebellum and neocortex, the brainstem contains vital nodes of all functional systems in the central nervous system, including the visual, auditory, gustatory, vestibular, somatic and visceral senses, and the somatomotor as well as autonomic nervous systems. While the brainstem has been extensively studied in animals using invasive methods, human studies remain scarce. Magnetic resonance imaging (MRI) as a non-invasive and widely available method is one possibility to access the brainstem in humans and measure its structure as well as function. The close vicinity of the brainstem to large arteries

and ventricles and the small size of the anatomical structures, however, place high demands on imaging as well as data analysis methods. Nevertheless, the field of brainstem-(f)MRI has significantly advanced in the past few years, largely due to the development of several new tools that facilitate studying this critical part of the human brain. Within this scope, the goal of this Research Topic is to compile work representing the state of the art in functional and structural MRI of the human brainstem.

# Table of Contents

Luke A. Henderson and Vaughan G. Macefield

Christian Lambert, Rumana Chowdhury, Thomas H. B. FitzGerald, Stephen M. Fleming, Antoine Lutti, Chloe Hutton, Bogdan Draganski, Richard Frackowiak and John Ashburner

*70 Change in Brainstem Gray Matter Concentration Following a Mindfulness-Based Intervention is Correlated With Improvement in Psychological Well-Being*

Omar Singleton, Britta K. Hölzel, Mark Vangel, Narayan Brach, James Carmody and Sara W. Lazar

*77 The Ascending Reticular Activating System From Pontine Reticular Formation to the Thalamus in the Human Brain*

Sang Seok Yeo, Pyung Hun Chang and Sung Ho Jang

*82 Imaging White Matter in Human Brainstem* Anastasia A. Ford, Luis Colon-Perez, William T. Triplett, Joseph M. Gullett, Thomas H. Mareci and David B. FitzGerald

## Investigating the human brainstem with structural and functional MRI

## **Florian Beissner 1,2\* and Simon Baudrexel 3**

1 Department of Neuroradiology, Somatosensory and Autonomic Therapy Research, Hannover Medical School, Hannover, Germany

2 Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital, Charlestown, MA, USA

3 Department of Neurology, Goethe University Frankfurt, University Hospital, Frankfurt am Main, Germany

## **Edited and Reviewed by:**

Hauke R. Heekeren, Freie Universität Berlin, Germany

\*Correspondence<EMAIL_ADDRESS>
**Keywords: brainstem, fMRI, MRI, physiological noise, autonomic nervous system, neuromodulatory systems, pain modulation, reticular formation**

The brainstem is one of the least understood parts of the human brain despite its prime importance for the maintenance of basic vital functions. Owing to its role as a relay station between spinal cord, cerebellum, and neocortex, the brainstem contains vital nodes of all functional systems in the central nervous system, including the visual, auditory, gustatory, vestibular, somatic, and visceral senses, and the somatomotor as well as autonomic nervous systems. The brainstem also contains cholinergic, dopaminergic, noradrenergic, and serotonergic nuclei whose cortical and subcortical projections are essential to the regulation of arousal, behavior, and cognition. Despite this indisputable importance, the brainstem is still largely neglected in attempts to measure or model brain function, especially in human neuroscience. One reason for this neglect is that the anatomical characteristics of the brainstem, specifically its close vicinity to large arteries and ventricles, and the small size of its anatomical substructures, present inherent challenges to neuroimaging analysis. These properties make the brainstem a difficult structure to study with non-invasive methods like magnetic resonance imaging (MRI), as they place high demands on image acquisition as well as data analysis methods. Nevertheless, the field of brainstem-(f)MRI has significantly advanced in the past few years, largely due to the development of several new tools that facilitate studying this critical part of the human brain.

Within this scope, the current goal of this research topic is to compile work representing the state of the art in functional and structural MRI of the human brainstem. We have assembled articles from a number of scientists who have made important contributions to this evolving field, and continue to shape it.

The functional section starts with a review by Brooks et al. (2013) that covers the central problem of physiological noise and presents strategies to suppress it. Ritter et al. (2013) have studied the nociceptive system and show its differential reactions to painful skin heating at different slopes. Ress and Chandrasekaran (2013) use the advantages of ultrahigh magnetic field strengths to study the substructure of the inferior colliculus and its tonotopic organization. The functional section concludes with Henderson and Macefield (2013) who provide a review of their extensive research on somatosensory and autonomic centers in the lower brainstem.

The structural section begins with an article by Deistung et al. (2013) introducing quantitative susceptibility mapping as a new means to boost the identification of anatomical details in structural MRI images. Lambert et al. (2013) use quantitative MRI and tensor based morphometry in a large study sample to characterize aging in the human brainstem. Singleton et al. (2014) apply volumetric methods to demonstrate gray matter changes related to meditation and mindfulness-based intervention. Yeo et al. (2013) use probabilistic fiber tracking on diffusion-weighted images to delineate the ascending reticular activating system. Lastly, by using diffusion tensor imaging at ultra high field strengths, Ford et al. (2013) demonstrate precise tractography results of the human brainstem.

The wealth of methods and applications covered by the authors indicates that functional and structural brainstem-MRI methods have developed to a point where they can be applied to study of a wide range of neuroscientific problems. It is the hope of the editors that the brainstem will soon lose its label of a *terra incognita* and become a region of major interest in the neuroscience community.

## **REFERENCES**

*Received: 23 January 2014; accepted: 17 February 2014; published online: 28 February 2014.*

*Citation: Beissner F and Baudrexel S (2014) Investigating the human brainstem with structural and functional MRI. Front. Hum. Neurosci. 8:116. doi: 10.3389/fnhum.2014.00116*

*This article was submitted to the journal Frontiers in Human Neuroscience.*

*Copyright © 2014 Beissner and Baudrexel. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Physiological noise in brainstem fMRI

#### **Jonathan C.W. Brooks 1\*, Olivia K. Faull 2 , Kyle T. S. Pattinson2 and Mark Jenkinson2**

1 Clinical Research and Imaging Centre, University of Bristol, Bristol, UK

2 FMRIB Centre, Nuffield Department of Clinical Neuroscience, University of Oxford, Oxford, UK

## **Edited by:**

Florian Beissner, Martinos Center for Biomedical Imaging, USA

## **Reviewed by:**

Florian Beissner, Martinos Center for Biomedical Imaging, USA Julien Cohen-Adad, Ecole Polytechnique de Montreal, Canada

## **\*Correspondence:**

Jonathan C. W. Brooks, Clinical Research and Imaging Centre, 60, St Michael's Hill, University of Bristol, Bristol BS2 8DX, UK e-mail<EMAIL_ADDRESS>
The brainstem is directly involved in controlling blood pressure, respiration, sleep/wake cycles, pain modulation, motor, and cardiac output. As such it is of significant basic science and clinical interest. However, the brainstem's location close to major arteries and adjacent pulsatile cerebrospinal fluid filled spaces, means that it is difficult to reliably record functional magnetic resonance imaging (fMRI) data from.These physiological sources of noise generate time varying signals in fMRI data, which if left uncorrected can obscure signals of interest. In this Methods Article we will provide a practical introduction to the techniques used to correct for the presence of physiological noise in time series fMRI data.Techniques based on independent measurement of the cardiac and respiratory cycles, such as retrospective image correction (RETROICOR, Glover et al., 2000), will be described and their application and limitations discussed. The impact of a physiological noise model, implemented in the framework of the general linear model, on resting fMRI data acquired at 3 and 7T is presented. Data driven approaches based such as independent component analysis (ICA) are described. MR acquisition strategies that attempt to either minimize the influence of physiological fluctuations on recorded fMRI data, or provide additional information to correct for their presence, will be mentioned. General advice on modeling noise sources, and its effect on statistical inference via loss of degrees of freedom, and nonorthogonality of regressors, is given. Lastly, different strategies for assessing the benefit of different approaches to physiological noise modeling are presented.

**Keywords: brainstem, physiological noise, fMRI, imaging, 7T**

## **INTRODUCTION**

The human brainstem occupies a key position in the central nervous system (CNS), and is an integral part of the parasympathetic and sympathetic nervous system. The brainstem receives convergent input from spinal and supra-spinal fibers and the cranial nerves, and effectively integrates these signals and coordinates behavioral and physiological responses. As such, this area is of key interest to clinicians and neuroscientists. However, the brainstem is located in an area that suffers from inherently poor signal to noise, due to the close proximity of bone and air-filled cavities, and potent sources of physiological noise. Consequently, functional neuroimaging of this area is problematic. One approach to improving our ability to detect the small blood flow changes associated with neuronal activity, is to account (and correct) for signal changes associated with physiological process that drive signal variation, but may not be related to the neuronal signals of interest.

## **PHYSIOLOGICAL NOISE**

Physiological noise is generally defined as signal changes in an image that are caused by the subject's physiology but excludes brain activity of interest (Jezzard, 1999). This excludes scannerrelated artifacts, such as ghosting or drift, but can include changes related to rigid-body head motion, although we will not consider these motion changes in the present article. The bulk of the physiological noise signals that we will consider arise from cardiac and respiratory processes.

## **SOURCES OF PHYSIOLOGICAL NOISE**

Both the cardiac cycle and the respiratory cycle induce changes in the brain that are detected by MRI. The mechanisms include those due to the cardiac cycle, which induce changes in cerebral blood flow (CBF), cerebral blood volume (CBV), arterial pulsatility, and CSF flow (Greitz et al., 1993; Purdon and Weisskoff, 1998; Dagli et al., 1999; Krüger and Glover, 2001). Mechanisms related to the respiratory cycle include induced changes in the main magnetic field (*B*0) (Raj et al., 2001), and changes in arterial CO2 partial pressure (Wise et al., 2004). In addition, there are interactions between the cardiac and respiratory systems, such as increased cardiac output during inspiration (the "respiratory pump"; Lin, 1999; Hayen et al., 2012) and changes in CSF flow (Klose et al., 2000).

The way in which these mechanisms affect the MRI signal also varies. Some changes are due to real movements, such as pulsatility, while some others lead to apparent movement due to varying geometric distortion, associated with changing *B*0 fields that are induced by the variable air cavity in the chest. In either case the displacements of voxels in the image cause changes that are particularly strong and noticeable near boundaries between areas with substantially different intensities in the image, such as at the edge of the brainstem (see **Figure 1**). Other mechanisms create changes in local blood susceptibility, either through oxygenation changes or blood volume changes. These changes will induce signal variations via the BOLD effect, in the same way as those induced by brain activations of interest. There are also mechanisms that affect the signal via changes in tissue composition or via inflow effects.

Many of these noise sources are stronger in the brainstem than in any other part of the brain, due to the increased capacity for movement, due to CSF flow and blood pulsatility, and the closer proximity to the lungs, and hence the stronger induced changes in the *B*0 field. Therefore the magnitude and composition of physiological noise in the brainstem is quite different from that encountered in other areas of the brain. The scope of this problem can be easily visualized when looking at a map of the temporal signal to noise ratio (tSNR), which is the ratio of the mean of the time course signal intensity divided by its temporal standard deviation (Parrish et al., 2000). It is clear from **Figure 1** that there is an appreciable reduction of tSNR in the brainstem compared to the cerebrum and cerebellum, caused by lower mean signal and/or increased signal variance in the brainstem.

## **CHARACTERISTICS OF PHYSIOLOGICAL NOISE**

## **Field strength dependence**

Physiological noise increases with the square of the main field strength, whilst signal to noise only increases linearly (Parrish et al., 2000; Krüger and Glover, 2001; Triantafyllou et al., 2005; Hutton et al., 2011). This means that for higher field scanners (e.g., at 7 T) physiological noise can become the dominant source of noise, since it scales in the same way as BOLD contrast. Hence the temporal SNR for functional magnetic resonance imaging (fMRI), and consequently statistical power, does not necessarily improve on higher field scanners, especially in areas where physiological noise is already strong, such as in the brainstem. However, higher field scanners do offer other advantages when scanning the brainstem, such as the ability to have increased spatial resolution, which can be invaluable when investigating small structures.

## **Tradeoffs versus thermal noise and resolution**

In MRI, thermal noise is an ever-present source of noise generated by thermal fluctuations within the subject, and to a lesser extent within the receiver electronics, that gives rise to primarily Gaussian-distributed, additive fluctuations in the received signal. When reconstructing images, the distribution of this noise is usually altered, becoming Rician or non-central Chi distributed (Aja-Fernández et al., 2011), although the noise is still reasonably well approximated by a Gaussian noise process in areas of sizeable mean signal magnitude and good tSNR, which are the areas of interest in most fMRI studies. However, whilst the use of array coils and parallel imaging may alter these noise processes, e.g., to a non-central Chi distribution (Breuer et al., 2009), in many cases the Gaussian approximation remains sufficiently good for statistical modeling of the fMRI time series. Further discussion of these issues is beyond the scope of this article.

Although the field strength dependence of thermal-noisedriven SNR is weaker than that of physiological-noise-driven SNR, it is common to scan with higher spatial resolution on high field scanners, which then increases the relative contribution of thermal noise, since thermal SNR scales linearly with the voxel volume. For example, consider 3 mm isotropic fMRI scans acquired at 3 and 7 T. In the 7 T scan the ratio of physiological noise to thermal noise will be 5.4 times higher than in the 3 T scan. However, reducing the voxel volume by a factor of 5.4 times, making it ∼1.7 mm will keep the ratio of physiological noise to thermal noise the same as it was in the 3 mm scan at 3 T. Increasing the resolution further, say to a 1 mm voxel size at 7 T, would then make the ratio (physiological noise to thermal noise) five times smaller in the 7 T scan. Therefore, combined changes in resolution and field strength could mean that the relative ratio of thermal noise and physiological noise was quite different, in either direction, depending on the actual resolutions. As a rough rule of thumb, the contributions of thermal and physiological noise are similar for *cortical* voxels with a 1 mm × 1 mm × 3 mm resolution at 7 T (Triantafyllou et al., 2005). However, in the brainstem, which typically suffers from reduced coil sensitivity due to the distance from (and size of) individual coil elements, there may be an increased contribution of thermal noise at high resolution (as can be seen on **Figure 3**). This makes it important to check on pilot scans (e.g., by measuring temporal SNR) with different resolutions, even at the same field strength.

## **Effect of imaging sequence**

The above comparisons of physiological and thermal noise are all based on using BOLD sensitive fMRI sequences. However, for ASL-based sequences the contributions are likely to be substantially different, since ASL seeks to separate out contributions of blood oxygenation from CBF and CBV. The different sources of physiological noise would therefore contribute quite differently in ASL (see Restom et al., 2006), with oxygenation-related changes being suppressed, blood flow/volume changes being similar or enhanced, and movement-related changes, either real movement due to pulsatility or field-induced image-motion from *B*0 changes, having a similar strength. Currently BOLD fMRI has been exclusively employed in studies focused on the brainstem, due to the inherently stronger basic SNR, but the use and characterization of ASL for brainstem imaging deserves further investigation.

## **CORRECTION METHODS FOR PHYSIOLOGICAL NOISE ACQUISITION STRATEGIES**

Several different approaches have been proposed to minimize the influence of physiological noise (including movement) on acquired MR images using acquisition-based strategies. The different techniques broadly fall into three classes: (i) those that attempt to capture images at a fixed point in the physiological process generating noise, e.g., cardiac gating, cardiac-gated multi-echo (ME) acquisitions; (ii) those that attempt to correct images for intensity variation induced by physiological fluctuation, e.g., un-gated ME acquisitions; and (iii) those that acquire additional scans to help identify physiological noise sources from an independent acquisition, i.e., calibration scans. These are discussed below.

## **Cardiac (or respiratory) gating**

Cardiac signals appear to be the dominant source of physiological noise in fMRI data obtained from the brain, brainstem, and spinal cord (Dagli et al., 1999; Piché et al., 2009). Previously researchers have attempted to minimize the influence of such signals by recording fMRI data at a fixed point relative to the cardiac cycle, i.e., "gating" (Guimaraes et al., 1998; Malinen et al., 2006). An unavoidable consequence of gating is that the time between consecutive fMRI volume acquisitions is no longer governed by a fixed repetition time (TR) and instead has, e.g., a cardiac cycle dependence, implying that there will be different amounts of longitudinal (T1) relaxation between samples. One approach to this problem is to attempt to correct for the amount of partial saturation of MR signal, via an adjustment based on the time between recorded volumes – the *effective* TR (Guimaraes et al., 1998; Malinen et al., 2006). These methods relying on taking a measurement of the apparent T1 relaxation time for each voxel, and the time between adjacent samples (effective TR) to increase or decrease the measured MR signal appropriately. Whilst this is a potentially attractive technique to use when dealing with structures that suffer from significant cardiac pulsatility, the model relies on perfect registration of adjacent images when determining the apparent T1, and assumes a single relaxation time per voxel, which may be erroneous in the case of partial voluming and imperfect realignment. Another potential confound with these approaches is that as the amount of T1 relaxation between volumes (and therefore measured MR signal) is now dependent on the heart rate (HR), so any stimulus that changes this directly or indirectly (e.g., pain, arousal, fear, hypercapnia, blood pressure) may produce a systematic change in signal intensity that may be independent of the neuronal mechanisms associated with BOLD. The potential for introducing bias in the measurement due to imperfect T1 correction, stimulation-coupled HR changes, and the additional processing required prior to analysis, have led to limited use of this approach.

An alternative approach to correcting for T1 changes induced by cardiac gating is to acquire gated imaging data with multiple echoes per TR (Zhang et al., 2006). By acquiring a minimum of two echoes per TR it is possible to derive an image whose intensity is independent of variable T1 saturation effects, and instead reflects T2∗ . The image is created by either computing the quotient of the two images (per TR) or by calculating the apparent T2∗ using a simple mathematical operation (Beissner et al., 2010). When using this technique, best results have been obtained with careful pre-processing of individual images (e.g., pre-smoothing input data) and optimal normalization to standard templates (Beissner et al., 2011). One potential drawback of this approach is that when modeling evoked BOLD activity obtained using this acquisition method, the time between samples (which is normally fixed) is variable and dependent on the HR, thus an effective (mean) TR is typically used (Beissner et al., 2011). Note that this is primarily a limitation of the statistical modeling packages used for fMRI analysis, which expect data to be sampled at regular intervals. The use of an effective TR and its effect on modeling the BOLD response, and the impact on temporal autocorrelation correction, has yet to be determined. However, it is worth noting that the errors introduced by using an effective TR are expected to be small when used in conjunction with a block design fMRI experiment, but are likely to be more problematic with event related designs.

## **Acquisition-based k- and image-space corrections**

By altering sequence parameters, or using modified pulse sequences it is possible to reduce or correct for the presence of physiologically induced image artifacts. It is worth noting that at the relatively long echo times (TE) required for BOLD imaging [TE ≈ T2∗ deoxy(blood)], imaging data are sensitive to signal dropout and BOLD-like signal changes induced by physiological processes. It is difficult to mitigate the problems related to signal dropout whilst maintaining sensitivity to BOLD. The echo planar (EPI) acquisitions typically used for BOLD imaging are also sensitive to accumulation of phase errors due to variation in static magnetic field (*B*0), which is affected by physiological processes (e.g., movement of the lungs). One approach to minimizing these influences was proposed by Pfeuffer et al. (2002), who used a navigator echo to correct for differences in zero-order phase caused by movement of the thorax and abdomen. The changing volume of air in the lungs and position of the tissues surrounding them will induce a time varying (at the breathing rate) change in the main magnetic field (*B*0), and the effect of this is to produce a shift in the phaseencoding direction of EPI data (Frank et al., 2001; Raj et al., 2001; Windischberger et al., 2002). By acquiring a navigator echo before the image read-out, and comparing it to the central line of the main k-space data, it is possible to adjust the phase of the acquired data and remove geometric shifts in the phase-encoding direction. It should be noted that this correction will apply uniformly to the imaged slice, and the same effect could be achieved by retrospective motion correction. However, in the case of a segmented EPI acquisition, such a correction will dramatically reduce ghosting in the image.

By combining recent developments in parallel imaging and multi-channel array coils, it is also possible to increase the contrast to noise of EPI data by acquiring ME EPI data (Poser et al., 2006). This approach makes use of the reduced echo train lengths offered with parallel imaging, to (1) reduced susceptibility induced distortions, (2) acquire additional echoes for the same TR, and (3) reduce the apparent through-plane dropout by optimally combining (through weighting each TE image by the contrast to noise ratio) the signal from, in this case, four echoes acquired with an acceleration factor of 2. The advantage of this technique is that non-BOLD contributions to physiological noise, that do not demonstrate strong correlations between the signal measured at each TE, will be reduced via signal averaging. Whilst one might expect physiological processes with typical frequencies of ∼1 Hz (cardiac) and ∼0.3 Hz (respiratory) to demonstrate

a strong temporal autocorrelation between the samples in a ME sequence (typically ∼15 ms apart), their contribution appears to not be strongly correlated, perhaps due to local dephasing and blood flow velocity mechanisms, hence echo averaging will be of benefit (Poser and Norris, 2009). Indeed, ME imaging techniques have been used to good effect to minimize physiological signal fluctuation in and around the brainstem (Kundu et al., 2012).

## **Calibration scans**

Whilst it may be possible to minimize contributions from physiological noise to acquired fMRI data, through alternate acquisition strategies, it may also be beneficial to identify physiological signals from "calibration scans" and remove them during post-processing of task fMRI scans. Such an approach was proposed by de Zwart et al. (2008), which aims to identify non-task related correlations within a region of interest (e.g., the brainstem) from an additional scan where no task is performed.  During the "real" analysis the time course of signal from the previously determined mask is extracted, orthogonalized (Jezzard et al., 2003) to the experimental design, and included in the GLM to determine activity within the region of interest. Whilst this approach was shown to increase *t*-scores by around 10% in the target brain area, the main disadvantages to such an approach are (1) the need to acquire an additional calibration scan, and (2) the assumption that if the task does produce correlated activity within the predefined mask, then it is safe to consider these signals as noise and remove them from the analysis. Different approaches to achieve the same goal have been proposed, e.g., by using independent component analysis (ICA, see Pre-Processing Strategies) to identify noise sources within fMRI data prior to model estimation using the GLM (Xie et al., 2012).

## **PRE-PROCESSING STRATEGIES**

There are a number of strategies for reducing physiological noise that occur after the acquisition but before statistical estimation is performed. These are the strategies that we are classifying as pre-processing strategies.

## **Retrospective correction of k-space data ("RETROKCOR")**

This strategy, originally introduced by Hu et al. (1995), works directly with k-space data, as opposed to the related technique RETROICOR, that works in image-space. RETROKCOR involves correcting the k-space data after acquisition, but prior to reconstruction, by regressing out signals related to the timing and amplitude of the respiratory and cardiac cycles. Therefore it requires the ability to save the k-space data, process this and then apply post-acquisition reconstruction to the processed data.

The principle used here, and in the related method RETROICOR, is to take independent physiological measurements (typically from a pulse-oximeter and respiratory bellows) and create regressors that are used to remove correlated signals from the data. The regressors are based on the timing of the physiological cycles and assuming a periodic but flexible shape. In addition, the respiratory regressors include amplitude modulation, based

on the depth of respiration. More information on the creation of regressors can be found in Section "RETROICOR."

It has been found that RETROKCOR did not perform as well as RETROICOR for standard acquisitions (Glover et al., 2000), but this has not been tested over all types of acquisitions and brain areas. Therefore it is possible that this may be more suitable for methods such as volume-based acquisitions where different k-space lines are associated with distinct times, as opposed to multi-slice images where different slices are associated with different acquisition times. Further investigation of the relative merits of this technique in such circumstances is warranted.

## **Temporal filtering**

The cardiac cycle and respiratory cycle have fairly well defined frequencies, being around 1 Hz for cardiac and 0.2 to 0.3 Hz for respiratory processes. Although these are quite variable between subjects they are usually relatively stable within a given subject. Interaction terms also have frequencies defined by the differences of these frequencies, for example, at 0.7 and 1.3 Hz (based on 1.0 ± 0.3 Hz). Therefore it is theoretically possible to remove signals related to these processes by filtering out these frequencies, or a band of candidate frequencies.

In practice, such a temporal filtering approach suffers from a number of problems: (i) the physiological signals also contain harmonics of the base frequencies; (ii) activations of interest may also contain similar frequencies, especially considering the harmonics; and (iii) at typical TR values the sampling rate is too low to uniquely distinguish these frequencies, and they are effectively aliased to different frequencies. The latter problem is the biggest one for typical sequences, as aliasing can only be avoided if the sampling frequency (1/TR) is at least twice as high as that of the highest frequency in the signal, which is the cardiac signal for physiological noise. For example, it requires the TR to be less than 0.4 s in order to avoid aliasing for cardiac signals up to 75 beats per minute.

Aliasing is a fundamental limitation when taking discrete samples and cannot be avoided (see Aliasing or Nyquist Frequency in any standard textbook on signal processing). This limitation effectively means that at least two samples are required within every period to avoid that signal's frequency being aliased. Otherwise, the signal is undersampled and will have the appearance of a different frequency in the sampled data, which could either be a high frequency, or a low frequency, or something in between. Therefore aliased frequencies are likely to overlap with frequencies from the signals of interest and, consequently, separation by temporal filtering cannot work. However, the emergence of accelerated fMRI acquisitions, such as multi-band fMRI (Feinberg et al., 2010), are making short TR acquisitions feasible without sacrificing brain coverage. With these sequences the use of temporal filtering for short TR acquisitions is becoming a practical option.

Another issue that arises when using temporal filtering, and other pre-processing filtering techniques such as denoising, is correctly accounting for the filtering in the subsequent statistical estimation. Failing to do anything can be problematic if the subsequent statistics are utilizing a form of parametric distribution and making assumptions about the frequency content of the noise. The filtering done here does modify the frequency content of the signal and, in cases where there is substantial filtering (e.g., low pass filtering), this can lead to large parts of the frequency spectrum containing zero power – leading to underestimation of the underlying variance in the data. This is sometimes modeled correctly in parametric statistics, such as is done in certain flexible pre-whitening methods, but would not fit less flexible models such as low-order auto-regressive (AR) models. The implications of this, and also the incorrect statistical Degrees of Freedom (DOF) associated with the estimation, are that the thermal noise influence tends to be underestimated (due to the removal of true signal variance, e.g., with low pass filtering), leading to an increase in false positives. However, it is very difficult to predict the magnitude of the effect on the final results, which may be negligible or substantial, depending on many different aspects of the experimental acquisition, SNR, length, design, etc. To avoid these problems it is often better to employ other, more flexible, statistical methods, such as permutation-based statistics for the higher-level analysis or to perform the equivalent filtering within the statistical estimation, as will be discussed for RETROICOR below (see RETROICOR).

## **Denoising with ICA**

Independent Component Analysis is a method of decomposing a dataset into its constituent sources, taking into account the spatial and temporal structure of the sources. In fMRI it is used in a way that enforces spatial independence between component maps but allows time courses in different components to be arbitrarily similar, or different. It has been shown (Beckmann and Smith, 2004) that ICA decompositions are capable of separating out sources of scanner artifact, physiological noise, and brain activation within fMRI. Furthermore, aliasing is not a problem as ICA is still able to distinguish components based on different spatiotemporal patterns in long TR data (Brooks et al., 2008), and thus has an advantage over many of the other strategies, which rely on temporal information alone.

Denoising with ICA is a strategy that relies on identifying unwanted components, such as physiological noise or scanner artifact, and removing them from the dataset prior to further analysis. The difficulty in applying this in practice lies in identifying which components are unwanted or not. The classification of the components can be done manually, or with automated classification tools.

Manual classification of components is subjective and relies on the experience of the person doing the classification. Typically the person inspects the set of components for each subject, which can range from a small number to hundreds, and decides which are noise/artifact that will then be removed. Invariably there are some difficult components where it is hard to tell, and it is known that ICA can produce components that are a mixture of different true sources when the SNR is limited. In such cases it is best to take the conservative approach and leave such components in the data; that is, do not classify them as noise. This ensures that most of the signals related to brain activation remain in the data, minimizing the risk of classifying the signals of interest as noise.

An alternative to this approach, which is particularly useful when considering physiological noise sources, is to compare the spatial location of putative noise components identified from a short TR resting calibration scan (see Calibration Scans), to those estimated from a long TR experimental acquisition (Brooks et al., 2008; Xie et al., 2012). The study by Xie et al. builds on earlier work (Piché et al., 2009) using spatial ICA to identify cardiac components from resting EPI data acquired from the spinal cord, which were classified as noise components on the basis of their location and correlation with recorded physiological traces. Similarly, Xie et al. recorded resting data at short TR (250 ms) to unambiguously identify cardiac and respiratory components, and then compared the obtained spatial maps with ICA results from task fMRI data (acquired with long TR). The components overlapping the noise signals identified on short TR data, could then be included as nuisance regressors in the GLM. This approach was seen to increase both the sensitivity and specificity of spinal fMRI responses to painful electrical stimulation. The main disadvantages to this approach are (1) the need to acquire resting data, and (2) the additional processing steps required to analyze short (resting) and long TR (task) data, and compare the obtained component maps before inclusion in the GLM.

Automated classification methods (Thomas et al., 2002; Tohka et al., 2008; Churchill et al., 2012; Smith et al., 2013a) tend to be based on machine learning approaches that use a training set to learn how to discriminate between different categories of components, such as physiological noise, motion effects, scanner artifacts, and brain activation. As with manual classification, these methods will not be perfect and some errors in classification will occur; for example, the recent method in (Smith et al., 2013a) had an accuracy of over 99% on data from the Human Connectome Project. If the methods provide some control or measure of confidence, it is better to err on the side of not removing certain components when possible. Furthermore, the automated methods rely heavily on the training data matching well with the acquisition being analyzed, and differences in acquisitions (e.g., FOV, field strength, resolution, etc.) can cause the number of mis-classifications to dramatically increase. Therefore, it is important that the classification be carefully monitored, especially when run for the first time on new data. However, these methods do show strong promise in substantially reducing physiological noise.

As mentioned above, care must be taken in the subsequent statistics to account for any denoising done in the pre-processing to avoid falsely inflated statistics. For moderate amounts of denoising this is unlikely to cause a problem for most higher-level analysis methods, where the between-subject variance dominates, but for more extensive denoising this may cause a bias when parametric statistics are used in the higher-level analysis. The precise quantification of these effects is currently unknown and therefore, non-parametric statistics, because they are robust to distributional changes, are currently the preferred option for higher-level analysis of denoised data.

## **PHYSIOLOGICAL NOISE MODELING**

As already suggested in the discussion of RETROKCOR (Hu et al., 1995), one method for estimating and removing physiological noise sources from time series fMRI data is to acquire independent physiological measurements from the subject and base the correction on these data. Glover et al. (2000) proposed that this correction was optimally performed in image-space, based on 2D

single shot acquisitions using EPI, and was termed RETROspective Image CORrection (RETROICOR).

## **RETROICOR**

Typically cardiac and respiratory processes will be monitored using a pulse-oximeter and respiratory bellows, respectively. These physiological waveforms may be recorded on a separate computer, along with scanner triggers to indicate the timing of each volume acquisition. The task is then to determine the phase (cardiac and respiratory) associated with the timing of each acquired slice in the imaged volume. This process is illustrated in **Figure 2**, see legend for full description.

The phase assigned (varying from zero at the arbitrarily chosen starting point in the cycle, to 2π at the end of the cycle) to each slice, may then be entered into a low-order Fourier expansion (see Glover et al., 2000; Brooks et al., 2008; Harvey et al., 2008), to derive time course regressors that attempt to explain signal changes, which are driven primarily by cardiac or respiratory processes (or their interaction). These signals can then be used to "regress out" physiological signals from the raw time course data as per Glover et al. (2000). Alternatively, the time course of regressors can be included in the GLM – as first suggested by Josephs et al. (1997), whereby the weighting ("beta" or "parameter estimate") of each component will be adjusted to produce the optimal fit to the data. One advantage of the GLM approach is that it explicitly accounts for the loss of DOF that will occur when including large numbers of nuisance regressors. Equally, if the

calculated physiological time courses are not orthogonal to the experimental design (which is likely to be the case), the GLM will apportion the variance between the different regressors, and thus provide a conservative estimate of appropriate statistical quantities, guarding against false positives. This conservative approach is recommended and will automatically happen provided that no user-enforced orthogonalizations of regressors are performed. In relation to the brainstem, Harvey et al. (2008) used hierarchical *F*tests to determine which regressors explained significant amounts of variance in resting data acquired at 3 T. They concluded that three orders of cardiac, four orders of respiratory, and a single set of interaction terms were sufficient to explain variance in their brainstem fMRI data.

## **Heart rate, respiratory rate, and carbon dioxide**

Whilst the cardiac and respiratory phases can explain significant sources of noise in fMRI data, it is also possible to model second order changes associated with variation in the cardiac and respiratory cycles.

Alterations in metabolic rate and ventilation (the product of tidal volume and respiratory rate) have opposing effects upon the partial pressure of carbon dioxide (CO2) present in the arterial blood (PaCO2). CO2, due to its vasodilatory effect upon the cerebral vasculature, has a measurable effect on the BOLD signal (Wise et al., 2004, 2007). One approach to correcting for this signal change is to measure the partial pressure of CO2 in expired breath (PETCO2), and include this measure as an explanatory variable (EV) at the analysis stage (Pattinson et al., 2009a,b). PETCO2 varies tightly with PaCO2 in people with healthy lungs, and is thus considered a valid approximation.

Although it is preferable to directly measure PETCO2, in the absence of the necessary monitoring equipment Birn et al. (2006) have proposed a method to approximate changes in PETCO2 via a surrogate measure calculated from a respiratory bellows trace: the "respiratory volume per unit time" (RVT, Birn et al., 2006). This method assumes that chest expansion represents measured tidal volume, although this is partly correct, respiratory bellows do not account for lung expansion in the vertical plane (i.e., downwards into the abdomen). Furthermore, changes in metabolism (e.g., induced by drugs) may affect PaCO2 in a way not fully explained by changes in respiration.

Additionally, variation in HR has also been shown (Shmueli et al., 2007; Chang et al., 2009) to explain "noise" in brain imaging data, however, as with all regressors care should be taken when removing associated signals from fMRI data that may be correlated with the experimental design. For example, in experiments using painful stimuli, HR and the stimulation timing have been shown to be correlated (Tousignant-Laflamme et al., 2005) and in this case including HR in the model has been shown to impact activation statistics (Kong et al., 2009).

## **Alternative model-based approaches**

Clearly the RETROICOR approach, and its derivatives, depends on the model chosen to approximate the effect of physiological noise on the measured fMRI signal (although it is still reasonably flexible). To avoid constraining the physiological noise model to a particular set of basis functions (e.g., Fourier basis functions in the case of RETROICOR), it is also possible to use finite impulse response (FIR) functions to model physiological time series data (Deckers et al., 2006). The main assumption with this method is that the physiological processes are quasi-periodic and constant amplitude, and that their effects can be modeled by only using the timing of slice acquisition relative to the peaks detected from the pulse and respiratory waveforms. By finding the peaks in the cardiac and respiratory waveforms, each slice is assigned to a particular "bin" (i.e., a fraction of the cardiac or respiratory cycle). By including separate regressors for each time interval (or bin), whose weight is one (1) for those images falling into the relevant time window, one can build up a complete set of regressors which aims to model the cardiac and respiratory signals. Performance was found to be "at least equivalent to the RETROICOR method" (Deckers et al., 2006). However, there are some issues that arise when using this technique, such as determining the optimal total number of bins (which is somewhat arbitrary and can depend on the fMRI acquisition length, see Kong et al., 2012) also, unlike RETROICOR, FIR-based approaches do not account for the depth of breathing, which has been shown to have a dramatic effect on the induced artifact in EPI data (Raj et al., 2001).

## **IMPROVEMENTS IN fMRI TIME SERIES MODELING**

One way to assess possible improvements of the different correction approaches for physiological noise, is to compare the tSNR calculated from the raw data (motion-corrected only) to that after removing physiological noise.

To illustrate the sort of improvement that might be obtained in the brainstem, we acquired resting fMRI data from two healthy subjects, scanned separately at 3 T (Siemens Verio, manufacturer's 32-channel head coil) and at 7 T (Siemens Magnetom 7 T, Nova Medical 32-channel Rx with single channel birdcage Tx). Data were acquired using a rectilinear EPI sequence with the following voxel sizes 3 T (3 mm isotropic), 7 T (2 and 1 mm isotropic), and with 100 time points in each case. The 3 T data were acquired with an axial interleaved slice order, number of slices = 56, field of view = 192 mm, TE/TR = 30/3000 ms, flip angle = 90°, phaseencoding R/L, iPat acceleration factor = 2 (GRAPPA reconstruction), and bandwidth = 2111 Hz/Px. The 7 T data (2 mm) were acquired with coronal ascending slice order, number of slices = 27, field of view = 192 mm, TE/TR = 24/2500 ms, flip angle = 90°, phase-encoding S/I, iPat acceleration factor = 2 (GRAPPA), and bandwidth = 1132 Hz/Px. The 7 T 1 mm data were acquired with an increased number of slices (*n* = 54) to retain coverage of the brainstem, iPat acceleration factor = 3 (GRAPPA), and TR = 5000 ms, all other parameters remained the same. All resting fMRI data were analyzed with and without physiological noise modeling (PNM) within the framework of the general linear model in FEAT (part of FSL software, Jenkinson et al., 2012).

Physiological recordings taken with a pulse-oximeter and respiratory bellows, were logged along with scanner volume triggers at 50 Hz sampling using a BIOPAC MP150 (Goleta, CA, USA), and stored on a computer running Acqknowledge 4.2. Cardiac, respiratory, interaction, HR, and RVT regressors were computed for each slice independently using PNM software (part of FSL), giving a total of 34 regressors. Note that no slice timing correction

was applied to the data, as the PNM takes into account that each slice is acquired at a specific time, and so slice timing correction, together with the additional interpolation that it performs, is not needed. These PNM regressors were input to the GLM along with a single dummy regressor (a requirement for FEAT), and the model fitted. No smoothing, brain extraction, or temporal filtering was applied, however data were pre-whitened (using FILM, part of FSL) as this was expected to have an effect on physiological sources producing temporal autocorrelation within fMRI time series data (Woolrich et al., 2001). Finally the tSNR was calculated for the raw (motion-corrected only) resting fMRI data and for the residuals following GLM estimation (i.e., with physiological components removed). One important consideration is the loss of DOF incurred when modeling with the PNM, which would be expected to change temporal smoothness and reduce variance even if using a set of randomly generated regressors. To address this the temporal standard deviation (tSTD) was calculated by normalizing with the true DOF, which is equal to the number of time points minus the number of regressors minus one (i.e., *N*−1−*N*reg; where *N*reg is equal to 35 in this case – 34 PNM regressors plus 1 dummy regressor).

To visualize where in the brain tSNR (temporal mean divided by temporal standard deviation) is increased using the PNM, the absolute tSNR computed for each acquisition before ("raw") and after correction ("corrected") is shown in **Figure 3**, along with the average value within a hand drawn brainstem mask (white outline shows the position of the brainstem). Improvement in tSNR was found for all three acquisitions, but is more easily visualized on the 1 mm isotropic 7 T data, where variance was particularly reduced around the surface of the pons and near the fourth ventricle.

To provide an indication of the spatial localization and magnitude of improvement in tSNR, the ratio of the corrected to the raw tSNR was computed, and is shown using the voxel-wise maps on the right hand side of **Figure 3**. Whilst there are clearly a significant number of voxels which do not benefit from correction with the PNM, the *average* improvement in tSNR within the brainstem was ∼13% across all three acquisitions.

The improvement seen in 3 and 7 T may be compared to that previously reported in the brainstem at 3 T (Harvey et al., 2008), see **Figure 4**. Harvey et al. studied 12 healthy subjects with a brainstem optimized coronal oblique EPI acquisition (voxel size 2.5 mm × 2.5 mm × 3 mm), TE/TR = 30/1000 ms, flip angle = 70°, to acquire 1130 volumes with simultaneous physiological monitoring (cardiac and respiratory). A series of hierarchical *F*-tests were performed to indicate which physiological regressors (cardiac, respiratory, and interaction) explained significant amounts of noise in their resting data. The final model included three cardiac, four respiratory, and one interaction term (Brooks et al., 2008). The improvement obtained with this model was demonstrated by comparing the temporal coefficient of variation (temporal standard deviation divided by temporal mean, multiplied by 100) between the uncorrected and corrected data. Across the group, temporal CV was reduced most in the medulla, but also near the surface of the pons, near the floor of the fourth ventricle and in mid-brain (particularly around the periaqueductal gray matter, PAG).

**FIGURE 3 | Representative improvements in temporal signal to noise (tSNR) obtained through modeling physiological noise (single subject data)**. On the left the absolute tSNR maps for each acquisition are provided for both pre- ("raw") and post-correction ("corrected"), and the corresponding average tSNR within a brainstem mask is given in brackets. Voxels overlapping with the CSF filled spaces around the brainstem benefit most from physiological noise correction. Maps on the right demonstrate the ratio of the

corrected to raw tSNR, and are thresholded at 1.1 to indicate where one might reasonably expect to see improvement of greater than 10% in the measured tSNR. Clearly the cortex appears to benefit most from this correction with increases in tSNR of 100% frequently observed at all resolutions and field strengths. The improvement in the brainstem is more modest, but nonetheless improved on average by at least 12.5% in this area (for all acquisitions).

## **DISCUSSION**

There are many available options for dealing with physiological noise, as detailed in the previous section. Deciding which options are best depends very much on the particular experiment being undertaken, the subject cohort, the available equipment, the available expertise, the available software, etc. So one recommendation will not suit all situations. Therefore, in this section we will contrast the various advantages and disadvantages of the different options and suggest a method of piloting and comparing various options in practice.

## **EXPERIMENTAL REQUIREMENTS**

There are many possible impacts of physiological correction methods on an experiment. One is the extra time in the scanner session required for setup or acquisition. This includes time taken to attach a pulse-oximeter and respiratory bellows (or equivalent devices) to the subject, and verify that signals are recorded correctly. In particular, subjects' hands should be kept warm to maintain circulation in the fingers, and any nail varnish removed to increase light transmittance; the bellows should be placed carefully near to the diaphragm on the ribcage to maximize changes in recorded lung volume. However, all this can normally be achieved in only a few minutes and has minimal impact on scanning sessions unless timing is very tight. Another impact on scanning time is the addition of new sequences to the session, such as would be required for a resting-state calibration or increased time for certain acquisitions, such as cardiac-gated sequences or certain ME choices. This tradeoff of scan time versus benefit in data quality and scope is a standard problem faced when planning MRI scanning sessions, and will be familiar to experimenters. Evaluating specific benefits of scanning time versus data quality is something that can usually be done with pilot data, as discussed below.

Another requirement for some physiological noise correction methods is extra equipment, such as physiological measuring devices (e.g., pulse-oximeter, respiratory bellows, expired gas analyzer) and digitizing equipment (e.g., BIOPAC, National Instruments, ADInstruments, Cambridge Electronic Devices). When such equipment is not already available, the experimenter must decide whether the additional cost of buying this equipment warrants the benefits to the fMRI analysis. However, given the typical cost of scanning subjects, these pieces of equipment are normally a fairly minor purchase and the monetary cost is rarely a major consideration. Indeed, many MR research centers will typically have patient monitoring apparatus (e.g., *In vivo* Precess, MEDRAD Veris), which can normally output some, if not all, of the required signals (e.g., cardiac triggers) for physiological monitoring. Equally, respiratory monitoring can be achieved with non-magnetic rubber bellows (e.g., from Lafayette Instrument) and commercially available pressure sensors, and data logged with an inexpensive analog to digital converter (e.g., from National Instruments). One essential requirement is that the timing of scan acquisition (i.e., slice or volume triggers) is recorded along with the physiological waveforms using a single computer with a common sampling rate.

Some of the correction methods discussed above also require access to different sequences, such as cardiac-gated sequences, multi-band fMRI, or ME fMRI. These may be more difficult to obtain, being dependent on the type of MR system, requiring input from MR radiographers/operators/physicists and typically a research agreement with the scanner manufacturer. Similarly, other methods require access to raw k-space data and reconstruction algorithms (e.g., RETROKCOR), which can also be difficult to acquire. Therefore, for certain cases some of these correction methods may not be feasible, but given the wide array of options that require nothing more than standard sequences, this is not a major problem.

## **ANALYSIS CHOICES AND IMPLICATIONS**

Given the availability of physiological data and analysis routines capable of generating suitable regressors, which should you choose? In our experience the performance of techniques based on "binning" time series data according to slice and physiological waveform timing (e.g., Deckers et al., 2006), depend critically on the choice of total number of bins (Kong et al., 2012), and do not address issues relating to depth of breathing, which can be critical. ICA approaches offer the possibility to automatically identify sources on physiological noise, however, their usefulness must be judged against the possibility of true signal being identified as a noise source and removed from the data. One approach to mitigate this outcome is to acquire resting and task fMRI data, and use the resting data to define physiological noise components, for subsequent removal from the task data (e.g., Xie et al., 2012).

Concerning the model-based techniques, such as RETROICOR (Hu et al., 1995; Josephs et al., 1997; Glover et al., 2000) and PNM (Brooks et al., 2008; Harvey et al., 2008) these will always be subject to the limitations of the basis functions used to describe the physiological processes they are attempting to model. Further developments that attempt to capture variance induced by changes in HR (e.g., Shmueli et al., 2007; Chang et al., 2009) or respiratory rate and depth (e.g., RVT, Birn et al., 2006) will be necessary to account for effects not accounted for in RETROICOR or PNM. These additional regressors are straightforward to compute, and mostly likely account for low frequency signal fluctuations in time series fMRI data. There are however, additional considerations one must take into account when using these slice dependent model-based approaches. For example, slice timing correction when applied to time series fMRI data could potentially break the association between the physiologically induced signal changes and signal changes predicted by the model. Equally, with increased use of multiplexed parallel imaging, e.g., multi-band (Feinberg et al., 2010) or simultaneous multi-slice imaging (Setsompop et al., 2012), researchers will need to be careful in assigning cardiac and respiratory phases to subsets of their slices acquired at the same time. For 3D or segmented acquisitions, it may be more appropriate to perform corrections on raw k-space data (Tijssen et al., in press). Finally, one must be aware of the penalty in terms of loss of DOF incurred when using these model-based techniques. However, the loss of 40 DOF within the context of a time series of 100 or more time points will constitute a relatively small change in statistics, with potentially greater ability to detect signal against the background noise in the experiment.

## **OPTIMIZING YOUR EXPERIMENT**

Given the many ways in which different experiments and experimental setups can vary it is often necessary to decide on the best strategy for acquisition and analysis for each experiment separately. Collecting specific pilot data is usually the best way to make this decision; however evaluating small pilot experiments is not straightforward. Here we offer some advice on how to collect and evaluate pilot data in order to evaluate different acquisition and analysis choices. It cannot be guaranteed that this will lead to the best possible experiment, but following a procedure like this is very likely to improve the experimental results.

Certain types of experiment can be more demanding in terms of finding effects of interest, such as pain experiments where subjects often move more and there may be correlations between the applied stimulus and physiological changes in respiratory and cardiac cycles (e.g., Tousignant-Laflamme et al., 2005; Kong et al., 2009). Furthermore, for brainstem investigations that involve neural mechanisms that play a role in respiratory control (Smith et al., 2013b) or autonomic regulation (Macefield and Henderson, 2010), it will be very difficult to disentangle BOLDrelated activations and physiological noise and therefore extra care and attention to the analysis is called for (Pattinson et al., 2009a,b).

Once the limitations and difficulties have been thought about and the various pros and cons of the correction methods have been weighed up, it is normally necessary to evaluate a limited number of alternative options. For example, using a cardiac-gated sequence or not, or applying a measurement-driven method like RETROICOR versus ICA denoising. If multiple options are being considered, it is usually beneficial to rank comparisons in order

### **Table 1 | Summary of available physiological noise correction/removal tools.**

Please note that this list should not be considered to be definitive. It is provided to give some indication of the available software, and does not constitute a recommendation or endorsement.

of importance and then conduct a number of separate pairwise comparisons, rather than trying to evaluate the full factorial set of options, since the latter usually ends up taking an enormous amount of time.

When comparing different acquisition methods, it is sufficient to compare the tSNR directly, provided that the basic contrast mechanism is the same; for example, both using BOLD contrast. However, the comparison should not be as simple as "larger is better" when there are changes in pre-processing (e.g., different amounts of spatial smoothing) or statistical estimation. This is because it would lead to the conclusion that increasing the filtering, or number of regressors, or number of components in denoising, is always beneficial, even though this is not the case. Instead, the relative benefit of including extra regressors, or changing the preprocessing options, should be evaluated using techniques such as *F*-tests or information-theory measures like the Bayesian Information Criterion (BIC). These techniques are based on probability theory and take into account the fact that random noise is always soaked up by extra filtering, denoising, or regression, but that removing too much of the thermal noise makes it more difficult to estimate the statistics reliably, thus decreasing the statistical DOF, which would lead to weaker, not stronger, statistics.

Calculation of the *F*-test is straightforward for methods that are embedded within the same statistical estimation framework. For example, testing whether a set of additional regressors in RETROICOR is beneficial can be done using an *F*-test over the extra regressors (e.g., as done in Brooks et al., 2008; Harvey et al., 2008). This will give a statistical result at each voxel that tests whether these extra regressors fit an amount of variance that is beyond what would be expected from thermal noise alone. The overall benefit can then be judged subjectively, based on whether a sufficient number of voxels in the area of interest show a significant statistical result or not.

When comparing differing pre-processing steps it is more difficult to use *F*-tests, and this is where the BIC becomes more useful. For example, the BIC for k regressors and N time points is BIC(*k*, *N*) = *N* × log(RSS/*N*) + *k* × log(*N*), where RSS is the residual sum of squares. Therefore, comparing models with k1 and k2 regressors, which can be from very different sets, yields the condition that log[RSS(*k*1)/RSS(*k*2)] < (*k*2 − *k*1) × log(*N*)/*N*

## **REFERENCES**

*Neuroimage* 52, 524–531. doi:10. 1016/j.neuroimage.2010.04.243

when the model with k1 regressors is considered superior. This test again needs to be performed separately on each voxel and then the overall result can be judged subjectively, based on the pattern of results. In cases where the correction method is not a straightforward regression, an approximate number of regressors can be used, based on the number of denoised components or the rank of the matrix used to perform filtering.

## **CONCLUSION**

In this Methods Article we have presented several different techniques for estimating and removing physiological noise artifacts from time series fMRI data. It is clear that the brainstem suffers from both intrinsically low signal to noise, and increased contributions from non-neuronal physiological sources, which serve to reduce the temporal SNR. The recorded temporal SNR in any given brain region is an important consideration as it will, to a large extent, dictate the length of an fMRI experiment required to measure signals of interest (Murphy et al., 2007). Physiological noise models, and denoising techniques in general, aim to reduce signal variation in the fMRI time series and increase sensitivity to detect activation. From considerations of the effects of field strength and voxel size on recorded tSNR, it is clear that to benefit from the increased BOLD signal afforded by higher field strength MRI systems (3 T and above),one should make a considered choice about the resolution of acquired images. Minimizing voxel size will reduce the influence of physiological noise, but the "intrinsic" signal to noise will be low. This places a large emphasis on choosing the correct voxel dimensions for your experiment, which will inevitably be a compromise between visualizing the small brainstem structures, whilst retaining sensitivity to neuronally induced BOLD changes.

**Table 1** provides the URLs for several different packages used for removing or correcting for the effects of physiological noise.

## **ACKNOWLEDGMENTS**

Kyle T. S. Pattinson is supported by the Medical Research Council (UK), and the National Institute for Health Research (NIHR) Oxford Biomedical Research Centre based at Oxford University Hospitals NHS Trust and University of Oxford.

General formulation for quantitative g-factor calculation in GRAPPA reconstructions. *Magn. Reson. Med.* 62, 739–746. doi:10.1002/mrm. 22066

noise correction on fMRI at 7 T.*Neuroimage* 57, 101–112. doi:10.1016/j. neuroimage.2011.04.018

subjects. *Hum. Brain Mapp.* 31, 539–549. doi:10.1002/hbm.20885

inhomogeneity-desensitized fMRI. *Magn. Reson. Med.* 55, 1227–1235. doi:10.1002/mrm.20900

acquisitions. *Neuroimage.* doi:10. 1016/j.neuroimage.2013.08.062

Fuchsjäger-Mayerl, G., Schmetterer, L., et al. (2002). On the origin of respiratory artifacts in BOLD-EPI of the human brain. *Magn. Reson. Imaging* 20, 575–582. doi:10.1016/ S0730-725X(02)00563-5

linear modelling of FMRI data. *Neuroimage* 14, 1370–1386. doi:10.1006/ nimg.2001.0931

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 17 June 2013; accepted: 09 September 2013; published online: 04 October 2013.*

*Citation: Brooks JCW, Faull OK, Pattinson KTS and Jenkinson M (2013) Physiological noise in brainstem fMRI. Front. Hum. Neurosci. 7:623. doi: 10.3389/fnhum.2013.00623*

*This article was submitted to the journal Frontiers in Human Neuroscience.*

*Copyright © 2013 Brooks, Faull, Pattinson and Jenkinson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Human brain stem structures respond differentially to noxious heat

## **Alexander Ritter \*, Marcel Franz, Caroline Dietrich,Wolfgang H. R. Miltner and ThomasWeiss**

Department of Biological and Clinical Psychology, Friedrich Schiller University, Jena, Germany

## **Edited by:**

Florian Beissner, Martinos Center for Biomedical Imaging, USA

### **Reviewed by:**

Vaughan G. Macefield, University of Western Sydney, Australia Clas Linnman, Boston Children's Hospital, USA

### **\*Correspondence:**

Alexander Ritter, Department of Biological and Clinical Psychology, Friedrich Schiller University, Am Steiger 3/Haus 1, D-07743 Jena, Germany

e-mail<EMAIL_ADDRESS>
Concerning the physiological correlates of pain, the brain stem is considered to be one core region that is activated by noxious input. In animal studies, different slopes of skin heating (SSH) with noxious heat led to activation in different columns of the midbrain periaqueductal gray (PAG). The present study aimed at finding a method for differentiating structures in PAG and other brain stem structures, which are associated with different qualities of pain in humans according to the structures that were associated with different behavioral significances to noxious thermal stimulation in animals. Brain activity was studied by functional MRI in healthy subjects in response to steep and shallow SSH with noxious heat. We found differential activation to different SSH in the PAG and the rostral ventromedial medulla (RVM). In a second experiment, we demonstrate that the different SSH were associated with different pain qualities. Our experiments provide evidence that brainstem structures, i.e., the PAG and the RVM, become differentially activated by different SSH. Therefore, different SSH can be utilized when brain stem structures are investigated and when it is aimed to activate these structures differentially. Moreover, percepts of first pain were elicited by shallow SSH whereas percepts of second pain were elicited by steep SSH. The stronger activation of these brain stem structures to SSH, eliciting percepts of second vs. first pain, might be of relevance for activating different coping strategies in response to the noxious input with the two types of SSH.

**Keywords:A-delta fiber, C-fiber, second pain, pain descriptors, PAG, RVM, periaqueductal grey, rostral ventromedial medulla**

## **INTRODUCTION**

Nociceptive stimulation evokes activity in a number of brain structures including the brain stem. Thereby differential nociceptive stimulation in animals leads to differential activity in the brain stem. Studies in rats indicate that brain stem structures related to nociception like the periaqueductal gray (PAG) and the nucleus Raphe magnus (NRM) in the rostral ventromedial medulla (RVM), are activated differentially by different slopes of skin heating (SSH) (Lumb et al., 2002; Lu et al., 2004; Parry et al., 2008). So the dorsolateral PAG was shown to be preferentially activated in response to steep SSH while activation of the ventrolateral PAG was observed preferentially to shallow SSH (Lumb et al., 2002; Parry et al., 2008). Furthermore, Lu et al. (2004) revealed in rats that activation of the NRM (and nocifensive effects) were different for steep vs. shallow SSH.

Moreover, brain stem activity is directly associated with modulation of pain intensity. So, electrical stimulation of the PAG, one of the brain stem areas usually activated by nociceptive input, has been shown to produce analgesia (Basbaum and Fields, 1984). In animals, distinct brain stem structures have been shown to be associated with distinct behavioral and cardiovascular components of nociceptive reaction. In humans, studies already revealed the importance of brain stem structures for the modulation of pain (Bromm and Treede, 1987b; Behbehani, 1995; Bandler et al., 2000). Recently, placebo analgesia was directly associated with the activity of PAG and RVM (Eippert et al., 2009). This might be of

clinical importance because a specific activation of the brain stem could be associated with a reduction of pain perception. Such a pain-modulation would be interesting especially for chronic pain patients.

Taking into account the above-mentioned activations of PAG and RVM in response to different SSH in animals, we aimed to investigate brain stem activation to two different SSH in humans. However, the paradigm used in animals is difficult to realize in human due to at least two reasons: first, nociceptive stimulation in animals was realized with temperatures up to 60°C for a longer period of time. Such stimulation would cause serious injury in humans (Lumb, 2002; Lu et al., 2004). Second, there were single heat ramps with a delay of 8 min between two stimulations in the animal experiments. This delay is too long even for a block design in functional MRI (fMRI). To our knowledge, there are no studies using different SSH in humans to investigate brain stem activation. However, there are human studies using trains of thermal stimuli with different frequencies (Price et al., 1977; Staud et al., 2007). In these studies, trains of thermal stimuli with different intervals between noxious heat stimulations were applied. Modulating frequency of heat stimulation, humans reported different pain percepts (Price et al., 1977; Price, 1988; Staud et al., 2007) that might be assigned to two distinct conceptual entities, i.e., "first pain" and "second pain." The so called "first pain" can be clearly localized, feels pricking, and occurs fast and first after nociceptive stimulation (Bromm and Treede, 1987a; Magerl et al., 1999; Beissner et al., 2010). First pain is considered to inform the individual about the location of an injury at and within the body and about the sensory quality of the injury. The so called "second pain" can less clearly be localized. Second pain is described as dull or pressing and occurs later after nociceptive stimulation than first pain (Price, 1988; Miltner, 1989; Magerl et al., 1999; Beissner et al., 2010). The prolonged second pain is considered to pull the individuals attention to the injury and to convey information to the brain that provides the basis for pain-related affect, arousal, and behavioral responses to limit further injury and to optimize recovery. Concerning the two types of pain, it has been shown that second pain is enhanced and first pain is suppressed when moderately painful heat is presented with a frequency of greater than 0.3 Hz (Price et al., 1977; Staud et al., 2007). When painful heat is presented with frequencies below 0.17 Hz, first pain is not suppressed and no enhancement of second pain occurs (Price et al., 1977; Staud et al., 2007). This is in line with Bromm's and Treede's suggestion (Bromm and Treede, 1987a) that second pain is perceived when first pain is reduced and vice versa.

In humans, different SSH have not been investigated to evoke different activation in brain stem structures so far. With the first (fMRI) experiment, we aimed at finding a method to test whether noxious heat stimulations with different SSH does activate brain stem structures differentially. According to animal studies, we expect a differential activation of PAG and RVM to different SSH. With the second experiment, we tested whether the different SSH used in the fMRI environment are associated with different pain qualities as stimulation with different SSH were associated with different behavioral responses in animal studies (Lumb, 2002).

## **MATERIALS AND METHODS**

We conducted two experiments, one inside and one outside the fMRI scanner. Both experiments used the same thermal stimulation with steep and shallow SSH. Subjects were informed about the procedure and provided written informed consent. To familiarize the participants with the experimental procedure and the stimulus types, each subject received a brief demonstration of the thermal stimulation prior to the experiment. Participants were otherwise naive about the purpose of the experiments. No subject had any history of neurological, psychiatric, or pain disorder. They were free to withdraw from the experiment at any time. The procedure was approved by the local ethics committee of the Friedrich Schiller University of Jena.

## **DETERMINATION OF THE PAIN SENSITIVITY**

Thermal stimuli were applied by a fMRI-compatible Peltier thermode (Medoc Advanced Medical systems; Ramat Yishai, Israel). The thermode had a surface area of 9 cm2 . Subjects were instructed to rate a series of thermal stimuli applied to the thenar eminence of their left hand using a modified Ellermeier scale (Ellermeier et al., 1991). This scale starts with 0 for "no pain" with an open scale with verbal description 1–10 ="just perceived," 11–20 ="clearly perceived but not painful," 21–30 ="very slightly painful," 31– 40 ="slightly painful," 41–50 ="medium pain," 51–60 ="strong pain," 61–70 ="very strong pain." It was explained that pain should be rated with values higher than 70 if pain becomes worse. The original Ellermeier scale has good psychophysical properties

(Ellermeier et al., 1991). We included 20 additional steps at the lower end to represent non-painful perceptions. Subjects were instructed to make a judgment regarding the categories first and subsequently rate the pain intensity within the range defined by this category. The rating requested was just the discrete number of rating that was monitored for further analysis.

To determine the pain sensitivity, thermode temperature was increased to a maximal stimulation temperature of 44°C–49°C in steps of 1°C. The procedure was as follows: a starting temperature of 34°C was established. Then, one full ramp with rise and fall of 10°C (2.5°C/s) was applied providing a maximal temperature of 44°C. Subject's intensity rating was recorded. The next starting temperature with a step of 1°C was established (up to a maximum of 39°C), followed by the next ramp with similar parameters (increase of 10°C with 2.5°C/s rise and fall). The procedure was finished at either 49°C maximal temperature or before maximal temperature of 49°C when subjects reported a rating of 51 or higher on the scale described above. This procedure allows the fitting of a stimulus-response curve presenting subjective ratings in dependence of the maximal temperature used for stimulation. It also familiarized the subjects with the kind of stimulation of the main experiments. For the succeeding main experiments, a Thot was determined as the temperature where the subject reported a value of 50 on our modified Ellermeier scale. Thot of all participants varied between 46.5°C and 49°C. Another maximal temperature of stimulation was used (Twarm = 40°C) providing ratings in the range below 20 on our scale.

## **THERMAL STIMULATION**

Trains of thermal stimuli with different SSH are used for the experiments. This is an ecologically valid procedure to induce pain percepts. Thermal stimuli were applied to the thenar eminence of the right hand. Subject received two different types of heat pulse trains (steep vs. shallow SSH) applied with two different temperature levels (Thot vs. Twarm). Heat pulse trains were balanced to control for order effects.

A design with four conditions was used, steep SSH with Twarm, steep SSH with Thot, shallow SSH with Twarm, and shallow SSH with Thot (**Figure 1**). The four conditions were presented in stimulation blocks. Five stimulation blocks containing one of each condition were presented throughout the whole experiment (**Figure 1**). Within a stimulation block, a baseline of at least 20 s (see below) was introduced between conditions (**Figure 1**). Each condition consisted of five heating ramps of identical type. During *WARM* conditions, temperature rose to 40°C (Twarm), whereas in the *HOT* conditions temperatures rose to Thot. The baseline temperature before stimulation was set to 10°C below Twarm/Thot and rose to these target temperatures with two different slopes: *steep SSH* runs had a slope of 7.5°C/s and *shallow SSH* runs had a slope of 2.5°C/s (**Figure 1**). Thus, painful heat peaks of the steep SSH stimuli were applied with a frequency of 0.3 Hz, whereas painful heat peaks of the stimuli with shallow SSH were applied with a frequency of 0.17 Hz in the *HOT* conditions. There was a baseline interval between stimulation blocks of 30 s.

A 30 s time interval with a constant baseline temperature (10°C below Twarm in the *WARM* and 10°C below Thot in the *HOT* conditions) was introduced between the steep and the shallow SSH

conditions of each temperature. The baseline temperature rose over the course of 20 s from 10°C below Twarm to 10°C below Thot for a change in stimulation from*WARM* to a succeeding*HOT* condition, or decreased from 10°C below Thot to 10°C below Twarm for a change in the stimulation from *HOT* to a succeeding *WARM* condition, respectively (**Figure 1**). This temperature was kept for another 20 s before the next heating ramps began. The sample was split concerning the order of shallow and steep SSH heating to control for order effects of conditions.

was presented once during each SB. Painful heat peaks of the steep SSH

## **Experiment 1**

Experiment 1 investigated the activation of the brainstem by means of fMRI for the different SSHs.

Sixteen healthy, right-handed subjects (seven male, nine female, 19–28 years) volunteered in the fMRI experiment. Subjects were paid C12 for completing the experiment. Prior to the experiment, the stimulus-response function to thermal stimulation of the subjects was examined as outlined above.

## **FUNCTIONAL IMAGE ACQUISITION**

Scanning was performed with a 3T magnetic resonance scanner (Tim Trio, Siemens Medical Systems, Erlangen, Germany). The experiment started with a high-resolution T1-weighted scan of the brain (192 slices, TE = 5 ms, FOV: 256 mm × 256 mm, resolution: 1 mm × 1 mm × 1 mm) for anatomical referencing and visualization. A shimming procedure preceded the succeeding functional MR scanning. The first four volumes were discarded in order to improve field homogeneity. In the experimental fMRI run, 650 volumes were acquired using a T2∗ weighted echo-planar sequence (TE = 75 ms, TR = 1.8 s; FOV = 192 mm × 192 mm). Each volume comprised 24 slices (2 mm thickness and 2 mm × 2 mm in-plane resolution) (see **Figure A1** in Appendix) which were prescribed parallel to the brainstem. The FOV covered the *a priori*-defined region of interest which was centered around the PAG and enclosed the upper brainstem and the midbrain (**Figures 2B**, **A1** in Appendix).

## **fMRI PREPROCESSING**

Preprocessing and analysis of fMRI data was performed using BrainVoyagerQX 2.1 (Brain Innovation, Maastricht, Netherlands). Primarily, all volumes were realigned to the first volume in order to minimize effects of head movements on data analysis. Further data preprocessing comprised spatial (6 mm full-width halfmaximum isotropic Gaussian kernel) and temporal smoothing (high pass filter: 15 cycles per run; low pass filter: 2.8 s; linear trend removal). The anatomical and functional images were co-registered (**Figure A2** in Appendix) and normalized to the Talairach space (Talairach and Tournoux, 1988).

Experiment 2, shallow and steep SSH were applied only for the hot temperature.

## **fMRI STATISTICAL ANALYSIS**

Statistical analyses were performed by multiple linear regression of the signal time course at each voxel. The expected blood oxygenlevel-dependent (BOLD) signal change for each of the four conditions (predictors) was modeled by a canonical hemodynamic response function. A random-effects General Linear Model was used to identify associated brain activity in all acquired slices. To minimize false-positive results (Straube et al., 2008) we tested whether the detected clusters survived a correction for multiple comparisons. We used the approach as implemented in Brain Voyager (Goebel et al., 2006), which is based on a 3D extension of the randomization procedure described by Forman et al. (1995). This procedure is based on the estimate of the map's spatial smoothness and on an iterative procedure (Monte Carlo simulation) for estimating cluster-level false-positive rates. After 1000 iterations, the minimum cluster size threshold that yielded a cluster-level falsepositive rate of 5% was applied to the statistical maps. Clusters reported here survived this control of multiple comparisons. For subsequent visualization of activated brain regions, the location of significantly activated regions was assessed by superimposing the results from group analysis on an averaged brain.

As the intent of this study was to characterize the changes in blood oxygen level dependent (BOLD) response to the painful stimulation, we compared the two different SSHs in the *HOT* conditions. The brain stem clusters found for this contrast constitute the basis for the analysis of further effects. The coordinates of the peak voxels were allocated to the anatomical structures with the assistance of an atlas of the human brain stem (Paxinos and Huang, 1995). For this comparison we also conducted repeated measures *t*-tests for the peak voxel of each structure.

Additionally, conditions we used repeated measures ANOVAs to assess the differential effects of the four conditions.

## **Experiment 2**

Experiment 2 was conducted to prove whether the different SSH are able to elicit different pain percepts. This is an important question with respect to the discussion concerning different types of afferents possibly involved during this type of stimulation.

Prior to the experiment, the stimulus-response function of the subjects was assessed analogously to Experiment 1. Stimuli with different SSH were then applied similarly to Experiment 1 with two exceptions: first, there was no *WARM* condition included. Second, the *HOT* condition was presented 10 times (5 times with steep and 5 times with shallow SSH). Directly after each *HOT* condition, subjects were requested to indicate how the stimuli of different SSHs were perceived. We used a (restricted) three-item verbal descriptor list which has been shown that it can reliably indicate whether the pain sensation evoked by the physical stimulus is the result of predominantly Aδ (first pain) or C-fiber activity (second pain) (Beissner et al., 2010). According to Beissner et al. (2010),"pricking" is an indicator for the first pain, while "pressing" or "dull" are indicators for the second pain. Thus, subjects were requested to choose the appropriate perception(s) from this list of three adjectives for the previous stimulation.

Twenty-three healthy right-handed subjects (3 male, 20 female, 19–28 years) volunteered in Experiment 2. Subjects were paid C5 for completing the experiment.

For the analysis of the data of Experiment 2, odds ratios (OR) were calculated separately for each of the three descriptors according to Beissner et al. (2010) as (A·D)/(B·C). The capital letters have the following meaning:

If "pricking" will be chosen more often for shallow SSH (OR < 1), then we might conclude that this stimulation preferentially activates Aδ-fibers. Accordingly, if "pressing" and/or "dull" will be chosen more often for steep SSH (OR > 1), then we might conclude that this stimulation preferentially activates C-fibers. Ninety-five percent confidence intervals, calculated as OR ± 1.96·(1/A + 1/B + 1/C + 1/D)0.5, were utilized to evaluate the significance of the respective ORs.

## **RESULTS**

## **EXPERIMENT 1**

First, we tested whether brain stem regions show activation during the noxious thermal stimulation compared to baseline. We found activation to both SSH in the painful *HOT* conditions compared to baseline in an inferior part of the PAG, probably including nucleus Raphe dorsalis (NRD) according to (Paxinos and Huang, 1995) [*t*(15) = 3.354, *p* < 0.005, *x*, *y*, *z*: 0, −29, −20] (**Figure 2A**, No. 3). More specifically, stimulation with shallow SSHs (vs. baseline) led to higher activation in a more superior part of the PAG [*t*(15) = 3.24, *p* < 0.01, *x*, *y*, *z*: −3, −28, −6, **Figure 2B**], whereas the stimulation with the steep SSHs yielded higher activation in the PAG/NRD complex [*t*(15) = 3.77, *p* < 0.005, *x*, *y*, *z*: 1, −25, −19, **Figure 2C**]. Furthermore, the steep SSHs in the Thot condition showed activation in a brain stem cluster that probably represents the RVM according to (Paxinos and Huang, 1995) [*t*(15) = 3.12, *p* < 0.01, *x*, *y*, *z*: 3, −33,−43, **Figure 2C**]. More importantly, we investigated whether these brain stem regions show differential responses under the two painful experimental conditions (*HOT*), i.e., with steep vs. shallow SSHs. This contrast revealed significantly stronger activation in the inferior part of the PAG [*t*(15) = 3.54, *p* < 0.005, 16 voxel, *x*, *y*,*z*: 1, 31, 17], NRD [*t*(15) = 4.93, *p* < 0.001, 40 voxel, *x*, *y*, *z*: 2, 25, 18], and RVM [*t*(15) = 3.82, *p* < 0.005, 16 voxel, *x*, *y*, *z*: 2, 30, 42]. No significant differences were found for the comparison between steep and shallow SSH in the*WARM* conditions (PAG [*t*(15) = 0.18, *p* > 0.1], NRD [*t*(15) = 1.57, *p* > 0.1], and RVM [*t*(15) = 0.53, *p* > 0.1]). The clusters of higher activation to steep vs. shallow SSHs are shown in **Figure 2D**; β-values of BOLD responses are shown in **Figure 3**. In contrast to the higher activation observed for steep SSHs, we did not find any statistically significant difference for the comparison shallow SSHs vs. steep SSHs for Thot.

Repeated measures ANOVAs for main effects and interactions of all other conditions were performed for β-values of peak voxels

where significant effects were found for the contrast of interest, i.e., steep vs. shallow SSH in the *HOT* condition. ANOVA revealed a main effect of temperature for the NRD cluster with stronger activation for Thot, a main effect of type of SSH in PAG with stronger activation for steep SSH, and an interaction for temperature × SSH in the NRD (**Table 1**). The significant interaction term for the NRD was investigated with a contrast analysis. There is a higher activation for steep SSH vs. shallow SSH in the *HOT* condition [*t*(15) = 4.932, *p* < 0.001], but no significant difference in activation between steep and shallow SSH in the *WARM* condition. Conversely, we found a higher activation for steep SSH with Thot as compared to with Twarm [*t*(15) = 4,058, *p* = 0.001] whereas no difference was found for the shallow SSH between Thot and Twarm.

Second, we used our restricted field of view (remember that the acquired slices concentrated on brainstem activation and did not allow the mapping of the whole "neuromatrix of pain" (Treede et al., 1999; Tracey and Mantyh, 2007; Iannetti and Mouraux, 2010; Schweinhardt and Bushnell, 2010)) to prove our experimental manipulation. We found higher activations to both SSH in the painful Thot conditions compared to baseline in the posterior cingulate cortex [PCC, *t*(15) = 3.498, *p* < 0.005, 252 voxel, *x*, *y*, *z*: 1, −19, 35; **Figure 2A**, No.1], in left amygdala [*t*(15) = 3.44, *p* < 0.005, 312 voxel, *x*, *y*,*z*: 21,−10,−12; **Figure 2A**, No. 2], and in the medial thalamus [*t*(15) = 3.19, *p* < 0.01, 40 voxel, *x*, *y*, *z*: −1, −22, 1; **Figure 2**]. These results indicate that our paradigm with different SSH is suitable to activate brain regions that have been found to process noxious thermal stimuli in other studies (Peyron et al., 2000).

## **EXPERIMENT 2**

Odds ratios (and 95% confidence intervals) were calculated separately for the three descriptors (**Figure 4**). Clearly, the descriptor "pricking" which is associated with first pain was chosen significantly more often for the stimulation with the shallow SSH, whereas the descriptor "dull" which is associated with second pain, was chosen significantly more often for the stimulus with the steep SSH. The descriptions for "pressing" did not reach a significant discrimination between the different SSH (**Figure 4**).

## **DISCUSSION**

Using different SSH, our primary finding is a stronger BOLD activity in response to trains of painful heat stimuli with steep SSH as compared to trains of painful heat stimuli with shallow SSH. Higher activation was found in the inferior part of the PAG, probably including the NRD according to (Paxinos and Huang, 1995), and a cluster that probably represents the RVM according to (Paxinos and Huang, 1995). We did not find any stronger activation in the brain stem for the contrast shallow vs. steep SSH. Thus, this

**Table 1 | Main effects and interactions of parameter estimates for the factorsTemperature and SSH within the brain stem.**

method is able to differentially activate structures in the human brain stem. To our knowledge, this is the first time that differential activation to peripheral noxious stimulation is demonstrated in humans.

The stronger activation in the PAG/NRD complex and the RVM observed to steep vs. shallow SSH is surprising with respect to animal studies. In animal experiments, the shallow SSH, and not the steep SSH yielded more activation. Several reasons might account for this result. First, the maximal temperatures of stimulation differed between the present experiment (49°C) and the animal studies [e.g., 55°C (Lumb, 2002)]. It is well known that the characteristics of nociceptors differ in their response behavior for this temperature range (Behbehani, 1995). Second, the differences might be due to our spatial resolution. The resolution of brain scans with 2 mm × 2 mm × 2 mm in the present study might not have been sufficiently high to detect activations in different columns of the PAG. Third and probably most important, single ramps with shallow SSH were employed in the animal experiments.  In contrast to the animal experiments, a train of five succeeding heating ramps without gaps was employed in our study. Thereby, the steep SSH resulted in a frequency of heat peaks of 0.3 Hz. A stimulus frequency of 0.3 Hz is known to produce the phenomenon of temporal summation of second pain, i.e., TSSP (Price et al., 1977; Herrero et al., 2000). TSSP is considered to result from C-fiber evoked responses in dorsal horn neurons, termed"windup"(Herrero et al., 2000; Sarlani and Greenspan, 2005). Thus, the involved fibers activated by the steep SSH in animal studies are Aδ-fibers while the activation with repeated steep SSH in our experiment might preferentially involve C-fibers. Following this interpretation, the activation of the brain stem by steep SSH of stimulus trains in our study has to be compared to the single shallow SSH in the animal experiments. Considering this, both experiments yielded similar results.

The results of Experiment 2 are of crucial importance for the latter consideration. It investigated the quality of pain percepts that are elicited by different SSHs in humans. Indeed, we

found evidence that the steep SSH was associated with the percept "dull" whereas the shallow SSH was associated with the percept "pricking."In accordance with Beissner et al. (2010), these descriptors distinguish best between first pain and second pain. Thus, steep SSH might be associated with a predominant activation of C-fibers while the shallow SSH might be associated with a predominant activation of Aδ-fibers. This interpretation is in line with other studies using trains of painful thermal stimuli that were associated with different pain percepts (Price et al., 1977; Staud et al., 2007). Characteristics of first pain were associated with a frequency of 0.17 Hz between the painful heat peaks whereas characteristics of second pain were associated with a frequency of 0.3 Hz between the painful heat peaks (Price et al., 1977; Staud et al., 2007). However, it should be mentioned that heat is perceived as painful in humans only at temperatures above 43°C (Julius and Basbaum, 2001). Thus, just the top of each heating ramp can be considered as a painful heat peak. In the present study, the painful heat peaks of the steep SSH stimuli were applied with a frequency of 0.3 Hz whereas heat peaks of the stimuli with shallow SSH were applied with a frequency of 0.17 Hz so that these ramps fulfill the frequency criterion for painful stimulation. Taking together the two experiments, we suggest that the different SSH probably activate different types of peripheral input resulting in different pain percepts (first vs. second pain).

We found activation of the PAG and the RVM to stimulation with noxious heat to both types of SSH. Several brain imaging studies have observed activations in brainstem structures to nociceptive stimulation (Apkarian et al., 2005; Tracey and Mantyh, 2007; Eippert et al., 2009; Schweinhardt and Bushnell, 2010). Brainstem modulation of neuronal activity in the spinal cord has been reported since more than a century ago (Bernard, 1858) and is thought being involved in top-down control of pain. In particular, the midline PAG integrates input from the spinal cord, cerebral cortex, and numerous other brainstem nuclei (Apkarian et al., 2005; Eippert et al., 2009). Its stimulation in humans was shown to result in antinociception and analgesia (Hosobuchi et al., 1977). In contrast, its lesion might result in chronic pain (Basbaum and Fields, 1984). Our result of an increased activity in the inferior part of the PAG in response to steep as compared to shallow SSH might, therefore, be an important result. Based on the animal studies mentioned above, this higher activation might indicate that the steep SSH stimulation, but not the shallow SSH stimulation might trigger the nocifensive part of PAG. Moreover, the hypothesis that this type of activation might result in an activation of nocifensive reaction might be tested and, if true, possibly be used for chronic pain patients.

We found a higher activation within the inferior part of the PAG to stimulation with steep as compared to shallow SSH. This result is in line with previous studies in animals. Animal studies have found different activation patterns within distinct columns of the PAG to preferential C- and Aδ-fiber stimulation (Lumb, 2002; Lumb et al., 2002; Parry et al., 2008). It has also been proposed that the different columns of the PAG not only differ with respect to the influence of the ongoing nociceptive information processing as outlined above, but also mediate different coping strategies (Lumb et al., 2002). Preferential Aδ-fiber stimulation is

associated with activation of dorsolateral and lateral columns of the PAG which in turn result in active coping strategies (Lumb et al., 2002). The activation of these columns evokes sympathetic excitation. Passive coping strategies are closely linked to the ventrolateral columns of the PAG activated by preferential C-fiber stimulation. The activity of the ventrolateral PAG is associated with sympathoinhibition (Lumb et al., 2002). Correspondingly, the PAG mediates differential control of spinal nociception as part of a defensive response or as withdrawal. It has been established to carry out integrative functions for cardiovascular and respiratory regulation, for sensory modulation, and for different motor behaviors (Clement et al., 2000; Morgan and Carrive, 2001; Subramanian et al., 2008; Heinricher et al., 2009). Reflecting these results on our data, the findings of this study provide a hint that stimulation with steep SSH but not shallow SSH might engage nocifensive mechanisms as shallow SSH do not seem to influence the PAG in the same way as steep SSH.

Although our main focus in this study was the PAG, we found activations in two neighboring and functionally related areas, i.e., in NRD and RVM. The NRD is embedded in the ventromedial part of the PAG (Mantyh, 1982) and was shown to modulate responses caused by noxious stimulation of the spinal dorsal horn neurons by its descending projections (Yu et al., 1988). In addition, PAG and NRD project to the spinal cord indirectly via the RVM, which is situated centrally around the pontomedullary junction. It includes the NRM and the adjacent reticular formation. It is known to project diffusely to dorsal horn laminae, including superficial layers and deep dorsal horn structures (Fields and Heinricher, 1985). Similar to the PAG, the RVM has a dual role in pain control: it is as well able to inhibit and to facilitate nociceptive input and can thus be considered as the output of the midline pain-modulation system. Profound analgesia can be produced by stimulating the NRM which is due to a decrease in responsiveness of spinothalamic dorsal horn neurons to input from peripheral nociceptors (Besson and Chaouch, 1987). Alternatively, analgesia evoked by stimulation of the ventral sites of the PAG can be blocked by lesion of the RVM (Behbehani and Fields, 1979; Prieto et al., 1983). In the light of these data, activations found in the present study might mirror the descending pathway to the dorsal horn. These observations parallel the results of the present study that stronger activation in RVM can be observed in response to steep SSH.

We found stronger activation in response to steep SSH stimulation both in parts of the brain stem as well as in some other structures in the field of view, i.e., the PCC and the amygdala. As argued above, we suggest that stimulation with steep SSH might preferentially activate C-fiber input. In this sense, the fMRI results of our study are in line with previous studies that also found stronger activation to selective C-fiber stimulation compared to Aδ-fiber stimulation. Stronger activations have been reported in structures associated with the affective processing of nociceptive information (i.e., ACC (Qiu et al., 2006); anterior insula (Weiss et al., 2008)). Similar to the present study these authors (Qiu et al., 2006; Weiss et al., 2008) also did not find any stronger activation to selective Aδ-fiber stimulation when comparing it with selective C-fiber stimulation. This might be another hint to the correctness of our suggestion concerning preferential C-fiber activation by steep SSH stimulation.

Several limitations of our study have to be considered that might influence future research. First, the conditions steep and shallow SSH were determined by the different slopes of heating. Different slopes affect the frequency of the painful heat peaks that are essential for the paradigm. There are two possibilities to proceed further, with the same number of stimuli within a condition (i.e., 5× up and down) or the same duration within a condition, but a different number of stimuli within a condition. We decided to use the same number of stimuli to have the same number of painful events within a condition. However, this leads to different durations of stimulation between the two SSHs used. Future studies might explore the effects of the two types of SSHs using the same length but unequal number of painful events within a condition.

Second, the energy transmitted to the skin depends on the frequency and duration of stimulation within a condition. In our study, the transmitted energy (area under the curve) was higher in the shallow SSH condition. Future studies might utilize the same amount of transmitted energy. However, different slopes of heating will then request different durations of baseline between ramps. These segments in turn might rise additional percepts in difference to the heat stimulation that might influence the results. However, it should be mentioned that heat pain receptors start firing at about 43°C (Julius and Basbaum, 2001) so that it is quite difficult to produce ramps that have the same amount of energy in the painful range; moreover, it seems to be impossible to produce ramps with different SSH that have the same energy both in the noxious as well as in the innocuous temperature range. Taking this consideration into account, the difference in transferred energy above the temperature threshold of 43°C is smaller as compared to differences in transferred energy for the whole heating ramps.

Third (as mentioned earlier),the resolution of 2 mm × 2 mm × 2 mm might not be sufficient to detect further differentiation within the brain stem,especially within the PAG. Possibly,scanning with higher field strengths might identify a columnar organization of the human PAG.

In summary, we found stronger activation in the inferior part of the PAG and in the RVM in response to painful stimulation with steep SSH. These observations provide first evidence for selective activation of the midbrain structures PAG and RVM in the human brainstem by different SSH. Therefore, this stimulation can be used when human brainstem structures are in the focus of interest during nociception. The specific activation of the midbrain to steep SSH seems to be associated with the specific perception of second pain and might possibly be related to passive coping strategies.

## **ACKNOWLEDGMENTS**

Thanks to Maria Carl for her contribution in data acquisition and analyses. This study was funded by the German Federal Ministry of Education and Research (BMBF, No.: 01EC1003).

## **REFERENCES**

J., et al. (2009). Activation of the opioidergic descending pain control system underlies placebo analgesia. *Neuron* 63, 533–543. doi:10.1016/j. neuron.2009.07.014

the nucleus raphe magnus. *Anesth. Analg.* 98, 414–419. doi:10.1213/01. ANE.0000094334.12027.06

nociceptive response of dorsal horn neurons and efferent pathway analysis in rats]. *Sheng Li Xue Bao* 40, 231–239.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 31 May 2013; accepted: 15 August 2013; published online: 06 September 2013.*

*Citation: Ritter A, Franz M, Dietrich C, Miltner WHR and Weiss* *T (2013) Human brain stem structures respond differentially to noxious heat. Front. Hum. Neurosci. 7:530. doi: 10.3389/fnhum.2013.00530*

*This article was submitted to the journal Frontiers in Human Neuroscience.*

*Copyright © 2013 Ritter, Franz, Dietrich, Miltner and Weiss. This is an open-access article distributed under the terms of the* *Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## **APPENDIX**

**FIGURE A2 | Section (between green lines in mid – saggital view) of the functional mean BOLD – data coregistered to the anatomical scan (neurological convention left** = **left).**

## Tonotopic organization in the depth of human inferior colliculus

## **David Ress 1\* and Bharath Chandrasekaran2**

1 Neuroscience Department, Center for Perceptual Systems, Imaging Research Center, Institute for Neuroscience, The University of Texas at Austin, Austin TX, USA 2 Department of Communication Sciences and Disorders, Center for Perceptual Systems, Institute for Neuroscience, The University of Texas at Austin, Austin, TX, USA

## **Edited by:**

Florian Beissner, Martinos Center for Biomedical Imaging, USA

## **Reviewed by:**

Elia Formisano, Maastricht University, Netherlands Marc Schönwiesner, University of Montreal, Canada

## **\*Correspondence:**

David Ress, Neuroscience Department, Baylor College of Medicine, 1 Baylor Plaza, S104N MS BCM240, Houston, TX 77030, USA e-mail<EMAIL_ADDRESS>
Experiments in animal models indicate that inferior colliculus (IC), the primary auditory midbrain structure, represents sound frequency in a particular spatial organization, a tonotopy, that proceeds from dorsal and superficial to ventral and deeper tissue. Experiments are presented that use high-resolution, sparse-sampling functional magnetic resonance imaging (fMRI) at 3T to determine if tonotopic gradients can be reliably measured in human IC using high-resolution fMRI. Stimuli were sequences of bandpass-filtered noise with different center frequencies, presented sequentially while fMRI data were collected. Four subjects performed an adaptive frequency-discrimination task throughout the experiment. Results show statistically significant tonotopic gradients within both ICs of all subjects. Frequency gradients as a function of depth were measured using surface-based analysis methods that make virtual penetrations into the IC tissue. This organization was evident over substantial portions of the IC, at locations that are consistent with the expected location of the central nucleus of IC.The results confirm a laminar tonotopy in the human IC at 3T, but with a heterogeneous, patchy character. The success of these surface-based analysis methods will enable more detailed non-invasive explorations of the functional architecture of other subcortical human auditory structures that have complex, laminar organization.

**Keywords: audition, tonotopy, inferior colliculus, fMRI, midbrain**

## **INTRODUCTION**

Organization by frequency, *tonotopy*, is one of the central organizing principles of the auditory system. The inferior colliculus (IC) is the primary auditory midbrain structure that forms a major convergence zone within the auditory system, receiving obligatory bottom-up connections from the brainstem auditory structures, and direct corticofugal feedback connections from auditory cortex. The IC consists of three prominent and functionally distinct subdivisions: the tonotopically organized central nucleus, the multisensory external cortex, and the dorsal cortex, which receives corticocollicular connections, primarily from the auditory cortex. Invasive studies on animal models show that the CN of IC exhibits tonotopic organization along the dorsal-to-ventral direction, which roughly corresponds to the depth direction within this laminated tissue (Schreiner and Langner, 1997; Malmierca et al., 2008; Baumann et al., 2011; Cheung et al., 2012). Electrophysiological studies in animal models show finely tuned topographical mapping of frequency along the laminar axis of the IC (Schreiner and Langner, 1997; Baumann et al., 2010). Similarly, anterograde labeling studies from the cochlear nucleus also show evidence for laminar frequency selectivity, but these results suggest a more heterogeneous, patchy distribution of frequency selectivity for the bottom-up inputs to IC (Oliver and Morest, 1984; Hofstetter and Ehret, 1992; Loftus et al., 2008).

In humans, the small size of the IC (<9 mm), coupled with the limited resolution of conventional human neuroimaging studies, have precluded the non-invasive examination of frequency representation. Further, subcortical structures are harder to image because of their susceptibility to physiological noise (Guimaraes et al., 1998). Conventional functional neuroimaging studies employ 3–10 mm sampling and spatial blurring, making it impossible to resolve functional subdivisions within the IC. In addition, the IC has a complex organization that represents multiple auditory features, with tonotopic specificity shown to be largely resolved along the laminar axis (Zook et al., 1985; Malmierca et al., 1993; Schreiner and Langner, 1997). The complicated internal structure of the IC requires a more detailed examination using high-resolution imaging and depth-resolved analyses.

Two recent studies have examined tonotopic organization using a resolution (2-mm FWHM) that is higher than traditional neuroimaging studies. Baumann et al. (2011) examined tonotopic and rate-dependent gradients in rhesus monkeys using a 4.7 T MRI scanner and found an orthogonal relationship between frequency and rate representation in the IC. More recently, De Martino et al. (2013) examined tonotopic organization in the human IC using a 7 T MRI scanner. Both studies showed a dorsolateral-toventromedial tonotopic gradient in the IC preferentially responding to low-to-high frequencies, respectively. While these studies demonstrate the utility of very high-field strengths in resolving tonotopic gradients within the IC, it is unclear whether clear tonotopy is evident at smaller field strengths. Functional neuroimaging studies at 3 T have demonstrated tonotopy in the human auditory cortex, in a manner consistent with the organization found in other primates (Barton et al., 2012; Moerel et al., 2012), suggesting the possibility that tonotopy can also be evidenced at 3 T.

The goal of this study is to examine frequency representation in the human IC using high-resolution (1.2 mm) functional magnetic resonance imaging (fMRI) at 3 T, adapted for auditory stimulation. During fMRI, subjects listened to sequences of five sounds in the range of 0.25–8 kHz, while performing a frequencydiscrimination task (**Figure 1**). Surface-based analyses were then used to systematically evaluate tonotopy within human IC.

## **MATERIALS AND METHODS**

## **SUBJECTS**

Subjects were four male adults (age range 21–31), with normal hearing (air-conduction thresholds <20 dB nHL for octaves from 0.25 to 8 kHz). Informed consent was obtained from all subjects in accordance with procedures approved by the Institutional Review Board of The University of Texas at Austin.

## **STIMULUS AND TASK**

Auditory stimuli, presented through MRI compatible earphones, consisted of bandpass-filtered (0.6-octave bandwidth) white noise, presented with five different center frequencies, logarithmically spaced on 0.25–8 kHz (**Figure 1A**). The sounds were followed by a sixth quiet-period with no sound presented to act as a reference. This sequence was repeated five times during each 6-min run.

To ensure that participants paid consistent attention to the stimuli, they performed a challenging two-interval forced-choice frequency-discrimination task. Sounds were presented in pairs, each lasting 2.85 s with a 0.3-s gap. The sound in one of the two periods, randomly chosen, had a slightly higher center frequency than the other, and subjects used a MRI compatible button box after each sound pair to indicate their judgment as to which interval had the higher frequency. No response was given during the final quiet-period of each sequence. During stimulus periods, the frequency difference, typically 3–8%, between the two sounds was adjusted in real time using a pair of randomly interleaved two-up, one-down staircases. After each run, the subject's current frequency-discrimination threshold value was estimated as the average frequency difference over the latter half of the trials in that run. Frequency difference for the next run was then initialized to a value 30% greater than this estimate. This creates a stable pattern of performance in which subjects perform very well on the task at the start of each run, and then steadily return to their threshold value (**Figure 1B**). In addition, subjects practiced the task in advance of scanning so that their performance was stable during the fMRI. Across runs, subjects' performance was 75–85% correct.

## **fMRI PROTOCOL**

Stimulus presentations occurred every 12 s while high-resolution (1.2-mm cubic voxels) T2\*-weighted images were collected in a sparse-sampling approach. In sparse-sampling designs, one attempts to sample the hyperoxic peak of the BOLD response evoked by preceding stimulus, then allocate a sufficient delay to allow that response to decay away before the next sample is obtained. Moreover, the delay is adjusted so that the nuisance response from the scanner-acquisition noise is minimized (Hall et al., 1999). The purpose of the long presentation period is to drive the desired stimulus-evoked response into a fairly stable"flattop"so that it can be sampled by the fairly long 3-s-duration 3-shot acquisition that immediately follows each stimulation period. Our timing choice was guided by previous sparse-sampling approaches that have examined auditory processing in cortical and subcortical structures (e.g., Gaab et al., 2003; Wong et al., 2008; Chandrasekaran et al., 2012).

Imaging was performed on a 3 T scanner (GE Signa Excite HD) using an eight-channel head coil. Ten 1.2-mm-thick quasi-coronal slices were acquired sequentially over a 170-mm field-of-view during a 3-s period immediately following each sound pair. A set of T1-weighted structural images was obtained on the same prescription using a three-dimensional (3D) RF-spoiled GRASS (SPGR) sequence (15°flip angle, 0.78-mm pixels). These images produced good gray-white tissue contrast and were used to align the functional data to a segmented structural reference volume, described below.

Functional image acquisition used a 3-shot spiral trajectory (Glover and Lai, 1998) with TE = 40 ms and TR = 1 s that has provided good results in human midbrain (Katyal et al., 2010, 2012). Between acquisitions, T1 equilibrium was maintained by continual slice excitation, but acquisition and crusher gradients were suppressed to create quiet-periods for presentation of the audio stimuli. Noise was further reduced by setting slew limits for the excitation pulses to 30 T/m/s. The scheme resulted in a noise reduction of ∼15 dB during the quiet-periods. A similar approach was used previously with EPI acquisitions (Schwarzbauer et al., 2006).

The use of the three-shot spiral acquisition acts as a temporal anti-aliasing filter that substantially reduces high-frequency physiological nuisance effects, particularly cardiac pulsations. In single-shot acquisition, each fMRI slice is acquired very quickly (∼20–30 ms) compared to hemodynamic time scales, but the samples are widely spaced (e.g., 2–3 s). This results in aliasing of high-frequency fluctuations driven by cardiac pulse (>1 Hz) and respiration (∼0.25 Hz). In our multi-shot imaging, acquisition is now distributed over multiple brief sampling periods separated by the relatively short TR = 1 s. This preparation implements a low-pass filter upon the data. This avoided the use of cardiac gating, and also reduced respiratory artifacts (Ress et al., 2007; Katyal et al., 2010, 2012).

The same monotonically increasing sequence of five sound frequencies and quiet-period was repeated five times in each fMRI run (6-min duration, 120 fMRI volumes per run), and 10–12 runs were collected in each scanning session. The fMRI data were compensated for subject motion using a robust intensitybased scheme (Nestares and Heeger, 2000). Data were also corrected for spatial variations, introduced mostly by receiver-coil inhomogeneity, using a spatially smoothed homomorphic normalization to the temporal mean of the fMRI time series, and high-pass filtered to reduce the effects of slow image-intensity drifts. The anatomical images collected in each session were then used to align and transform the functional data to a structural 3D reference volume, which was acquired for each subject in a separate session (**Figure 2A**). The structural reference volume was T1-weighted with good gray-white contrast, acquired using a 3D, inversion-prepared, SPGR sequence (*T*I = 450 ms, 0.7-mm isometric voxel size).

The noise in fMRI data is known to have a non-Gaussian distribution (Biswal et al., 1995; Holmes et al., 1997; Glover et al.,

**FIGURE 2 | Surface-model construction**. **(A)** MRI volumes obtained with 0.7-mm voxels were segmented to delineate the boundaries of brainstem tissue (blue). **(B)** A smooth surface was constructed at the segmentation boundary; gray overlay shows curvature. **(C)** Enlargement showing IC boundaries (black).

2000; Kruger and Glover, 2001). Typically, fMRI data is therefore "whitened" by performing a mixture of spatial and temporal blurring (Worsley and Friston, 1995; Friston et al., 2000). To maintain high spatial resolution, we wanted to avoid any blurring. Therefore, we used "bootstrapping," a non-parametric method to estimate statistical significance. The bootstrapping procedure obtains a high-quality estimate of the statistical distribution of the observed data. The data, for example, the BOLD response from 55 presentations of a particular stimulus frequency, are resampled with repetition, and then averaged together. Each iteration of this procedure thus selects both a random subset of the original data, and a random set of weight factors. It can be shown analytically (Ephron and Tibshirani, 1998), that many repetitions of this procedure converges very nearly to the original distribution of the data. Thus, bootstrapping entirely avoids the need for assumptions of a particular, usually Gaussian, noise distribution. Instead, the actual distribution is estimated, and this distribution is used to directly calculate confidence intervals upon the data. For our data here, we resampled the data with replacement 1000 times; error bars then show the 68% confidence intervals on the bootstrapped distribution, which is similar to what would be presented for the standard-error-of-the-mean for a normal distribution.

## **IMAGE ANALYSIS**

We segmented the brainstem tissue (**Figure 2A**) using the ITK-SNAP application (Yushkevich et al., 2006). A smooth surface was then interpolated at cerebrospinal fluid-tissue interface of the IC (**Figure 2B**) using a deformable-surface algorithm (Xu et al.,2006). This surface provided vertices and outward normal vectors used as a reference for the depth calculations described below as well as a means to visualize the functional data (**Figure 2C**). A nearestneighbor Euclidean distance map was calculated between the IC tissue voxels and the vertices of the refined IC surface. Thus, each volume voxel was associated with a frequency response across the five tones, and a depth coordinate (Ress et al., 2007; Khan et al., 2011).

As a first form of analysis, we calculated the BOLD response amplitude differences between the highest and lowest stimulus frequencies (**Figure 3D**). This procedure is analogous to that used by Baumann et al. (2011).