# **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: researchtopics@frontiersin.org

# **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 <sup>3</sup>**

<sup>1</sup> Department of Neuroradiology, Somatosensory and Autonomic Therapy Research, Hannover Medical School, Hannover, Germany

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

<sup>3</sup> 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: coffeefellow@gmail.com

**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 articles have been divided into a functional (Brooks et al., 2013; Henderson and Macefield, 2013; Ress and Chandrasekaran, 2013; Ritter et al., 2013) and a structural section (Deistung et al., 2013; Ford et al., 2013; Lambert et al., 2013; Yeo et al., 2013; Singleton et al., 2014).

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 <sup>1</sup>\*, Olivia K. Faull <sup>2</sup> , Kyle T. S. Pattinson<sup>2</sup> and Mark Jenkinson<sup>2</sup>**

<sup>1</sup> Clinical Research and Imaging Centre, University of Bristol, Bristol, UK

<sup>2</sup> 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: jon.brooks@bristol.ac.uk

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 CO<sup>2</sup> 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*<sup>0</sup> 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*<sup>0</sup> 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*<sup>0</sup> 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<sup>∗</sup> . The image is created by either computing the quotient of the two images (per TR) or by calculating the apparent T2<sup>∗</sup> 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<sup>∗</sup> 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. Correlation between the seed region (e.g., brainstem) and the remainder of the brain is performed, and an arbitrary threshold used to define a mask of these areas (excluding the seed region). 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 CO<sup>2</sup> in expired breath (PETCO2), and include this measure as an explanatory variable (EV) at the analysis stage (Pattinson et al., 2009a,b). PETCO<sup>2</sup> varies tightly with PaCO<sup>2</sup> 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 PETCO<sup>2</sup> 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 PaCO<sup>2</sup> 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 k<sup>1</sup> and k<sup>2</sup> regressors, which can be from very different sets, yields the condition that log[RSS(*k*1)/RSS(*k*2)] < (*k*<sup>2</sup> − *k*1) × log(*N*)/*N*

## **REFERENCES**


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


when the model with k<sup>1</sup> 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: alexander.ritter@uni-jena.de

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 cm<sup>2</sup> . 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<sup>∗</sup> 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. Such stimulation was shown to preferentially activate C-fibers, whereas single steep SSH are able to preferentially excite Aδ-fibers (Lumb, 2002; Lumb et al., 2002). 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 <sup>1</sup>\* and Bharath Chandrasekaran<sup>2</sup>**

<sup>1</sup> Neuroscience Department, Center for Perceptual Systems, Imaging Research Center, Institute for Neuroscience, The University of Texas at Austin, Austin TX, USA <sup>2</sup> 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: ress@bcm.edu

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*<sup>I</sup> = 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).

As a second form of analysis, we characterized the frequency response for each voxel within the tissue of IC by calculating the centroid (first moment) of its frequency response (**Figure 4D**): *f*cent = Σ*<sup>i</sup> fiA*(*f* )/Σ*iA*(*fi*), where *A*(*fi*) is the fMRI response amplitude to sound at frequency *f<sup>i</sup>* , and *i* indexes across the five frequencies. Amplitudes for each sequence of five stimuli were made positive definite by subtracting the minimum value. Because the difference between sound pairs was small (typically <5%), the frequency increment was neglected in these calculations. Each sound was repeated typically 55 times in a session. Using a bootstrapping approach as described above, we estimated confidence intervals on the centroid frequencies. To estimate tuning width we used the second moment of the amplitude response: *f*tune = p Σ*i*(*f<sup>i</sup>* − *f*cent) <sup>2</sup>*A*(*fi*)/Σ*iA*(*fi*).

We had originally planned to use a phase-encoded analysis approach, similar to what has been used in the measurement of

**FIGURE 3 | Response differences between the highest and lowest frequency stimuli for subject 1**. **(A–C)** Response difference is shown as a color overlay upon axial high-resolution structural anatomy images of the IC. Transparency is modulated by statistical significance of the sound stimuli vs. quiet-periods. Insets show location of each slice on the IC surface model. **(D)** Responses as a function of depth to the lowest- (red) and highest-frequency stimuli (blue); difference is shown in black. These

responses were obtained from a single depth penetration with a strong tonotopic gradient shown by the black rectangle in **(A)**. Thin dashed lines show 68% confidence intervals (e.g., SEM). **(E)** Response difference vs. depth for all voxels in another penetration marked in **(B)**. Red line is a linear fit the data beginning at 1-mm depth into the tissue. **(F)** The slopes of these linear fits are shown as a color overly upon this surface rendering of the subject's IC.

blue near 3.5 kHz. **(E)** Centroid frequency for all voxels in the penetration marked by the black rectangle in **(B)**. A strong depth gradient that begins 1 mm below the surface is captured by the linear fit shown by the red line. **(F)** The linear fit slope data is marked as a color overlay upon the IC surface model.

visual retinotopy and tonotopy in cortex (Engel et al., 1997; Barton et al., 2012). However, the very slow period, 72-s, of each stimulus cycle made this analysis susceptible to low-frequency scanner drift noise. Although they were therefore noisier, the phase-encoding results were qualitatively similar to those obtained by the centroid-frequency approach described above.

We evaluated the gradients of both high-low response difference and centroid frequency with depth across the IC in a depth range of 0–5 mm. The notion that IC frequency preference is organized by depth is has been suggested by previous animal studies (e.g., Schreiner and Langner, 1997; Loftus et al., 2008), and the utility of such methods with high-resolution fMRI was demonstrated by our previous work (Katyal et al., 2010). The portion of the surface corresponding to IC, estimated based on positive tissue curvature, ranged from 110 to 182 mm<sup>2</sup> across the four subjects. At every vertex in this portion (328–433 locations) on each IC surface model, we created normally oriented, 1.2-mm-diam cylinders. Specifically, we start from each surface vertex and create a patch corresponding to all surface vertices with manifold distances <0.6 mm from the starting vertex. These surface vertices are then extended along the inward-directed mean surface normal of the patch to penetrate the full depth of the tissue with a cylinder of constant diameter. We use these cylinders to probe the variation of response difference (**Figure 3E**) or centroid frequency (**Figure 4E**) with depth in the tissue, a hemodynamic analog to making electrode penetrations into the tissue. We then recorded the centroid frequency as a function of depth along the virtual penetration, and fit the result with a straight line; the slope of the fitted line, *g*, was an estimate of tonotopic gradient with depth. Linear fits were attempted after removing a variable amount of superficial tissue that varied from 0 to 1.5 mm, and the starting depth that yielded the best linear fit was chosen. Significance of the slope values was again obtained by a bootstrapping procedure. The amplitude-difference or preferred-frequency data obtained from each of 1000 resampled averages was fit as a function of depth by a straight line. The fraction of the resampled data sets with negative slopes provides an estimate of the *p* value for the gradient.

The orientations of those penetrations that corresponded to significant tonotopic depth gradients were used to obtain a measure of the mean tonotopic depth gradient. Specifically, we formed a weighted average: h*n*i = Σ*igini*/Σ*ig<sup>i</sup>* , where *n<sup>i</sup>* is the direction of each significant penetration *i*, having depth gradient *g<sup>i</sup>* . The resulting unit vectors yield two angles, one in the sagittal plane, and the other in the axial plane. We adopt the convention that an angle of 0°points directly anterior, with positive sagittal angles rotating superior, and positive axial angles rotating to the right.

## **RESULTS**

The sparse-sampling scheme captured large amplitude changes that often exceeded 1% of the local mean. Response differences showed strong depth gradients in many portions of each IC that was studied (**Figure 3**). The gradients had a complex spatial distribution, but bands of low- and high-frequency responsive tissue could be seen on particular slices (**Figures 3A–C**). Virtual penetrations were made by constructing cylinders were normally oriented to the superficial surface of the tissue (shown schematically by black rectangles in **Figures 3A,B**). A thin band of overlaying tissue, consistent with the anatomic location of the external cortex, 0–1.5-mm thick, was observed in some regions, e.g., the lateral surface of the right colliculus in **Figure 3B**. Since tonotopy is not a feature of the external cortex of the IC, this region is ignored in the subsequent depth analysis.

Responses to the lowest and highest-frequency stimuli often showed opposing gradients with depth for many of the virtual penetrations, yielding a strong positive depth gradient in the high-low response difference (**Figure 3D**). The difference is initially noisy or constant as we penetrate through what we assume to be external cortex, then rises strongly. For each penetration, we fitted a straight line to the depth gradient after discarding the putative external cortex (**Figure 3E**). These virtual penetrations and linear fits were performed across the entire surface of both ICs in each subject. The slopes of the fits could then be visualized as color overlay upon the IC surface models (**Figure 3F**). Depth-gradient results for all subjects show regions of significant positive depth gradients across both ICs, with a patchy spatial organization very similar to that obtained from the frequency-centroid analysis presented below.

**Figures 4A–C** shows the spatial structure of centroid frequency on three axial slices in subject 1; examples from other slice orientations and subjects are shown in Figure S1 in Supplementary Material. Color overlay shows centroid frequency with the transparency modulated by the sound > quiet *p* value: voxels with *p* < 0.05 are opaque, logarithmically transitioning to transparent for voxels with *p* ≥ 0.5. Much like the amplitude-difference data, the preferred-frequency responses exhibit patchy bands of frequency selectivity.

Example voxel responses taken from a region with significant tonotopic depth gradients in centroid frequency illustrate their selectivity (**Figure 4D**). Superficial voxels often displayed low-frequency preference, with a distinct low-pass response characteristic. Similarly, deep voxels often had high-pass characteristics. Voxels at intermediate depths had mid-frequency preferences overall, but had more heterogeneous response patterns that often included a mixture of low- and high-frequency response.

The centroid frequency, *f*cent, was quite stable in all IC voxels. Bootstrapped confidence intervals had a mean value <0.5 octaves, and did not exceed 0.9 octaves. Thus, compared to the 5-octave frequency band of the stimuli, the experimental variability in centroid frequency was small, typically <±10%.

In portions of IC, the centroid frequency showed positive gradients with depth. Examples of virtual "penetrations" into the tissue clearly show the depth gradients (**Figure 4E**). The linear fit captured the depth gradient magnitude, which was then displayed upon the IC surfaces models (**Figure 4F**).

Significant (*p* < 0.05) patches of tonotopic representation with depth were evident in all subjects (**Figure 5**). Significant gradients varied from 0.3 to 1.2 octaves/mm. The transparency of the depthgradient color overlay is modulated by its statistical significance: the overlay is opaque for *p* < 0.05, and becomes more transparent as *p* approaches 0.5. The spatial distribution of the tonotopic gradients was complex and varied substantially amongst the four subjects. Gradients were typically stronger and more extensive on the left colliculus then on the right. Left-colliculus gradients were typically larger on the lateral and rostral surfaces, with less distinct gradients on the medial surface. Gradients were generally least significant over the very center of the colliculus. Weak but significant regions of reversed depth gradients were present around the edges of the ICs of subjects 1 and 2. The centroid-frequency depth gradients were somewhat stronger than those obtained by the high-low response differences described above, but both exhibited qualitatively similar spatial distribution across the IC surfaces.

Various features of the centroid-frequency data are summarized in **Table 1**. Total collicular area varied from 110 to 182 mm<sup>2</sup> . Areas of significant depth gradients varied from 19 to 32 mm<sup>2</sup> (mean across subjects of 25 mm<sup>2</sup> ), which was 18% of the total area of bilateral IC. About twice that area (53 mm<sup>2</sup> ) showed trends (*p* < 0.16) toward depth gradients. Mean gradient was 0.53 octaves/mm, about half of what would be required to linearly resolve the full stimulus gamut in a 5-mm depth of tissue (see Discussion). Linear fits within each penetration were typically of high-quality, *R* <sup>2</sup> ∼ 0.5. Mean tonotopic depth gradients for 3/4 subjects were inclined downward; subject 2 had a mean gradient pointing slightly upward. Three of the four subjects also displayed small rightward tilts to their mean gradients. The initial depth used to maximize the fit quality varied from 0.57 to 0.9 mm across subjects.

We examined the effects of spatial smoothing on our data by blurring the original fMRI data volumes with a 2-mm FWHM Gaussian kernel, then repeating the same depth analysis (Figure S2 in Supplementary Material). Significant positive depth gradients are still evident, with spatial structure largely similar to that observed in the original data. However, the blurring tends to reduce the absolute magnitude of the depth gradients, and increases the area of reversed gradients observed at the margins of the colliculi.

Frequency tuning was examined by plotting the tuning width as a function of centroid frequency (**Figure 6**). The results show an "inverted-U" shape for all individual subjects, with a similar profile for all subjects taken together. Each data point shows tuning width for a single IC voxel, while the red lines show a regridded average with bootstrapped 68% confidence intervals. On average, tuning widths are about 40% broader at intermediate frequencies than they are at the lower and higher ends of the stimulus frequency band. Bootstrapped reliability estimates indicate that tuning widths at both the lowest centroid frequencies (0.25–0.5 kHz) and highest centroid frequencies (4–8 kHz) are both significantly lower (*p* < 0.02) than tuning widths at intermediate frequencies (1.5–3 kHz) for each individual subject.

## **DISCUSSION**

We used high-resolution imaging and surface-based analyses at traditional 3 T field strength to assess frequency representation within the human IC. Maps of centroid frequency showed a consistent depth-wise frequency organization in portions of both ICs of all subjects. To examine the laminar topography of the IC, we computationally modeled the surface of the IC, which permitted electrode-like "penetrations" into the depth of tissue. We observed significant tonotopy in the depth of IC tissue, with superficial layers showing greater preference to low frequencies and deeper layers to high frequencies. The depth gradients were evident by examining a simple subtraction of response to highfrequency stimuli from that evoked by low-frequency stimuli, or by evaluating centroid frequencies for each voxel sampled.



While several studies have documented midbrain activity to auditory stimulation (Guimaraes et al., 1998; Hawley et al., 2005; Sigalovsky and Melcher, 2006; Rinne et al., 2007; Chandrasekaran et al., 2012), these studies could not determine the complex 3D functional topography of the IC because of their limited spatial resolution. More recently, neuroimaging studies have used very high-field fMRI to evaluate tonotopy at the IC. Baumann et al. (2011) showed a dorsolateral-to-medioventral tonotopic gradient in rhesus monkeys at 4.7 T. Similar results were found in humans for 7 T imaging (De Martino et al., 2013). In the current study, the use of an auditory-optimized fMRI pulse sequence at a resolution of 1.2 mm enabled frequency analysis within the IC at 3 T field strength. In addition, surface-based analysis allowed the examination of functional variation in frequency preference across the laminar depth of IC, which can be viewed as a hemodynamic analog of electrode depth penetrations.

The mean depth gradient for 3/4 subjects was fairly well aligned to a vector pointing toward the neuraxis. Subject 2, with perhaps the weakest data of the four, displayed a mean gradient that was almost horizontal. All subjects also exhibited small left-right asymmetries in the gradient direction, probably reflecting a combination of lateral-medial and left-right asymmetries. Thus, the notion of using tissue depth as a reference frame was proven to work fairly well, despite the patchy character of the data. However, IC has a great deal more curvature than SC, so it is not surprising that we see some regions of reverse gradients at the margins of the colliculi, particularly on caudal regions where curvature is largest. Another drawback of the depth-penetration approach is that voxels in deeper regions tend to get more heavily resampled by adjacent penetrations because of curvature effects.

Our results obtained from humans using non-invasive high (1.2 mm) resolution are largely consistent with invasive animal measurements of frequency preference in CN (Schreiner and Langner, 1997) as well as high-field neuroimaging studies, which indicate that individual neuronal frequency preferences are organized systematically by depth. In contrast to these studies, the detailed spatial structure of our preferred-frequency measurements shows a patchy character. The patchy nature of the topographical representation in our study may reflect the different neural sensitivity of BOLD fMRI. Specifically, various experiments indicate that fMRI is more sensitive to neuronal inputs (Logothetis et al., 2001), or to local-circuit activity (Angenstein et al., 2009). Accordingly, our results may more closely reflect the patchy, somewhat interdigitated tonotopic organization of bottom-up neuronal inputs to IC from the cochlear nucleus, as reported by anterograde labeling studies (Oliver and Morest, 1984; Hofstetter and Ehret, 1992; Loftus et al., 2008).

The patchy character of our results disagree with the much more continuous spatial gradient of tonotopic selectivity observed by Baumann et al. (2011) and De Martino et al. (2013). The observed gradients, by construction, are oriented normal to the tissue surface, but the aggregate weighted mean shows an angle that is oriented more perpendicular to the neuraxis than the nearly superior-inferior gradient direction observed by De Martino et al. The observed structure is unlikely to be consequence of noise; they are a real feature of the data. Accordingly, they are likely either to reflect one or more neural or hemodynamic effects. On the neural side, the observed patchiness could be the consequence of taskrelated recruitment of sub-regions of IC tissue. Task demands such as selective attention have been shown to affect IC (Rinne et al., 2008), but such effects have not yet been evaluated at these fine spatial scales. The observed patchiness could also be hemodynamic in character, as related to the particular vascular compartment sensitivities of the BOLD response evoked at 3 T field strength as compared to higher fields (Gati et al., 1997; Yacoub et al., 2001; Duong et al., 2003; van der Zwaag et al., 2009). Specifically, at 3 T, only about half of the BOLD contrast corresponds to the extravascular component that is dominant in the capillary parenchyma, while half of the signal comes from the intravascular component in larger blood vessels. Although the relative contribution of venous and capillary signals was estimated to be similar at 3 and 7 T in cerebral cortex (van der Zwaag et al., 2009), the situation could be different in subcortical tissue. If so, the observed

patchy structure may reflect stronger contrast from regions of IC where there are larger concentrations of draining veins, while little contrast is observed in other regions.

Our high-resolution data had comparatively high levels of thermal noise. Conventional fMRI studies operate where the thermal SNR is typically >100; our images had SNR ∼20–30. This noise could also have caused the observed patchiness of the results. To address the effects of thermal noise on our results, we performed a parallel analysis in which the data were first blurred with a 2-mm FWHM Gaussian spatial kernel. The post-blurring is expected to improve thermal SNR by the square root of the voxel-volume ratio, that is (2/1.2)1.5 = 2.2; this increase in SNR should greatly reduce the degradation of the thermal noise. The spatial structure of the results was qualitatively similar: a patchy distribution of tonotopically selective regions on the IC surface. This result indicates that the observed patchiness is a real feature of the data, and not the consequence of random noise. This similarity is also consistent with the bootstrapped evaluation of the data, which indicate the tonotopic gradients are significant only within the patchy regions presented in the results. The invariability of the tonotopy to the spatial resolution is also consistent with the results of higher-field fMRI studies (Baumann et al., 2011; De Martino et al., 2013).

In particular, we tend to observe weaker tonotopic depth gradients in the very center of the IC, an area expected to have a clearly organized depth tonotopy. There are at least three possible reasons for this lack of clear tonotopy. First, tonotopic inputs to this region may too finely interdigitated to be resolved by our measurement. If, in fact, the measured hemodynamic response were evoked more strongly by the inputs to this portion of the colliculus, then the heterogeneity of responses in each voxel would tend to blur the frequency selectivity. Second, it is again possible that local vascular characteristics modulate the fMRI responses, yielding weaker contrast in the central portions of each colliculus. Third, the task demands of our protocol may tend to selectively recruit IC circuits in the observed sub-regions. Further experiments will be necessary to investigate these issues. In particular, it will be important to investigate the repeatability of the observed patterns across sessions.

Our stimulus protocol utilized a repetitive, monotonic cycle of stimulus frequencies. This choice was motivated by our original interest in a phase-encoded analysis (Engel et al., 1997; Barton et al., 2012), but has several drawbacks when used in conjunction with the current approach. In particular, the undershoot of the HRF evoked by each stimulus will impact the hyperoxic peak of the subsequent stimuli. However, the monotonic sequence somewhat softens this effect: the 0.25-kHz response undershoot will reduce the response to the 0.59-kHz stimulus, the 0.59-kHz response undershoot will reduce the response to the 1.41-kHz stimulus, and so forth. Thus, the undershoot acts as a form of blurring in the frequency domain, with relatively benevolent effects.

We did not observe narrow frequency tuning in our voxels, particularly at intermediate frequencies. The most likely explanation for this result is the spatial resolution of our measurements. Each voxel encompasses a volume of ∼1.7µL containing several million neurons, so some heterogeneity of response can be expected. Interestingly, we observed sharper tuning at low and high frequencies, suggesting that the population response for these frequency bands

is more spatially homogeneous than at intermediate frequencies. Also, the observed depth gradients did not span the full range of the stimuli. Again, this is likely the consequence of the available spatial resolution. Sharply tuned responses from both the highest and lowest frequency stimuli cannot be resolved, so the preferredfrequency depth gradients show a smaller range. Nonetheless, we do resolve the presence of a tonotopic gradient as a function of depth, and bands of tissue with clearly evident preferences for low, intermediate, and high frequencies.

Tonotopic depth gradients were particularly evident in the left colliculus of all subjects, but were weaker in the right colliculus for three of our four subjects. There are two possible reasons for these left-right asymmetries. First, they may be a consequence of purely vascular effects. In our superior colliculus neuroimaging work (Katyal et al., 2010), we found that visual responses for a particular subject were frequently much greater in one colliculus or the other, without any obvious neural correlate. Second, such asymmetries (relatively stronger leftward bias) could be a consequence of the task used in this study. Auditory tasks that require attentive processing are known to show a left-hemispheric bias (Bryden et al., 1983) and can modulate IC activity via corticocollicular connections (Rinne et al., 2007, 2008). Corticocollicular connections are largely ipsilateral (from the auditory cortex) (Winer, 2006), and connect A1 neurons to neurons in the IC that have a similar frequency profile (Lim and Anderson, 2007). Thus, it is possible that the left-dominant collicular activity could reflect corticocollicular effects. Future studies need to compare the effects of task (e.g., compare passive listening to active task), and to evaluate the role of top-down corticocollicular pathways in determining frequency selectivity within the IC.

Examining the IC encoding of frequency information is of significant clinical relevance. The IC has been considered as an alternate site for implant in hearing-impaired individuals who have contraindications for cochlear implants. The IC implant was developed based on animal models that suggest tonotopic organization in the CN. In fact, a few patients have already received IC implants and their success has been varied (Lim and Anderson, 2007). A recent study showed that initial frequency resolution in an individual with an IC implant was poor, but more precise selectivity was found 4 months after implantation (Lim et al., 2013). While our data demonstrates tonotopic representation as a function of depth, we find that tonotopic gradients are patchy and variable across individuals, which could pose a significant challenge for the placement of midbrain implants. Accordingly, further studies using high-resolution fMRI could provide a critical tool for surgical planning.

Our high-resolution surface-based fMRI methods and results provide a foundation for further non-invasive experimental work in small auditory subcortical structures. The ability to localize frequency selectivity with the depth of IC opens up new frontiers in auditory neuroscience, permitting the quantification of functional response in humans using complex stimuli and experimental paradigms. Moreover, the high-resolution sparse-sampling also enables investigations in other brainstem nuclei such as the cochlear nucleus and superior olivary complex. For clinical uses, these results open up the use of readily available 3 T MRI scanners for surgical planning and other treatment options. In conclusion, our methods are a significant advance for high-resolution fMRI research of auditory function in general, and will permit a host of critical auditory neuroscience experiments in small subcortical auditory structures that have been inaccessible to conventional auditory fMRI studies.

## **ACKNOWLEDGMENTS**

We gratefully acknowledge grant support from NSF BCS 1063774, sequence development assistance from Rez Khan, data analysis assistance from Seth Koslov, and helpful discussions with Dr. Sari Andoni.

## **REFERENCES**


## **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/Human\_Neuroscience/10. 3389/fnhum.2013.00586/abstract

## **Figure S1 | Additional color overlays of centroid frequency on anatomy**

**slices**. **(A,B)** Sagittal and coronal views of subject 1, respectively. **(C–E)** Sagittal, coronal, and axial views of subject 3. Inset surface figures show slice locations on surface models of each subject's midbrain.

**Figure S2 | Depth gradients obtained from data blurred with 2-mm FWHM kernel**. Same format at **Figure 5**.

*Reson. Med.* 39, 361–368. doi:10. 1002/mrm.1910390305


processing in human auditory cortex but not the inferior colliculus. *Neuroreport* 18, 1311–1314. doi:10. 1097/WNR.0b013e32826fb3bb


mechanisms of speech perception in noise. *J. Speech Lang. Hear. Res.* 51, 1026–1041. doi:10.1044/1092- 4388(2008/075)


single unit properties and neuronal architecture. *J. Comp. Neurol.* 231, 530–546. doi:10.1002/cne. 902310410

**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: 01 May 2013; accepted: 30 August 2013; published online: 19 September 2013.*

*Citation: Ress D and Chandrasekaran B (2013) Tonotopic organization in the depth of human inferior colliculus. Front. Hum. Neurosci. 7:586. doi: 10.3389/fnhum.2013.00586*

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

*Copyright © 2013 Ress and Chandrasekaran. 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.*

## Functional imaging of the human brainstem during somatosensory input and autonomic output

## **Luke A. Henderson<sup>1</sup>\* and Vaughan G. Macefield<sup>2</sup>**

<sup>1</sup> Department of Anatomy and Histology, University of Sydney, Sydney, NSW, Australia

<sup>2</sup> School of Medicine, University of Western Sydney, Sydney, NSW, Australia

## **Edited by:**

Florian Beissner, Martinos Center for Biomedical Imaging, USA

## **Reviewed by:**

Alexander Ritter, Friedrich Schiller Universität Jena, Germany William Blessing, Flinders University, Australia

## **\*Correspondence:**

Luke A. Henderson, Department of Anatomy and Histology, University of Sydney, Anderson Stuart Building, F13, Camperdown, NSW 2006, Australia e-mail: lukeh@anatomy.usyd.edu.au

Over the past half a century, many investigations in experimental animal have explored the functional roles of specific regions in the brainstem. Despite the accumulation of a considerable body of knowledge in, primarily, anesthetized preparations, relatively few studies have explored brainstem function in awake humans. It is important that human brainstem function is explored given that many neurological conditions, from obstructive sleep apnea, chronic pain, and hypertension, likely involve significant changes in the processing of information within the brainstem. Recent advances in the collection and processing of magnetic resonance images have resulted in the possibility of exploring brainstem activity changes in awake healthy individuals and in those with various clinical conditions. We and others have begun to explore changes in brainstem activity in humans during a number of challenges, including cutaneous and muscle pain, as well as during maneuvers that evoke increases in sympathetic nerve activity. More recently we have successfully recorded sympathetic nerve activity concurrently with functional magnetic resonance imaging of the brainstem, which will allow us, for the first time to explore brainstem sites directly responsible for conditions such as hypertension. Since many pathophysiological conditions no doubt involve changes in brainstem function and structure, defining these changes will likely result in a greater ability to develop more effective treatment regimens.

**Keywords: trigeminal nuclei, dorsal column nuclei, pain, sympathetic nerve activity, rostral ventrolateral medulla**

## **INTRODUCTION**

Since the advent of human functional Magnetic Resonance Imaging (fMRI) in the early 1990s, thousands of scientific investigations have explored brain activation patterns during various challenges. However, only a few investigators have explored brainstem function in humans using fMRI techniques and a few have even extended this into exploring regional activation patterns within the spinal cord. There a numerous factors that hinder human brainstem exploration using fMRI, including the intricate parcellation of regional functionality so that small voxel sizes need to be employed, the resulting low signal to noise ratios of small voxels, susceptibility due to structures such as the pontine cistern, the pulsatile signal from cerebrospinal fluid, and respiratory related brainstem movement.

However, steady improvements in signal to noise ratios, due to improved scanner hardware, software, and increases in magnetic field strengths, the development of scanning techniques that lower signal dropout in susceptible regions, such as the basis pontis, as well as improved image analysis techniques – particularly software that allows improved spatial normalization of brainstem structures – have resulted in greatly improved accuracy and feasibility of brainstem fMRI imaging. Using a 3 T MRI scanner, we have found that, over the past 7 years, improvements in signal:noise have resulted in the ability to acquire images with significantly greater spatial resolution, i.e., raw voxel sizes have decreased from approximately 3 to 1.5 mm<sup>3</sup> . Of course the development of human MRI scanners with even higher field strengths, e.g., 7 T, has also resulted in improved spatial resolution, making it possible to collect images with raw voxel sizes of 1 mm<sup>3</sup> (**Figure 1**). In addition to greater spatial resolution, the ability to compare signal intensity changes between subjects relies heavily on the accuracy of spatial normalization software. The development of brainstem-specific normalization software has allowed for greater accuracy between subjects and greatly improved the feasibility and accuracy of brainstem fMRI, particularly in studies involving comparisons within and between groups of subjects. In a recent article, Beissner et al. (2011) compared various brainstem imaging analysis techniques and found that five methods were found to be sensitive for brainstem activation, including analysis using various spatial smoothing kernels and spatial normalization in various software packages.

Given these improvements in brainstem fMRI image acquisition and processing, we are now in a position where brainstem function in humans can be explored with greater accuracy and ease. This is extremely important, since it is likely that functional changes in specific brainstem regions underlie the expression of many significant clinical conditions. For example, investigations in experimental animals have revealed that changes within the dorsal horn of the spinal cord or trigeminal nuclear complex in the brainstem are important for the expression of chronic pain (D'Mello and Dickenson, 2008). Additionally, many animal studies have defined the brainstem circuitry responsible for the generation of spontaneous sympathetic drive, such as the nucleus of the solitary tract, caudal, and rostral ventrolateral medulla (Dampney, 1994; Colombari et al., 2001). It is possible that activity changes within

these brainstem regions are at least partially responsible for generating the increase in sympathetic drive associated with conditions like essential hypertension and obstructive sleep apnea. Although animal work has greatly enhanced our understanding of brainstem function, for the main part they have been performed in anesthetized, and often artificially ventilated preparations, and in disease models that do not precisely mimic those seen in humans. Given this, it is important that we begin to explore brainstem function in awake humans, to understand both normal physiology and the disturbed control seen in various pathophysiological states.

T1-weighted anatomical images also collected using a 7T MRI scanner.

Over the past decade, we and others have been exploring human brainstem function and structure in healthy individuals and in various patient populations. The following sections will describe some of the recent results we have obtained in healthy individuals in which we have explored brainstem processing of somatosensory stimuli and sympathetic motor output. These provide the basis for exploring changes in brainstem function associated with conditions such as chronic pain and hypertension.

## **SOMATOSENSORY PROCESSING**

## **NOXIOUS AND NON-NOXIOUS SOMATOSENSORY PROCESSING IN THE OROFACIAL SYSTEM**

Brainstem imaging provides us with the unique opportunity to investigate the processing of both noxious and non-noxious somatosensory information from the orofacial region. It is well known that primary afferent neurons carrying thermoreceptive and nociceptive somatosensory information from the orofacial region synapse in the spinal trigeminal nucleus (SpV), which extends throughout the lateral medulla and into the upper cervical spinal cord (Sessle, 2000). In contrast, orofacial discriminatory touch information is processed in the principle or chief sensory trigeminal nucleus, which lies within the dorsolateral pons. This differential afferent termination pattern allows for the investigation of nociceptive versus non-nociceptive processing within the human brainstem using high resolution fMRI.

Chronic pain, i.e., pain that extends beyond the expected period of healing and for more than 3 months, has a significant detrimental impact on an individual's life. Since most forms of chronic pain are resistant to currently available treatments, they often persist for years and even decades. The lack of effective treatment regimens results from the fact that little is known about the underlying mechanisms, particularly those that result from nervous system damage (neuropathic pain). Given that animal evidence suggests that changes at the level of the primary synapse are involved in the generation and/or maintenance of some forms of chronic pain, and that in humans some pain treatments are thought to exert their analgesic effects via actions on the primary synapse (e.g., opiates, transcutaneous electrical nerve stimulation), it is important that we begin to explore the processing of noxious information at the primary synapse so that the underlying mechanisms of chronic pain are understood and more effective treatment regimes are developed.

With this in mind, we and others have begun to use fMRI to explore neural activation patterns during acute noxious and non-noxious somatosensory stimulation in healthy individuals. In a recent series of investigations, we measured brain activity using a 3 T MRI scanner during non-noxious somatosensory activation evoked by lightly brushing various parts of the body in healthy individuals (Wrigley et al., 2009; Henderson et al., 2011; Gustin et al., 2012). Using Blood Oxygen Level Dependant (BOLD) imaging, we collected 130 fMRI volumes covering the entire brain, during which we brushed the big toe or thumb using a repeated "On-Off" paradigm (seven Off, six On periods). Although the aim of this investigation was to explore changes within the primary somatosensory cortex, in five individuals our fMRI scans also encompassed the caudal medulla, including the level at which non-noxious somatosensory afferents from the body terminate. In a preliminary investigation we explored medullary signal changes during thumb and toe brushing to determine if we could identify activation of the cuneate and gracile nuclei, respectively. Using brainstem-specific software (SUIT toolbox in SPM8) (Diedrichsen et al., 2011) we isolated and normalized the brainstem, applied global signal detrending (Macey et al., 2004) and then determined significant patterns of signal intensity increases during each brushing period using a repeated box-car model. A second-level conjunction analysis (*p* < 0.005, uncorrected for multiple comparisons) revealed the well-described somatotopic organization of non-noxious somatosensory termination patterns within the brainstem. That is, innocuous brushing of the big toe resulted in activation of the region of the ipsilateral gracile nucleus, whereas brushing of the thumb activated the region of the

ipsilateral cuneate nucleus (**Figure 2**). It can be seen in **Figure 2** that signal intensity increased during each brushing period and subsequently returned toward baseline levels during the rest periods. A similar activation pattern during innocuous stimulation of the upper limb has been shown previously (Ghazni et al., 2010). We also found signal intensity increases during both brushing paradigms in the region of the lateral reticular nucleus. It is known from experimental animal investigations that lateral reticular nucleus neurons receive ascending inputs from the dorsal horn and respond to noxious and non-noxious somatosensory stimulation (Kitai et al., 1967). Although the precise role of this region remains unknown, it projects to the cerebellar vermis and may play a role in co-ordinating movement, although a role for this region in cardiovascular regulation has also been described (Thomas et al., 1977). In addition to a somatotopic representation within the lower brainstem, we have previously shown a similar organization within the relevant recipient region of the thalamus and primary somatosensory cortex (Wrigley et al., 2009; Gustin et al., 2012).

As well non-noxious somatosensory activation patterns, we have explored brainstem activity evoked by noxious somatosensory stimulation involving two separate methods: acute muscle pain – evoked by a 0.5 ml injection of hypertonic (5%) saline into the belly of the right masseter muscle – and acute cutaneous pain, evoked by injection of hypertonic saline into the skin overlying the masseter muscle (Nash et al.,2009). We used the same image collection parameters and a similar analysis procedure as that employed for the experiments involving non-noxious somatosensory stimulation. We know from animal work that the SpV is divided into three major divisions: oralis, interpolaris, and caudalis, from rostral to caudal (Sessle, 2000). Although the roles of these subnuclei have begun to be explored in experimental animals, it remains

**FIGURE 2 | Medullary activation during repeated innocuous brushing of the thumb and big toe in five subjects**. Note that thumb brushing activates a region of the ipsilateral dorsal medulla that is lateral and marginally rostral to that region activated by brushing of the big toe. These activation patterns are consistent with the well-described termination patterns of non-noxious somatosensory afferent neurons, i.e., the cuneate and gracile nuclei. Note that during each of the six brushing period (vertical gray bars), the mean (±SEM) percentage change in signal intensity increases. The slice locations in Montreal Neurological Institute space are shown at the lower left of each section.

unknown if differences in noxious processing occur in these different subnuclei in humans. In our investigations, we found that although acute pain evoked significant increases in signal intensity within the ipsilateral SpV, the pattern of activation differed depending on whether the noxious stimulus was delivered to skin or muscle. Whereas both cutaneous and muscle pain activated the caudalis and oralis divisions of SpV, only cutaneous pain activated the interpolaris division (**Figure 3**) (Nash et al., 2009). Furthermore, we found that orofacial muscle pain – but not cutaneous pain – evoked increases in signal intensity in the region of the pons that was also activated by innocuous lip brushing, i.e., the principal sensory nucleus (Vp). It is possible that these differential patterns of activation reflect differences in higher processing of the perceptual, emotional, and motoric consequences of deep compared with superficial pain (Lewis, 1942).

**FIGURE 3 | Brainstem activation during noxious stimulation of the masseter muscle and the skin overlying the masseter muscle in 18 healthy subjects**. Note that whereas cutaneous and muscle pain evoked activation of both the caudalis (data not shown) and oralis divisions of the spinal trigeminal nucleus (SpVo), only cutaneous pain activated the interpolaris division of SpV (SpVi). Furthermore, only muscle pain activated the principle trigeminal nucleus (Vp). The slice locations in Montreal Neurological Institute space are shown at the lower left of each section. Plots of mean (±SEM) percentage change in signal intensity (SI) over time are shown for each region. The vertical dashed lines indicate the onset of either the cutaneous (green) or intramuscular (gray) injections of hypertonic saline.

In addition to differential activation patterns according to the type of tissue stimulated, it is well known that the caudalis division of SpV is somatotopically organized, with inputs from the three trigeminal nerve divisions – ophthalmic, maxillary, and mandibular – being arranged in a rostrocaudal fashion. Of course, this termination pattern is responsible for the unique "onion skin" pattern of pain and temperature perception loss that occurs in humans following a lesion to a discrete rostrocaudal level of the spinal trigeminal caudalis nucleus. Remarkably, somatotopic organization of the SpV, as well as the trigeminal ganglion, has been reported using fMRI. Over 10 years ago, Borsook et al. (2003) assessed changes in fMRI signal intensity within the trigeminal ganglion during innocuous brushing and noxious thermal stimulation of the face within the three trigeminal nerve receptive fields. Remarkably, they found that stimulation in these different receptive fields evoked differential activation patterns within the trigeminal ganglion, consistent with the known afferent locations derived from anatomical investigations.

Even more remarkable are the recent investigations that use fMRI to explore activity changes within the dorsal horn of the spinal cord. For example, using a 3 T scanner, Nash et al. (2013) collected BOLD images of the cervical spinal cord using a spiral inout image sequence with retrospective correction for physiological noise, to investigate activation patterns within the cervical spinal cord during noxious stimulation of the skin overlying the deltoid muscle and thenar eminence. They found that noxious stimulation evoked increases in signal intensity in the region of the ipsilateral dorsal horn at the relevant spinal cord level (deltoid: C4–C5; thenar eminence: C6–C7). These human fMRI studies show that it is now possible to investigate activation patterns throughout the entire "pain" pathway, from the sensory ganglion and dorsal horn/SpV, rostrally to the thalamus and up to cortical regions. This opens up the possibility to explore changes at all levels of the nervous system in individuals with, for example, chronic pain. Such an improvement in our understanding of the underlying neurophysiological processes responsible for the generation of acute pain may ultimately lead to better treatment of chronic pain.

## **AUTONOMIC OUTPUT**

The brainstem contains several nuclei involved in autonomic control, the most studied being those which regulate the cardiovascular system. Since hypertension is a significant risk factor for the development of a number of serious medical conditions, and neurogenic hypertension is brought about by an increase in sympathetically mediated constriction of muscle (and splanchnic) vascular beds, it is critical that we begin to explore the underlying circuitry responsible for generating sympathetic vasoconstrictor drive in humans. Although thousands of investigations have defined brainstem and higher circuitry involved in generating and/or modulating arterial pressure, heart rate, and sympathetic nerve activity in anesthetized experimental animals, relatively few have explored this circuitry in awake human subjects (Macefield et al., 2006; Macefield and Henderson, 2010; Sander et al., 2010). Additionally, despite much experimental focus on the causes of hypertension in humans, the underlying mechanisms responsible for most hypertensive conditions remains almost completely unknown. Indeed, essential hypertension, which defines

about 95% of all hypertensive patients, has no identifiable cause, although a number of risk factors have been identified. Furthermore, despite the development of antihypertensive drugs, in most individuals, hypertension and its associated risk factors remain uncontrolled (Messerli et al., 2007). If effective control of hypertension is to be obtained in all individuals with essential hypertension or in individuals with other conditions associated with high blood pressure, such as obstructive sleep apnea (Carlson et al., 1993),an understanding of the underlying circuitry responsible for generating baseline sympathetic drive needs to be obtained in normotensive individuals before we begin to explore potential changes in this circuitry in individuals with various forms of hypertension.

Muscle sympathetic nerve activity (MSNA) can be recorded directly in via tungsten microelectrodes inserted percutaneously into a peripheral nerve in awake human subjects (microneurography). Developed over 40 years ago the technique has provided a wealth of information on the human sympathetic nervous system and its role in the control of blood pressure (Vallbo et al., 2004). It is well established that there is a tight relationship between dynamic changes in arterial pressure and MSNA through the operation of the baroreflex. Spontaneous variations in arterial pressure are detected by high-pressure arterial baroreceptors in the aortic arch and carotid sinuses, which send signals to the medulla to evoke opposing changes in MSNA. This reflex acts to cushion the original change in pressure and maintain resting arterial pressure at a stable level. The neural circuitry responsible for the baroreflex has been extensively investigated over many decades in anesthetized experimental animals. Baroreceptor afferents, traveling in the glossopharyngeal and vagus nerves, project to the caudal region of the nucleus tractus solitarius (NTS) which in turn provides excitatory drive onto tonically active inhibitory neurons in the region of the caudal ventrolateral medulla (CVLM). This region then inhibits excitatory neurons within the region of the rostroventrolateral medulla (RVLM), which contains premotor neurons which descend down the spinal cord to synapse with sympathetic preganglionic neurons within the intermediolateral cell column throughout the thoracic and upper lumbar cord (Dampney,1981). These investigations have also revealed that activity of the RVLM is critical for the maintenance of resting vasomotor tone and arterial pressure, since RVLM destruction results in precipitous falls in resting arterial pressure (Dampney et al., 1982). Furthermore, it appears that the RVLM provides the major output nucleus from which almost all brain regions, including the cerebral cortex, control arterial pressure. For example, during physical and psychological stressors, higher brain regions such as the cingulate, insular, and prefrontal cortices project to the RVLM to alter arterial pressure and sympathetic drive (Dampney et al., 2002; Gabbott et al., 2005). Given this, the RVLM and its contacts are likely candidates in which altered structure and function may result in hypertension in various medical conditions.

Whilst the medullary circuitry responsible for the baroreflex has been extensively studied in experimental animals, for the vast majority of these studies the animal was anesthetized and often ventilated. Although this reflex may work the same way in awake humans, it is known that in at least some brain regions, anesthesia has a significant effect on the cardiovascular response to stimulation. For example, chemical stimulation of the nucleus raphe obscurus in the medulla evokes increases in arterial pressure and heart rate under halothane anesthesia, yet decreases under urethane anesthesia (Henderson et al., 1998; Heslop et al., 2002). Furthermore, given the lack of effective hypertensive animal models, at least ones that mimic hypertension in conditions such as obstructive sleep apnea and essential hypertension, it is important that we investigate this circuitry in normotensive awake humans such that we can eventually elucidate changes in brain structure and function in individuals with various forms of hypertension.

## **BRAINSTEM ACTIVITY DURING EVOKED CHANGES IN MUSCLE SYMPATHETIC NERVE ACTIVITY**

Over the past decade a few investigators have begun to explore brain sites responsible for changes in sympathetic drive. We and others have indirectly measured human baroreflex functioning by using maneuvers such as lower-body negative pressure (Kimmerly et al., 2005), Valsalva maneuvers (Harper et al., 1998; Henderson et al., 2002), static hand grip (Sander et al., 2010), and maximal inspiratory apneas (Macefield et al., 2006). Although these studies were not designed or analyzed with brainstem-specific techniques, some did find that increased sympathetic drive was associated with increased signal intensity in the region of the RVLM, as well as a number of higher brain regions. However, there are some considerable methodological issues associated with these maneuvers that increase sympathetic nerve activity. While the sympathetic responses to challenges such as maximal inspiratory apneas and Valsalva maneuvers are relatively large and robust, they are also associated with significant subject movement, and may require significant skeletal motor activity to perform the actual task and often large and rapid increases in arterial pressure. Even though these

issues can be somewhat overcome by various analytical procedures, they do complicate the interpretation of the data. Furthermore, in these investigations, MSNA (and often arterial pressure) changes were recorded during a separate session, a situation that does not

allow for subtle difference in brain function to be assessed. In an attempt to remove these potential confounding issues, we have recently managed to successfully record brainstem fMRI and MSNA (using microneurography) concurrently in awake normotensive human at rest (Macefield and Henderson, 2010). By using a scan repetition time of 8 s and implemented a 4 s image collection period followed by a 4 s rest period we were able to collect brainstem BOLD images and MSNA in concurrent 4 s periods. This protocol was chosen to take advantage of the temporal delays inherent in BOLD imaging. Since the peak of the hemodynamic delay is estimated to lag neural activity by approximately 5 s (Logothetis et al., 2001), and given that we need to allow approximately 1 s for a muscle vasoconstrictor volley to travel from the brainstem to the peripheral recording site at the knee (Fagius and Wallin, 1980), the changes in BOLD signal intensity within the brainstem would reflect changes in neural activity associated with emission of sympathetic volleys recorded in the previous 4 s epoch (**Figure 4**). Furthermore, we collected axial slices sequentially from caudal to rostral so that the timing of each slice could be related to the MSNA recording.

## **BRAINSTEM ACTIVITY DURING SPONTANEOUS FLUCTUATIONS IN MUSCLE SYMPATHETIC NERVE ACTIVITY**

Using the protocol described above, we explored the medullary circuits responsible for the baroreflex by measuring regional brainstem activity changes during small, spontaneous changes in MSNA during 10 recording sessions in 8 awake normotensive humans. Using brainstem-specific analysis techniques, we determined which brainstem regions displayed high or low signal intensity during periods of high MSNA. We found that small increases in MSNA were associated with fMRI signal intensity increases in the RVLM and signal decreases within the regions of the CVLM and NTS (**Figure 5**). That is, significant correlations between regional brainstem signal intensity and MSNA were found to occur in those regions previously described as being responsible for baroreflex functioning in anesthetized experimental animals, the RVLM, CVLM, and NTS (Macefield and Henderson, 2010). The increase in signal intensity within the region of the RVLM coincides with similar activation patterns evoked by maximal inspiratory capacity apneas and sustained hand grip, challenges which are also associated with significant increases in MSNA (refs). Furthermore, these activations lie in the same region defined anatomically as the RVLM by the binding of angiotensin receptors (Allen et al., 1998) (**Figure 6**).

The identification of sites within the brainstem, and in particular the medulla, are the first step in elucidating brainstem and higher changes associated with hypertensive conditions. It has been shown that essential hypertension and obesity (a risk factor for hypertension and obstructive sleep apnea, a condition with comorbid hypertension) are associated with altered brain structure including gray matter volume changes in areas such as the insula cortex and hippocampus (Macey et al., 2002; Taki et al., 2008). It is possible that in these and other conditions, changes in

**FIGURE 5 | Functional magnetic resonance imaging (fMRI) signal intensity changes correlated to spontaneous fluctuations in muscle sympathetic nerve activity (MSNA) in eight subjects**. Increases in signal intensity with increases in MSNA are coded by the hot color scale and signal decreases with the cool color scale and are overlaid onto a series of sagittal and axial fMRI slices from an individual subject. Myelin stained sections through the medulla are shown to the right. CVLM, caudal ventrolateral medulla; NTS, nucleus tractus solitaries; RVLM, rostral ventrolateral medulla; SI, signal intensity.

higher brain structures that occur either before or during the early pre-clinical stages of the condition, result in altered influences over the RVLM, i.e., reduced inhibition or increased excitation, resulting in increased RVLM activation and a resulting increase in arterial pressure and on-going MSNA. Future investigations should firstly examine changes within the region of the RVLM and subsequently changes in higher brain structures that may influence RVLM function.

## **BRAINSTEM ACTIVITY DURING EVOKED CHANGES IN SKIN SYMPATHETIC NERVE ACTIVITY**

Finally, in addition to work exploring brainstem sites responsible for evoked and spontaneous changes in MSNA, recent investigations have also explored changes in brain activity, including the brainstem, during changes in skin sympathetic nerve activity (SSNA) or sweating, a rather sluggish marker of SSNA. The sympathetic innervation of the skin in humans primarily subserves thermoregulation, by constricting or dilating cutaneous blood vessels and increasing or decreasing the release of sweat. In thermoneutral conditions, resting SSNA is primarily related to the level of arousal and emotional state (Hagbarth et al., 1972). We recently recorded whole brain fMRI and SSNA concurrently in 13 healthy individuals whilst viewing emotionally charged images (erotica and mutilation) (Henderson et al., 2012). We found that viewing emotionally charged images increased SSNA, and that this SSNA increase was associated with altered activity in brain regions such as the amygdala, pons, thalamus, nucleus accumbens, cerebellar

cortex, and precuneus. Additionally, we recently made concurrent SSNA and fMRI recordings in 13 awake human subjects to identify brain regions associated with the generation of spontaneous SSNA (James et al., 2013). Spontaneous SSNA increases were associated with signal changes in the thalamus, insula, cingulate cortex, and precuneus.

In recent publications by McAllen et al., brainstem activation patterns were explored during changes in skin temperature and during psychological stressors. They found that skin cooling evoked fMRI signal intensity increases in rostral medullary raphe region, a region shown to contain thermoregulatory neurons in rodents (McAllen et al., 2006). Further, they found that thermal (whole-body heating,11 subjects) and psychological stressors (Stroop test; 11 subjects), evoked increased sweating and were associated with increased signal intensity in the rostral lateral midbrain and in the rostral lateral medulla (Farrell et al., 2013). Although the site of activation within the lateral medulla

## **REFERENCES**


*Hypertension* 38, 549–554. doi:10. 1161/01.HYP.38.3.549


was located close to that which we previously identified as the RVLM, the authors propose that the activation is located in the region between the facial nuclei and pyramidal tracts, an area that corresponds to a neuronal group found to evoke sweating in experimental animals.

## **CONCLUSION**

While the brainstem – and especially the medullary component – present significant difficulties for neuroimaging, there is a clear need to understand the functional organization of this phylogenetically ancient structure. With the advent of higher field strengths (e.g., 7 T), and improvements in spatial resolution, our capacity to understand the operation of the brainstem will increase enormously. Nevertheless, it is remarkable that individual nuclei can be functionally identified even at 3 T, and that the microcirculation within the brainstem allows BOLD imaging in such a small volume.

Gizewski, E. R., et al. (2011). Imaging the deep cerebellar nuclei: a probabilistic atlas and normalization procedure. *Neuroimage* 54, 1786–1794. doi:10.1016/j. neuroimage.2010.10.035


(2005). Cortical regions associated with autonomic cardiovascular regulation during lower body negative pressure in humans. *J. Physiol.* 569, 331–345. doi:10.1113/jphysiol.2005. 091637


McKinley, M. J., et al. (2006). Human medullary responses to cooling and rewarming the skin: a functional MRI study. *Proc. Natl. Acad. Sci. U.S.A.* 103, 809–813. doi: 10.1073/pnas.0509862103


lateral reticular nucleus in central cardiovascular regulation in the cat. *Am. J. Physiol.* 232, H157–H166.


**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: 13 June 2013; paper pending published: 23 June 2013; accepted: 26 August 2013; published online: 17 September 2013.*

*Citation: Henderson LA and Macefield VG (2013) Functional imaging of the human brainstem during somatosensory input and autonomic output. Front. Hum. Neurosci. 7:569. doi: 10.3389/fnhum.2013.00569*

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

*Copyright © 2013 Henderson and Macefield. 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.*

## High-resolution MR imaging of the human brainstem in vivo at 7 Tesla

#### **Andreas Deistung<sup>1</sup>\*, Andreas Schäfer <sup>2</sup> , Ferdinand Schweser <sup>1</sup> , Uta Biedermann<sup>3</sup> , Daniel Güllmar <sup>1</sup> , Robert Trampel <sup>2</sup> , Robert Turner <sup>2</sup> and Jürgen R. Reichenbach<sup>1</sup>**

<sup>1</sup> Medical Physics Group, Institute of Diagnostic and Interventional Radiology I, Center of Radiology, Jena University Hospital – Friedrich Schiller University Jena, Jena, Germany

<sup>2</sup> Department of Neurophysics, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany

3 Institute of Anatomy I, Jena University Hospital – Friedrich Schiller University Jena, Jena, Germany

## **Edited by:**

Florian Beissner, Martinos Center for Biomedical Imaging, USA

## **Reviewed by:**

Markus Barth, Radboud University Nijmegen, Netherlands Christine Preibisch, Klinikum rechts der Isar der Technischen Universität München, Germany

### **\*Correspondence:**

Andreas Deistung, Medical Physics Group, Institute of Diagnostic and Interventional Radiology I, Jena University Hospital – Friedrich Schiller University Jena, Philosophenweg 3 – MRI Building, 07743 Jena, Germany e-mail: andreas.deistung@med. uni-jena.de

The human brainstem, which comprises a multitude of axonal nerve fibers and nuclei, plays an important functional role in the human brain. Depicting its anatomy non-invasively with high spatial resolution may thus in turn help to better relate normal and pathological anatomical variations to medical conditions as well as neurological and peripheral functions. We explored the potential of high-resolution magnetic resonance imaging (MRI) at 7 T for depicting the intricate anatomy of the human brainstem in vivo by acquiring and generating images with multiple contrasts:T <sup>2</sup>-weighted images, quantitative maps of longitudinal relaxation rate (R<sup>1</sup> maps) and effective transverse relaxation rate (R<sup>2</sup> <sup>∗</sup> maps), magnetic susceptibility maps, and direction-encoded track-density images. Images and quantitative maps were compared with histological stains and anatomical atlases to identify nerve nuclei and nerve fibers. Among the investigated contrasts, susceptibility maps displayed the largest number of brainstem structures. Contrary to R<sup>1</sup> maps andT <sup>2</sup>-weighted images, which showed rather homogeneous contrast, R<sup>2</sup> <sup>∗</sup> maps, magnetic susceptibility maps, and track-density images clearly displayed a multitude of smaller and larger fiber bundles. Several brainstem nuclei were identifiable in sections covering the pons and medulla oblongata, including the spinal trigeminal nucleus and the reticulotegmental nucleus on magnetic susceptibility maps as well as the inferior olive on R1, R<sup>2</sup> ∗ , and susceptibility maps. The substantia nigra and red nuclei were visible in all contrasts. In conclusion, high-resolution, multi-contrast MR imaging at 7T is a versatile tool to non-invasively assess the individual anatomy and tissue composition of the human brainstem.

**Keywords: brainstem, quantitative susceptibility mapping, effective transverse relaxation, longitudinal relaxation, diffusion tensor imaging, track-density imaging, brain, anatomy**

## **INTRODUCTION**

The brainstem is the primary relay center for afferent and efferent connections between the cerebral cortex, the cerebellum, and the spinal cord (Sutin and Carpenter, 1976). It plays an important role in the regulation of vital functions (e.g., respiratory function, cardio-vascular function, nausea) and coordinates motor control signals sent from the brain to the body. Anatomically, the human brainstem consists of three connected parts, medulla oblongata, pons, and midbrain, and is composed of a multitude of axonal nerve fibers as well as cranial and non-cranial nerve nuclei. Due to the spatial concentration of important neural structures in this relatively small brain region, pathological variations or lesions in the brainstem can introduce severe and multiple neurological effects (Donaldson et al., 2006). The brainstem is also known to be involved in a number of degenerative diseases, such as

Alzheimer's disease, Parkinson's disease, and Huntington's disease (Urban and Caplan, 2011). Hence, detailed non-invasive depiction of its morphology and anatomy may help to better relate normal and pathological anatomical variations to medical conditions as well as neurological and peripheral functions.

Magnetic resonance imaging (MRI) at ultrahigh magnetic field strength (*B*<sup>0</sup> ≥ 7 T) is a powerful means to assess non-invasively normal and abnormal brain tissue with high spatial resolution (Li et al., 2006; Duyn et al., 2007; Deistung et al., 2008; Trampel et al., 2011; Moser et al., 2012; Turner, 2012, 2013; Eichner et al., 2013; Marques and Gruetter, 2013). MRI provides a variety of qualitative and quantitative tissue contrasts, mainly reflecting nuclear relaxation (*T*1,*T*2, and *T* ∗ 2 , or equivalently *R*1, R2, and*R* ∗ 2 ), diffusion, and magnetic susceptibility. In contrast-sensitized MR images (e.g., *T*1-weighted, *T*2-weighted, or *T* ∗ 2 -weighted images), however, the measured signal is an intricate function of proton density, longitudinal (*T*1) and transverse relaxation (*T*2) and depends on the applied MRI sequence parameters. Thus, it may be difficult to interpret the underlying tissue composition and structure based on the intensity of these images. Quantitative MRI

**Abbreviations:** DTI, diffusion tensor imaging; GRE, gradient recalled echo; QSM, quantitative susceptibility mapping; ppb, parts per billion; *R* ∗ 2 , effective transverse relaxation rate; SHARP, sophisticated harmonic artifact reduction for phase data; TDI, track-density imaging.

techniques, on the other hand, provide information that is intrinsically more tissue specific and consequently less dependent on the chosen scan parameters. Hence, quantification of different physical tissue properties by using quantitative MRI techniques also increases intra- and inter-individual comparability and enables objective measurements of disease related brain changes (Tofts, 2003).

The most common approach of quantitative MRI is relaxometry, i.e., mapping of the longitudinal (*R*1) and/or transverse (R2, *R* 0 2 , *R* ∗ 2 ) relaxation rates. In the brain these relaxation rates are influenced predominantly by both myelin (Koenig, 1991; Bender and Klose, 2010; Langkammer et al., 2012a; Lee et al., 2012) and tissue iron (Ogg and Steen, 1998; Langkammer et al., 2010). An additional,recently introduced novel quantitative MR method that is unique in its sensitivity to both tissue constituents is quantitative susceptibility mapping (QSM) (Schweser et al., 2011; Reichenbach, 2012). QSM is a post-processing method that uses the phase signal of complex-valued gradient (recalled) echo (GRE) data to produce maps of tissue magnetic susceptibility. Quantification of brain tissue iron content (Wharton et al., 2010; Schweser et al., 2011; Bilgic et al., 2012; Langkammer et al., 2012b; Zheng et al., 2013) and blood oxygenation *in vivo* (Haacke et al., 2010) as well as assessment of brain myelination (Liu et al., 2011; Li et al., 2012a,b) has recently been reported by applying QSM. Excellent anatomical delineation of cortical and deep gray matter structures and substructures has been demonstrated with QSM (Schweser et al., 2011, 2012; Deistung et al., 2013a), particularly in concert with relaxation rate mapping (Sigalovsky et al., 2006; Fukunaga et al., 2010; Geyer et al., 2011; Lebel et al., 2012; Deistung et al., 2013a), thus rendering this combined approach highly promising for improved delineation of the anatomical structure of the brainstem.

Nerve fibers and their complex network are usually imaged and analyzed by exploiting the translational molecular diffusion of protons. So far, most studies have applied diffusion tensor imaging (DTI) (Basser and Pierpaoli, 1996) for anatomical imaging of the brainstem (Stieltjes et al., 2001; Nagae-Poetscher et al., 2004; Salamon et al., 2005). DTI, however, is commonly limited by the relatively coarse spatial resolution available, which makes detailed analysis of smaller fiber tracts difficult, and often fails in regions of crossing fibers (Jones et al., 2013). One recently suggested approach to overcome this limitation is track-density imaging (TDI), which combines high angular resolution diffusion imaging (HARDI) with fiber tracking information to delineate structures beyond the image resolution of the acquired diffusionweighted data (Calamante et al., 2010). TDI has been already applied successfully to delineate fiber pathways in the cerebrum and cerebellum, and the substructure of the thalamus (Calamante et al., 2010, 2013).

The present study explores the potential of these highresolution MRI approaches to resolve non-invasively the intricate anatomy of the human brainstem at 7 T *in vivo*. To this end, *T*2 weighted images, quantitative *R*1, *R* ∗ 2 , and magnetic susceptibility maps as well as track-density images were compared with respect to their ability of depicting nerve nuclei and nerve fibers. Furthermore, *R*1, *R* ∗ 2 , and magnetic susceptibility values are presented for selected anatomical regions. It is shown that the combination of multiple image contrasts provides a distinctly improved portrayal of the morphology of the brainstem.

## **MATERIALS AND METHODS**

## **DATA ACQUISITION**

The study was approved by the internal Institutional Review Board of the Max Planck Institute in Leipzig and written informed consent was obtained from all participating subjects. A total of six healthy volunteers (four male and two female; 27.7 ± 3 years) were examined on a 7 T human whole body MRI system (Siemens Healthcare, Erlangen, Germany) using a 24-channel head array coil (NOVA Medical Inc., Wilmington, MA, USA). Two dielectric pads containing deuterated water (99%, Sigma Aldrich GmbH, Germany) enriched with calcium titanate (Alfa Aesar GmbH and Co KG, Karlsruhe, Germany) were placed on the left and right side of the subjects' head to increase both field strength and homogeneity of the transmit radio frequency field, *B* + 1 , in the area of the brainstem (Teeuwisse et al., 2012).

Coronal *T*2-weighted images were acquired using a twodimensional gradient-echo and spin-echo (GRASE) sequence (Feinberg and Oshio, 1991; Oshio and Feinberg, 1991) with echo time (TE) = 35 ms, repetition time (TR) = 10,000 ms, bandwidth (BW) = 343 Hz/px, acquisition matrix = 384 × 384 and in-plane resolution of 0.53 mm × 0.53 mm. Thirty-five images were collected with a slice thickness of 0.6 mm and a gap between two adjacent slices of 0.3 mm. For each TR interval there were three 180°-refocusing pulses, and between each 180°-pulse there were three gradient recalled echoes, resulting in an effective echo train length of nine and an acquisition time (TA) of 7:22 min:s. The GRASE sequence was applied to account for the high specific absorption rate (SAR) that substantially restricts *T*2-weighted imaging with turbo-spin-echo sequences at ultra-high magnetic fields.

High-resolution *R*1-mapping was performed based on the data acquired with the MP2RAGE sequence (Marques et al., 2010),a magnetization-prepared rapid gradient-echo (MP-RAGE) sequence with two different inversion times (TI). The acquisition parameters of the coronal MP2RAGE scan included TI<sup>1</sup> = 900 ms, TI<sup>2</sup> = 2750 ms, TE = 3.7 ms, TR = 5000 ms, BW = 240 Hz/px, acquisition matrix = 320 × 260 × 240, voxel size = 0.6 mm × 0.6 mm × 0.6 mm, and TA = 17:40 min:s. Magnetization preparation was achieved with a tailored radiofrequency pulse to take into account the heterogeneity of the radiofrequency transmit field (Hurley et al., 2010).

Coronal, multi-echo, three-dimensional gradient-echo imaging was carried out to compute quantitative *R* ∗ 2 and susceptibility maps. To this end, three echoes with monopolar readout were recorded (TE<sup>1</sup> = 11 ms, TE<sup>2</sup> = 21 ms, and TE<sup>3</sup> = 31 ms), TR = 43 ms, flip angle (FA) = 12.5°, BW = 149 Hz/px, acquisition matrix = 448 × 364 × 104, and voxel size = 0.43 mm × 0.43 mm × 0.43 mm. Data were collected with 75% and 87.5% partial Fourier along phase and slice encoding direction, respectively, resulting in an acquisition time of 17:48 min:s. Choosing a readout BW of 149 Hz/px yielded single-echo images with a rather high signal-to-noise ratio (SNR), however, at the expense of a rather long inter-echo distance. The signal decay was sampled with three echoes only to keep acquisition time below 20 min.

Finally, diffusion-weighted imaging (DWI) data were acquired in sagittal orientation with two-dimensional, single-shot, spinecho echo-planar imaging (EPI), 60 diffusion encoding directions each with a b-value of 1000 s/mm<sup>2</sup> and seven volumes with a *b*-value of 50 s/mm<sup>2</sup> , TE = 64 ms, TR = 10,000 ms, BW = 1050 Hz/px, 86 contiguous slices with an acquisition matrix of 170 × 170,voxel size = 1.2 mm × 1.2 mm × 1.2 mm,and TA = 14:20 min:s. To reduce geometric distortions partial parallel under sampling [GRAPPA (Griswold et al., 2002)] with an acceleration factor of 3 and 45 reference lines was applied in the phase-encoding direction.

In addition, a field-map was acquired with a 2D doubleecho gradient-echo sequence to correct residual geometric distortions of the diffusion-weighted images caused by susceptibility differences between air-bone and air-tissue interfaces. The field-map was recorded with TE<sup>1</sup> = 6 ms, TE<sup>2</sup> = 7 ms, TR = 1000 ms, FA = 50°, 51 contiguous slices, matrix = 102 × 102, voxel size = 2 mm × 2 mm × 2 mm, and TA = 3:22 min:s.

## **DATA PROCESSING**

## **Relaxometry**

Maps of the longitudinal relaxation times were calculated directly at the MRI scanner by applying the MP2RAGE reconstruction framework provided by the manufacturer. The relaxation times were subsequently converted to relaxation rates according to the relation *R*<sup>1</sup> = 1/*T*<sup>1</sup> to facilitate better depiction of anatomical structures.

Maps of the effective transverse relaxation rate, *R* ∗ 2 , were obtained using the power method (McGibney and Smith, 1993; Miller and Joseph, 1993), i.e., the squared magnitude signal decay, *S*(E*r*, *TE*) 2 , of the 3D multi-echo GRE scan was used for the regression (McGibney and Smith, 1993; Miller and Joseph, 1993):

$$\left(\mathcal{S}(\vec{r}, TE)^2 = \mathcal{S}\_0(\vec{r})^2 \cdot \exp\left(-2 \cdot TE \cdot \mathcal{R}\_2^\*\left(\vec{r}\right)\right),\tag{1}$$

where *S*0(E*r*) is the signal intensity at *TE* = 0 and E*r* is the position vector. This approach reduces contamination of Rician noise to the fit and provides more accurate relaxation rates than fitting a mono-exponential model to the magnitude signal decay *S*(E*r*, TE) (van der Weerd et al., 2000).

## **Quantitative susceptibility mapping**

Single-channel GRE magnitude images were combined using the sum-of-squares method (Roemer et al., 1990), whereas singlechannel GRE phase images were combined by taking into account the channel-dependent phase offset, which was estimated from the single-channel images at different echo times by voxel-wise linear fitting of the phase evolution (Robinson et al., 2011). Phase aliasing in the combined GRE phase data was resolved using a 3D Laplacian-based phase unwrapping algorithm (Schofield and Zhu, 2003) and phase images of different echo times were then combined in a CNR-optimized manner according to (Wu et al., 2012). Background phase contributions were eliminated with sophisticated harmonic artifact reduction for phase data (SHARP) (Schweser et al., 2011) (regularization parameter: 0.05) combining 10 different spherical kernels with varying radii ranging from 0.43 to 4.3 mm (Li et al., 2011). Susceptibility mapping was performed using homogeneity enabled incremental dipole inversion

(HEIDI), which incorporates *a priori* information extracted from the complex GRE signal to address the ill-posed nature of QSM (Schweser et al., 2012). Since the calculated susceptibility values represent relative rather than absolute values (Schweser et al., 2011), susceptibility differences were specified with respect to a homogenous region of normal appearing white matter (NAWM, semioval center).

## **Track-density imaging**

Geometric distortion of diffusion-weighted images was corrected using FUGUE (FSL toolbox, FMRIB, Oxford, England) based on the additionally acquired field-map. The unwarped DWI data were then processed using the MRtrix software package (Brain Research Institute, Melbourne, VIC, Australia) to produce track-density images as proposed in Calamante et al. (2010). Processing steps performed with MRtrix included the calculation of tensor data as well as fractional anisotropy maps, estimation of the response function, spherical deconvolution, probabilistic tractography and track-density mapping. Only voxels with fractional anisotropy, derived from the tensor model, that exceeded a threshold of 0.7 were taken into account for estimating the coefficients of the response function based on DWI data (Tournier et al., 2004). These coefficients and a maximum harmonic order of eight were used for constrained spherical deconvolution of DWI data (Tournier et al., 2007). This approach allowed modeling of multiple fiber populations within an imaging voxel, thereby overcoming the well-known limitation of the diffusion tensor model in regions with crossing fibers. Probabilistic streamline tracking was carried out using the second order integration over fiber orientation distributions (iFOD2) algorithm (Tournier et al., 2010) by seeding randomly throughout the predefined mask covering the whole brainstem and using the following relevant parameters: step size = 0.1 mm, number of fibers = 800,000, minimum curvature radius = 0.4 mm and a maximum number of 2000 trials for each point. In order to generate the track-density images, a virtual regularly spaced grid was superposed with a resolution of 0.43 mm × 0.43 mm × 0.43 mm on the tractography results and the total number of fiber tracks present in each grid element was calculated. Hereby, the superresolution property is achieved by utilizing the additional information provided by modeling the fiber tracking results (Calamante et al., 2010). The track-density images were calculated on a virtual isotropic resolution of 0.43 mm to facilitate comparison with the *R* ∗ 2 and susceptibility maps. Finally, the track-density images were direction-encoded by assigning each grid element an RGB color that represents the local fiber orientation as given by averaging the colors of all the fiber tracking (streamline) segments contained within each grid element.

## **DATA ANALYSIS**

All MRI contrasts were converted into an identical space for visual inspection. To this end, the *R*<sup>1</sup> maps, track-density images and *T*2-weighted images were linearly transformed into the space of the susceptibility maps based on orientation information encoded in the DICOM images. The susceptibility maps were chosen as target to maintain the high spatial resolution of the susceptibility and *R* ∗ <sup>2</sup> maps. Manual reformatting of quantitative susceptibility maps in multi-planar orientations was performed with Freeview of the Freesurfer software library<sup>1</sup> and the resulting rigid transformation matrix was applied to the previously registered MR image contrasts and to the *R* ∗ <sup>2</sup> maps. The MR images were averaged across three adjacent slices for signal-to-noise improvement, resulting in an effective slice thickness of 1.3 mm.

For qualitative comparison horizontal histological sections obtained from a 36-year-old male stained for myelin and for cells were consulted. These histological images were courtesy of the Brain Biodiversity Bank of the Michigan State University2,3 with support from the US National Science Foundation. An experienced neuroanatomist (UB; experience more than 15 years) assessed the MR images and identified anatomical structures and substructures if they coincided with a histoarchitectonic atlas (Bergman et al., 1989; Naidich et al., 2009; Paxinos et al., 2012) and if they were identifiable with respect to the surrounding tissue. For this qualitative analysis the window and level settings of the MR images were adjusted freely.

For quantitative characterization the contrast-to-noise ratio (CNR) of selected brain structures [red nucleus, substantia nigra, central tegmental tract, superior colliculus, inferior colliculus, reticulotegmental nucleus, middle cerebellar peduncle (MCP), superior cerebellar peduncle (SCP), and transverse pontine fibers] was determined across the different MR image contrasts (*R*1, *R* ∗ 2 , susceptibility, and *T*2-weighted contrast) according to:

$$\text{CNR}\_{\text{cy}} = \left| \frac{\text{S}\_{\text{cy}} - \text{S}\_{\text{Cys}}}{\sigma} \right| \cdot \eta\_{\text{c.}} \tag{2}$$

CNRcy denotes the CNR for MR image contrast *c* and anatomic region *y*. Mean image intensities measured in a volume-of-interest (VOI) of the investigated structure *y* and its surrounding tissue are denoted by *S*cy and *S*cys, respectively. The standard deviation of the signal intensities measured in a VOI of normal appearing white matter (semioval center) was used as an estimate of noise, σ. To account for variations in spatial resolution of the different investigated contrasts, the CNR was normalized to a voxel volume of 1 mm × 1 mm × 1 mm as suggested in (Nölte et al., 2012). To this end, the normalization factor η<sup>c</sup> is introduced in Eq. 2, which is defined as the ratio between 1 mm<sup>3</sup> and the acquired voxel volume of image contrast *c*. CNR was not determined on diffusion-weighted data because of their incompletely compensated geometric distortions.

To calculate the CNR and to measure both the relaxation rates and magnetic susceptibility, VOIs were identified based on magnetic susceptibility, *R*1, and *R* ∗ <sup>2</sup> maps in the coordinate space of the susceptibility maps for each subject in both hemispheres. The VOIs were then transferred to the space of the *R*<sup>1</sup> maps and the *T*2 weighted data, respectively, by employing orientation information encoded in the DICOM images and nearest-neighbor interpolation. Mean values and standard deviation of CNR and of the quantitative tissue parameters (*R*1, *R* ∗ 2 , and magnetic susceptibility) were then computed across hemispheres and subjects.

## **RESULTS**

All qualitative findings presented in this study were consistent across all subjects unless otherwise specified.

**Figure 1** presents representative MR images and corresponding histological sections of the midbrain. The red nuclei and the substantia nigra were most strikingly discernible on the *R* ∗ 2 (**Figure 1B**) and susceptibility maps (**Figure 1C**) as well as on the direction-encoded track-density (**Figure 1D**) and *T*2 weighted images (**Figure 1E**), but were barely visible on the *R*<sup>1</sup> map (**Figure 1A**). This qualitative finding was also supported by the CNR measurements presented in **Figure 2**. On the trackdensity image (**Figure 1D**) both nuclei showed lower anisotropy as reflected in the reduced color saturation (arrows b and d). Different fiber bundles traversing the crus cerebri [corticobulbar fibers (arrow c1), corticospinal fibers (arrow c2), and corticopontine fibers (arrow c3)] could be identified on the track-density images. In the midbrain, interestingly, the central tegmental tract (arrow f) and the medial lemniscus (arrow e, five of six subjects) were only identifiable reliably on susceptibility maps, whereas the mammillary body (**Figure 1**; arrow a) was discernible on all MR image contrasts except the track-density images. The superior (not shown) and inferior colliculi (**Figure 1**; arrow g) were distinguishable across all subjects only on the *R*<sup>1</sup> maps. Both of these structures could be assessed on the other MRI contrasts in at least four of six subjects (superior colliculus on susceptibility maps in six of six subjects). The superior colliculus and inferior colliculus presented highest CNR on susceptibility maps and *R* ∗ <sup>2</sup> maps, respectively (**Figure 2**). Anatomic regions with low magnetic susceptibility (**Figure 1C**; arrows c, e, f), such as the hypointense rim around the red nuclei coincided with regions of high myelin content (**Figure 1G**; arrows c1–3, e, f) and had no unique imaging correlate on any of the other MR images. The histological stains (**Figures 1G,H**) showed increased myelin and decreased cell density in the red nucleus compared to the substantia nigra.

MR images of the rostral and the middle part of the pons are depicted in **Figures 3** and **4**, respectively. The *R* ∗ 2 and magnetic susceptibility maps clearly displayed fibers traversing the pons in medial-lateral direction and fibers of the middle cerebellar peduncle (MCP) (arrow a in **Figure 3**) as also indicated by the myelin stains (**Figures 3G** and **4G**). These fiber structures were barely visible on the *R*<sup>1</sup> maps [**Figures 2** (TPF), **3A**, and **4A**]. The direction of the transverse pontine fibers could not be reliably separated from the superior-inferior running corticospinal tract on the track-density images (**Figures 3D** and **4D**). In contrast, the susceptibility maps revealed even smaller fiber tracts [e.g., central tegmental tract (arrow f, two of six subjects), medial longitudinal fasciculus (arrow e)] and nerve nuclei [reticulotegmental nucleus (arrow c)] that were indiscernible on *R*<sup>1</sup> and *R* ∗ <sup>2</sup> maps as well as on the *T*2-weighted images. In the direction-encoded trackdensity images larger [e.g., MCP (arrow a), superior cerebellar peduncle (arrow d), corticospinal tract (arrow h), pontocerebellar fibers (white solid outline in **Figure 3D**)] and smaller [e.g., central tegmental tract (arrow f)] fiber tracts were clearly identifiable. The MCP was also discernible on the *R*1, *R* ∗ 2 and susceptibility maps, whereas the superior cerebellar peduncle displayed on the *R*<sup>1</sup> and the susceptibility maps on all subjects and in two of six subjects on the *R* ∗ <sup>2</sup> maps. Both of these structures (MCP, SCP)

<sup>1</sup>http://surfer.nmr.mgh.harvard.edu/

<sup>2</sup>http://www.brains.rad.msu.edu

<sup>3</sup>http://brainmuseum.org

**FIGURE 1 | Images of the midbrain**. R1, R ∗ 2 , and susceptibility maps as well as direction-encoded track-density and T <sup>2</sup>-weighted images of the same brain region are presented in axial orientation from **(A–E)**, respectively. Red, blue, and green in the track-density image represent anisotropy along medial-lateral, superior-inferior, and anterior-posterior directions, respectively. Note that the track-density image is slightly distorted compared to the other image contrasts **(A–C,E)**. The sectional plane of the axial images is indicated

by the dashed line overlaid on the sagittal R<sup>1</sup> map shown in **(F)**. Axial

histological sections from a different subject stained for myelin and cells are illustrated in **(G)** and **(H)**, respectively. The arrows in the images indicate: (a) mammillary body, (b) substantia nigra, (c) crus cerebri, (c1) corticobulbar fibers, (c2) corticospinal fibers, (c3) corticopontine fibers, (d) red nucleus, (e) medial lemniscus, (f) central tegmental tract, (g) inferior colliculus, and (h) medial longitudinal fasciculus. [The histological stains **(G,H)** were adapted with permission from http://www.brains.rad.msu.edu and http://brainmuseum.org, supported by the US National Science Foundation.]

#### **FIGURE 2 | Contrast-to-noise ratio (CNR) of selected anatomical regions**. Mean values of CNR across six subjects and both hemispheres were extracted from R<sup>1</sup> (blue bars), R ∗ 2 (red bars), and susceptibility maps (green bars) as well as T <sup>2</sup>-weighted images (black bars). Brain regions include: RN red nucleus vs. surrounding tissue, SN - substantia nigra vs. tissue between red nucleus and substantia nigra, CTT - central tegmental tract vs.

surrounding tissue (without red nucleus), sup col - superior colliculus vs. adjacent tissue, inf col - inferior colliculus vs. adjacent tissue, RTN reticulotegmental nucleus vs. surrounding white matter tissue, MCP - middle cerebellar peduncle vs. adjacent cerebrospinal fluid, SCP - superior cerebellar peduncle vs. reticular formation, TPF - transverse pontine fibers vs. corticospinal tract. The error bars indicate the 95%-confidence interval.

were also barely visible on *T*2-weighted images (MCP in three subjects, SCP in one subject; **Figure 2**). The susceptibility maps showed the highest CNR for nearly all investigated structures within the pons (**Figure 2**). Only the CNR value of MCP with respect to adjacent cerebrospinal fluid were outperformed by the *R* ∗ <sup>2</sup> maps. Although pontine veins were only seen on both *R* ∗ 2 and susceptibility maps (**Figures 4B,C**; arrow n), their identification compared to adjacent tissue structures (e.g., transverse pontine fibers) was substantially improved on susceptibility maps. The solid white outline in **Figures 4A–D** encloses the caudal part of the pontine reticular nucleus, the facial nucleus, salivatory nucleus and nucleus abducens. Although *R*1, *R* ∗ 2 , and susceptibility maps displayed heterogeneous signal intensities in this area, discrimination of these nuclei was not possible. The interruption of the fiber tracts between the MCP seen on the track-density image (dashed outline in **Figure 4D**) is most likely caused by data acquisition and TDI post-processing issues (incompletely corrected geometric distortions in the unwarped diffusion-weighted images that impeded reliable TDI post-processing) because no such interruption is observed on the susceptibility map (**Figure 4C**) and on the histological stains (**Figures 4G,H**). The *T*2-weighted images (**Figures 3E** and **4E**) displayed rather homogeneously with low contrast (**Figures 2**–**4**), which prevented reliable identification of any substructures across all subjects.

**Figures 5** and **6** display the rostral and middle part of the medulla oblongata, respectively. The pyramid (arrow a) is located anterior to the inferior olive and was clearly visible on all contrasts (four of six subjects on *R* ∗ <sup>2</sup> maps) except the *T*2-weighted images. It could be well discriminated from the inferior olive (arrow b) which was also clearly discernible on the *R*<sup>1</sup> (five subjects) and susceptibility map (four subjects), and less obvious on *R* ∗ <sup>2</sup> maps (three subjects). The inferior cerebellar peduncle was clearly visible on track-density images and was distinguishable on *R*1, *R* ∗ 2 , and susceptibility maps in four subjects at least. Only the susceptibility maps enabled identification of the spinal trigeminal nucleus (arrow c in **Figure 5**; four subjects). Although the histological stains (**Figures 5G,H**) suggested the presence of several cranial nuclei, further nuclei were not visible on the MR images.

Finally, *R*1, *R* ∗ 2 , and susceptibility differences are listed in **Table 1** for selected anatomical regions. The differences of the susceptibility and the *R* ∗ 2 values between fibers running nearly parallel (transverse pontine fibers) and nearly perpendicular (corticospinal tract) to the static magnetic field were (−49 ± 11) ppb and (8.9 ± 6) s−<sup>1</sup> , respectively.

2 susceptibility maps as well as track-density (same color-encoding as in **Figure 1**; slightly distorted compared to the other axial image contrasts) and T <sup>2</sup>-weighted images of the same region are presented in axial orientation from **(A–E)**, respectively. The sectional plane of the axial images is indicated by the dashed line overlaid on the sagittal R<sup>1</sup> map shown in **(F)**. Axial histological sections from a different subject stained for myelin and cells are illustrated in **(G)** and **(H)**, respectively. The arrows in the images indicate: (a) middle cerebellar peduncle, (d) superior cerebellar peduncle, (e) medial longitudinal fasciculus, (h) corticospinal fibers,

## **DISCUSSION**

We have investigated the anatomy of the human brainstem *in vivo* by applying MR imaging at 7 T with high spatial, isotropic resolution, and using multiple image contrasts. QSM and TDI have been applied jointly for the first time to identify anatomical substructures of the brainstem. Both MR methods directly reflect subtle variations in tissue composition that were found to be consistent with histology as demonstrated with corresponding histological stains selected from a data base. In contrast,longitudinal relaxation rate mapping and *T*2-weighted imaging both performed inferiorly in their ability to delineate anatomical details (**Figure 2**).

*R* ∗ 2 relaxation comprises both intrinsic microscopic transverse relaxation (R2) and relaxation due to heterogeneous magnetic susceptibility distributions on a mesoscopic or macroscopic scale (*R* 0 2 ). In contrast to the *R* ∗ 2 relaxation rate, which increases with the concentration of both iron (Langkammer et al., 2010) and myelin (Lee et al., 2012), magnetic susceptibility shows a different dependence on the concentration of these substances, i.e., higher magnetic susceptibility values with increasing iron concentration (l) abducens nucleus, (m) solitary nucleus, (n) vein, (o) fiber bundle that traverses the pons in medial to lateral direction, and (p) fourth ventricle (cerebrospinal fluid). The solid white outline in **(A–D)** summarizes the caudal part of the pontine reticular nucleus, the facial nucleus, salivatory nucleus, and abducens nucleus.The dashed white outline in **(D)** indicates a region between the middle cerebellar peduncle that could not be resolved with TDI. [The histological stains **(G,H)** were adapted with permission from http://www.brains.rad.msu.edu and http://brainmuseum.org, supported by the US National Science Foundation.]

(Langkammer et al., 2012b), and lower magnetic susceptibility values with increasing myelin density (Liu et al., 2011; Schweser et al., 2011; Li et al., 2012a,b). Hence, image contrast on *R* ∗ 2 and susceptibility maps appears complementary. The longitudinal relaxation rate, *R*1, in turn is heavily influenced by myelin content, but shows little effect with respect to tissue iron concentration (Bridge and Clare, 2006).

Both *R* ∗ <sup>2</sup> maps and susceptibility maps provided detailed depiction of the non-cranial nuclei and nerve fibers in the midbrain in accordance with recent studies (Schweser et al., 2011; Deistung et al., 2013a). In particular, the red nuclei and the substantia nigra were clearly discernible on susceptibility and *R* ∗ <sup>2</sup> maps because both structures contain substantially higher iron concentrations compared to the surrounding brainstem tissue (Morris et al., 1992). Interestingly, the highly myelinated rim around the red nuclei, which was clearly seen on the myelin stain and on the susceptibility map was not visible on the *R*<sup>1</sup> map, potentially due to the lower spatial resolution of the MP2RAGE acquisition compared to the gradient-echo acquisition.

**FIGURE 5 | Images of the rostral part of the medulla oblongata**. R1, R ∗ 2 , and susceptibility maps as well as direction-encoded track-density (same color-encoding as in **Figure 1**; slightly distorted compared to the other axial image contrasts) and T <sup>2</sup>-weighted images of the same region are presented in axial orientation from **(A–E)**, respectively. The sectional plane of the axial images is indicated by the dashed line overlaid on the sagittal R<sup>1</sup> map shown in **(F)**. Axial histological sections from a different subject stained for myelin and cells are illustrated in

**(G)** and **(H)**, respectively. The arrows indicate: (a) pyramid, (b) inferior olive, (c) spinal trigeminal nucleus, (d) inferior cerebellar peduncle, (e) medial longitudinal fasciculus, (f) medial lemniscus, (h) reticular formation, (i) accessory cuneate nucleus, (j) solitary nucleus, (k) dorsal motor nucleus of the vagus nerve, and (l) hypoglossal nucleus. [The histological stains **(G,H)** were adapted with permission from http://www.brains.rad.msu.edu and http://brainmuseum.org, supported by the US National Science Foundation.]

**FIGURE 6 | Images of the middle part of the medulla oblongata**. R1, R ∗ 2 , and susceptibility maps as well as track-density (same color-encoding as in **Figure 1**; slightly distorted compared to the other axial image contrasts) and T <sup>2</sup>-weighted images of the same region are presented in axial orientation

from **(A–E)**, respectively. The sectional plane of the axial images is indicated by the dashed line overlaid on the sagittal R<sup>1</sup> map shown in **(F)**. The arrows in the images indicate: (a) pyramid, (b) inferior olive, and (d) inferior cerebellar peduncle.

Brainstem nuclei were barely seen on the *R*1, *R* ∗ 2 and susceptibility maps sectioning the pons and medulla oblongata; however, the inferior olive (**Figures 5** and **6**, arrow b) could be identified and

both the spinal trigeminal nucleus (**Figure 5C**, arrow c) and the reticulotegmental nucleus (**Figure 3C**, arrow c) were discernible on the susceptibility maps. The inferior olive is the sole source **Table 1 | Mean values of longitudinal relaxation rate (R1), effective transverse relaxation rate (R** ∗ **2 ), and volume magnetic susceptibility difference (**∆χ**) with respect to normal appearing white matter over six subjects and both brain hemispheres for selected anatomical regions**.


The values are presented as mean ± standard deviation.

of the climbing fibers to the Purkinje cells in the cerebellar cortex and is involved in learning and timing of motor behavior (De Zeeuw et al., 1998). Lesions in the inferior olive can thus introduce restrictions in performing specialized motor tasks (Martin et al., 1996). The spinal trigeminal nucleus is the largest of all cranial nerve nuclei and one of the three paired sensory nuclei associated with the trigeminal nerve. It receives input from pain, temperature, and some tactile afferents in the trigeminal nerve and is involved in migraine headache (Kaube et al., 2002). The reticulotegmental nucleus of the pons is an important component of the oculomotor circuit that regulates horizontal saccades and smooth pursuit movements of the eyes (Keller and Crandall, 1981; Büttner-Ennever and Horn, 1997). Damage to the reticulotegmental nucleus has been observed in spinocerebellar ataxia (Rüb et al., 2004) and it most probably contributes to hypometric horizontal saccades and slowing of smooth pursuits that characteristically develop in patients suffering from Alzheimer's disease (Rüb et al., 2001). Therefore, non-invasive assessment of these nuclei may help to improve diagnosis and understanding of these diseases.

These nuclei (inferior olive, spinal trigeminal nucleus, reticulotegmental nucleus) were visible on the susceptibility maps because of their different iron concentration and/or myelin content with respect to surrounding tissue. A histological study by Morris et al. (1992), for instance, demonstrated iron accumulation in the nuclei of the dorsal parts of the brainstem, the pyramids, the spinal trigeminal nucleus, and the inferior olive (although with substantially lower iron concentration than in the red nuclei and substantia nigra) and reported only little or no intrinsic iron reactivity for many other nuclei and tracts of the brainstem (e.g., solitary nucleus, hypoglossal nucleus, raphe pallidus nucleus, cuneate fasciculus, and gracile fasciculus). Brainstem nuclei have a similar, relatively low myelin density as seen from the histological myelin stains. Only the inferior olive (**Figure 5G**, arrow b) and the spinal trigeminal nucleus (**Figure 5G**, arrow c) and reticulotegmental

nucleus [see (Paxinos et al., 2012) and **Figure 3C**, arrow c] are embedded in tissue with substantially increased myelin density.

In addition to myelin's bulk diamagnetic nature, a dependence of the magnetic susceptibility on the axons' alignment with respect to the external magnetic field has been recently reported (Liu,2010; Li et al., 2012a; Schweser et al., 2012; Wharton and Bowtell, 2012). Similar direction dependent observations have been made for *R* ∗ 2 (Bender and Klose, 2010; Denk et al., 2011; Lee et al., 2011; Sati et al., 2012). We measured differences of the susceptibility and the *R* ∗ 2 values between fibers running nearly parallel (transverse pontine fibers) and nearly perpendicular (corticospinal tract) to the static magnetic field of (−49 ± 11) ppb and (8.9 ± 6) s−<sup>1</sup> , respectively. This susceptibility difference is nearly threefold as large as the values that have been derived for a hollow cylinder model of myelin (−18 to −16 ppb) (Li et al., 2012a; Wharton and Bowtell, 2012). *In vivo*, however, Li et al., 2012a also reported susceptibility differences between parallel and perpendicular fibers of up to −40 ppb in individual white matter regions,a value which is almost the same as reported here. The result for the orientation dependence of *R* ∗ 2 is in good agreement with the recently reported *ex vivo* value of −6.44 ± 0.15 s−<sup>1</sup> at 7 T (Lee et al., 2011). The orientationdependent differences of magnetic susceptibility and *R* ∗ 2 in our study may be caused by the presence of additional contributors other than myelin, such as pontine nuclei and neuropil within the investigated volumes-of-interest. Nevertheless, consistent with the theory of myelin's orientation dependency, the contrast between fibers running in medial-lateral direction (oriented nearly perpendicular to the magnetic field; e.g., transverse pontine fibers, arrow o in **Figure 4**) and fibers running inferior-superiorly (oriented nearly parallel to the magnetic field; corticospinal tract) is enhanced on the *R* ∗ 2 and susceptibility maps. More sophisticated approaches, such as susceptibility tensor imaging (STI) (Liu, 2010; Li et al., 2012a) or multipole susceptibility tensor imaging (MSTI) (Liu and Li, 2013) may further improve discrimination between fibers and nuclei of the pons and medulla, particularly in combination with high spatial resolution.

Direction-encoded TDI enabled discrimination of the directionality of nerve fibers within the brainstem and yielded complementary information to the relaxation and susceptibility maps. The depiction of larger fiber bundles was consistent with the results of other DTI studies conducted *in vivo* at lower magnetic fields (Nagae-Poetscher et al., 2004; Salamon et al., 2005). We were, however, unable to resolve unambiguously the transverse pontine fibers, which were clearly visible on the *R* ∗ 2 and susceptibility maps (**Figure 4**). This may be due to the lower base voxel resolution of 1.2 mm of the DWI data, which was chosen as a compromise between spatial resolution, SNR, and acquisition time. To overcome this limitation we interpolated the DWI data using a super-resolution approach (TDI) to obtain an isotropic virtual resolution of 0.43 mm. It should be noted, however, that the TDI technique is not able to fully recover the whole information that would be present if the data would be acquired with such a high spatial resolution (Calamante et al., 2011). Due to the larger errors in tractography by performing fiber tracking on data with a lower resolution, the higher resolved track-density images appeared slightly blurred and lost some fine details; explaining the difficulties in identifying the transverse pontine fibers in the presented track-density images (**Figures 3D** and **4D**). As might be expected, the level of blurring in the track-density images depends on the desired virtual resolution (Calamante et al., 2011). Higher spatial resolution is thus needed to delineate transverse pontine fibers unambiguously, which requires sophisticated modifications of the diffusion MRI sequence such as combining zoomed imaging and parallel imaging (Heidemann et al., 2012; Eichner et al., 2013). Karampinos et al. (2009) were recently able to resolve transverse pontine fibers in direction-encoded fractional anisotropy maps with a resolution of 0.8 mm × 0.8 mm × 3 mm at 3T by applying a dedicated acquisition scheme (self-navigated, multi-shot, variable-density, spiral-imaging with outer volume suppression).

One particular problem of brainstem imaging with MRI are the induced field inhomogeneities caused by nearby air-tissue or bonetissue interfaces. In the vicinity of such interfaces significant image distortions of EPI-based diffusion-weighted images occur and additional spin dephasing arises in gradient-echo images. Despite unwarping the diffusion-weighted images based on gradient-echo field-map information, it was not possible to completely compensate geometric distortions, thus impeding accurate depiction of ventral fiber structures of the middle pons and medulla oblongata on track-density images. These induced field inhomogeneities manifest themselves as steep gradients in gradient-echo phase images and lead to signal reductions in gradient-echo magnitude images and decreased SNRs, affecting both susceptibility maps [due to increased phase noise (Gudbjartsson and Patz, 1995)] and *R* ∗ <sup>2</sup> maps. In future studies, the impact of spin dephasing due to airtissue or bone-tissue interfaces on *R* ∗ <sup>2</sup> maps may be minimized by taking into account the macroscopic field inhomogeneity reflected in the phase images (de Leeuw and Bakker, 2012). *R*<sup>1</sup> maps and*T*2 weighted images, on the other hand, are distinctly more immune against these field inhomogeneity effects.

The quantitative analysis revealed similar values for the longitudinal relaxation rates (*R*1) for nuclei and myelinated fibers (**Table 1**), whereas susceptibility and *R* ∗ 2 values varied more substantially. This broader spread of *R* ∗ 2 and susceptibility values

## **REFERENCES**


concentration in normal aging using quantitative susceptibility mapping. *Neuroimage* 59, 2625–2635. doi:10.1016/j.neuroimage.2011.08. 077


is reflected in the large CNR values (**Figure 2**) underlining the good discrimination of nerve fibers on susceptibility maps and the exquisite delineation of midbrain nuclei on susceptibility and *R* ∗ <sup>2</sup> maps. It should, however, be noted that the discrimination of brain tissue with respect to cerebral spinal fluid (CSF) on susceptibility maps is not as striking as on the *R*1, *R* ∗ 2 , and *T*2-weighted images [**Table 1**; **Figure 4** (arrow p)]. Due to the quantitative nature of *R*1,*R* ∗ 2 , and magnetic susceptibility maps, these contrasts may be combined by projections along support vectors observed from discriminant analysis or support vector machines to generate composite images with improved depiction of anatomical features while providing improved discrimination of CSF (Deistung et al., 2013b).

In conclusion, maps of magnetic susceptibility displayed the largest number of brainstem structures, including larger and smaller fiber pathways as well as several nerve nuclei. Usage of multiple image contrasts enables a detailed non-invasive view into tissue structure and composition. Hence, multi-contrast MR imaging of the brainstem at ultra-high magnetic fields that utilizes relaxation, diffusion, and magnetic susceptibility information is a versatile tool to assess anatomy individually in great detail.

## **ACKNOWLEDGMENTS**

We would like to thank Andrew Webb (Gorter Center for High Field MRI, Department of Radiology, Leiden University Medical Center) for providing the dielectric pads as well as Penny Gowland and Andrew Peters (Sir Peter Mansfield Magnetic Resonance Centre, School of Physics and Astronomy, University of Nottingham) for providing the tailored RF pulse. This research was funded by a German Research Foundation research grant (DFG, RE 1123/9-2), by the Max-Planck society, and by seed grants awarded by the Interdisciplinary Center for Clinical Research (IZKF) in Jena (to Andreas Deistung), by the International Society of Magnetic Resonance in Medicine (to Ferdinand Schweser), and by the Friedrich Schiller University Jena (to Ferdinand Schweser).


Koekkoek, S. K., and Ruigrok, T. J. (1998). Microcircuitry and function of the inferior olive. *Trends Neurosci.* 21, 391–400. doi:10.1016/ S0166-2236(98)01310-1


*Magn. Res. Med.* 34, 910–914. doi:10.1002/mrm.1910340618


*Neuroimage* 59, 1413–1419. doi:10. 1016/j.neuroimage.2011.08.045


of brain dysmyelination by quantitative mapping of magnetic susceptibility. *Neuroimage* 56, 930–938. doi:10.1016/j.neuroimage.2011.02. 024


for deep brain stimulation using a standard installation protocol at 3.0 Tesla. *Acta Neurochir. (Wien)* 154, 481–494. doi:10.1007/s00701- 011-1242-8


tensor imaging. *Neuroradiology* 47, 895–902. doi:10.1007/s00234-005- 1439-8


*Caused by Disease*. Chichester, UK: John Wiley & Sons, Ltd. doi:10.1002/ 0470869526


**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: 30 May 2013; paper pending published: 25 June 2013; accepted: 07 October 2013; published online: 29 October 2013.*

*Citation: Deistung A, Schäfer A, Schweser F, Biedermann U, Güllmar D, Trampel R, Turner R and Reichenbach JR (2013) High-resolution MR imaging of the human brainstem in vivo at 7 Tesla. Front. Hum. Neurosci. 7:710. doi: 10.3389/fnhum.2013.00710*

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

*Copyright © 2013 Deistung , Schäfer, Schweser, Biedermann, Güllmar, Trampel, Turner and Reichenbach. 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.*

## Characterizing aging in the human brainstem using quantitative multimodal MRI analysis

#### **Christian Lambert 1,2\*, Rumana Chowdhury <sup>2</sup> ,Thomas H. B. FitzGerald<sup>2</sup> , Stephen M. Fleming3,4 , Antoine Lutti 2,5, Chloe Hutton<sup>2</sup> , Bogdan Draganski <sup>5</sup> , Richard Frackowiak <sup>5</sup> and John Ashburner <sup>2</sup>**

<sup>1</sup> Clinical Neuroscience, St George's University of London, London, UK

<sup>2</sup> Wellcome Trust Centre for Neuroimaging, UCL Institute of Neurology, University College London, London, UK

<sup>3</sup> Center for Neural Science, New York University, New York, NY, USA

<sup>4</sup> Department of Experimental Psychology, University of Oxford, Oxford, UK

<sup>5</sup> LREN, Département des Neurosciences Cliniques, CHUV, Université de Lausanne, Lausanne, Switzerland

## **Edited by:**

Simon Baudrexel, Goethe University Frankfurt, Germany

### **Reviewed by:**

Pierre-Louis Bazin, Max Planck Institute for Human Cognitive and Brain Sciences, Germany Berkin Bilgic, Martinos Center for Biomedical Imaging, USA Lijuan Zhang, SIAT; Chinese Academy of Sciences, China

## **\*Correspondence:**

Christian Lambert, St George's University of London, Cranmer Terrace, SW17 0RE, London, UK e-mail: clambert112358@gmail.com Aging is ubiquitous to the human condition.The MRI correlates of healthy aging have been extensively investigated using a range of modalities, including volumetric MRI, quantitative MRI (qMRI), and diffusion tensor imaging. Despite this, the reported brainstem related changes remain sparse.This is, in part, due to the technical and methodological limitations in quantitatively assessing and statistically analyzing this region. By utilizing a new method of brainstem segmentation, a large cohort of 100 healthy adults were assessed in this study for the effects of aging within the human brainstem in vivo. Using qMRI, tensor-based morphometry (TBM), and voxel-based quantification (VBQ), the volumetric and quantitative changes across healthy adults between 19 and 75 years were characterized. In addition to the increased R2<sup>∗</sup> in substantia nigra corresponding to increasing iron deposition with age, several novel findings were reported in the current study. These include selective volumetric loss of the brachium conjunctivum, with a corresponding decrease in magnetization transfer and increase in proton density (PD), accounting for the previously described "midbrain shrinkage." Additionally, we found increases in R1 and PD in several pontine and medullary structures. We consider these changes in the context of well-characterized, functional age-related changes, and propose potential biophysical mechanisms. This study provides detailed quantitative analysis of the internal architecture of the brainstem and provides a baseline for further studies of neurodegenerative diseases that are characterized by early, pre-clinical involvement of the brainstem, such as Parkinson's and Alzheimer's diseases.

**Keywords: brainstem, quantitative MRI, tensor-based morphometry, voxel-based quantification, aging**

## **INTRODUCTION**

Aging is ubiquitous to the human condition. Previously, the cortical and subcortical correlates of normal aging have been extensively investigated using several modalities including volumetric MRI (Jernigan et al., 2001;Walhovd et al., 2005, 2011), quantitative MRI (qMRI) (Armstrong et al., 2004; Draganski et al., 2011; Bilgic et al., 2012), cortical thickness measures (Thambisetty et al., 2010), and diffusion tensor imaging (DTI) (Pfefferbaum et al., 2010; Vaillancourt et al., 2012). These studies have established that the normal aging process is associated with distinctive morphological changes including volumetric loss of the cortex, predominantly in the prefrontal, parietal, and temporal regions in addition to the amygdala, hippocampus, striatum, and cerebellum (Woodruff-Pak et al., 2010; Draganski et al., 2011; Walhovd et al., 2011). Corresponding with this, there is a reduction in cortical white-matter myelination, reflected by falling magnetization transfer (MT) values, and increasing iron deposition particularly in basal ganglia structures, notably the substantia nigra, as reflected by increasing R2<sup>∗</sup> values (Pfefferbaum et al., 2009, 2010; Haacke et al., 2010; Draganski et al., 2011; Bilgic et al., 2012). Despite these widespread

effects, reported brainstem related changes remain sparse due to technical limitations of imaging, segmenting, and statistically analyzing data from this region (Luft et al., 1999; Raz et al., 2001; Lee et al., 2009). Aging is known to negatively impact on several brainstem-mediated functions, for example the sleep-wake cycle (Hut and Van der Zee, 2011), sympathetic outflow (Samuels and Szabadi, 2008), vestibular-ocular reflexes (Baloh et al., 1993), and cardiovascular reflexes (Vita et al., 1986). Post-mortem reports indicate that changes within the brainstem sub-nuclei do take place during aging (Alvarez et al., 2000; Samuels and Szabadi, 2008), some of which may be early precursors for subclinical neurodegenerative disease (Tsopelas et al., 2011).

Quantitative MRI produces quantitative MR parameters that can be used as biomarkers of tissue microstructure (Tofts, 2003; Draganski et al., 2011). Examples of these parameters include MT, Proton Density (PD), R2<sup>∗</sup> , and R1. MT emerges from hydrogen in motionally restricted macromolecules and more directly relates to macromolecular content. Biologically, it is a reflection of the quantity of myelin within a voxel (Helms et al., 2008b). It is important to note that in the literature MT contrast is often reported as

"MT ratio" (Dousset et al., 1992), a parameter that shows residual T1 dependence. However, by using semi-quantitative parameter, the MT saturation, MT and T1 effect can be separated (Helms et al., 2008b). PD refers to the concentration of MRI-visible water (Tofts, 2003). R2<sup>∗</sup> (=1/T2<sup>∗</sup> ) is the relaxation rate of the transverse magnetization, and is linearly correlated with tissue iron concentration (Yao et al., 2009). Finally R1 (=1/T1) is the longitudinal relaxation rate, and arises from a mix of water content, iron, and tissue macromolecule fraction (Rooney et al., 2007). Increases in iron (Rooney et al., 2007), decreases in PD (Gelman et al., 2001), and increases in lipid content (Stanisz et al., 2005) all cause an increase in the measured R1 signal.

The aim of this cross-sectional study was to characterize volumetric and tissue parameter changes associated with aging within the human brainstem *in vivo* using qMRI, tensor-based morphometry (TBM), and voxel-based quantification (VBQ). This was motivated by the observation that certain neurodegenerative diseases, such as Parkinson's and Alzheimer's diseases, are characterized by early pre-clinical involvement of the brainstem (Simic et al., 2009; Hawkes et al., 2010). In order to further study these effects, changes associated with the normal aging process must first be better characterized.

## **MATERIALS AND METHODS**

## **SUBJECTS**

Multiparametric maps were acquired from previous studies (FitzGerald et al., 2012; Lambert et al., 2012; Chowdhury et al., 2013). In total, imaging data for 100 healthy adults (47 males, aged 40.3 ± 21.2 years; 53 females, aged 48.2 ± 22.7 years; age ranged 19–75 years) who had MRI scanning at the Wellcome Trust Centre for Neuroimaging was used. Due to the selection criteria, these predominately fell into a bimodal distribution (as is shown **Figure 5**), with a younger cohort less than 60 years [*n* = 58, mean age 25.8 years (SD 7.6 years)] and an older cohort above 60 years [*n* = 42, mean age 69.1 years (SD 3.5 years)]. Subjects above 60 years had a normal neurological examination performed by a physician (Rumana Chowdhury), MMSE score >28 and a normal performance (within 1.5 SD of age-related norm) on a range of neuropsychological tests. Involvement of human volunteers was approved by the local ethics committee, and each subject provided written informed consent prior to MRI examination.

## **IMAGE ACQUISITION**

MR imaging was performed on a 3T whole-body MRI system (Magnetom TIM Trio, Siemens Healthcare, Erlangen) operated with a whole-body transmit radio-frequency (RF) coil and a 32 channel RF receive coil. MR data of 21 subjects was acquired on a second identical MRI system located within the same department. Each participant underwent a multiparametric mapping (MPM) scanning protocol for quantitative mapping of multiple MR parameters. MT, R1 (=1/T1), and PD weighted images were acquired using 3D multi-echo FLASH (fast low-angle shot) acquisitions (Helms et al., 2008a,b). Full imaging parameters are summarized in **Table 1**. The image resolution was 1 mm isotropic. The total acquisition time was 19 min. For each subject quantitative MT, R1, PD, and R2<sup>∗</sup> (=1/T2<sup>∗</sup> ) maps were extracted from the acquired images using in-house MATLAB program. An additional dataset

was acquired on each subject for mapping of the RF transmit field B1+ over the brain (4 mm isotropic resolution, acquisition time 3 min) using the 3D EPI SE/STE method described in (Lutti et al., 2010, 2012). A B0-field map was also acquired to correct for the distortions of the EPI images acquired for the B1+ mapping acquisition (acquisition time 2 min). The resulting B1+ maps were used to correct the MPM maps for RF transmit field inhomogeneity effects (Helms et al., 2008a). The PD maps were corrected for the spatially varying sensitivity profile of the receive coil using the UNICORT algorithm (Weiskopf et al., 2011). The resulting "flattened" signal amplitude maps were converted into PD maps by scaling the voxels by the expected average PD of white-matter [69% (Tofts, 2003)].

## **PRE-PROCESSING**

The MT maps were initially segmented into gray, white, and CSF tissue classes whilst maintaining the native resolution (1 mm) by using the unified segmentation within SPM8 (http://www. fil.ion.ucl.ac.uk/spm/) (Ashburner and Friston, 2005). For each individual, the total intracranial volume (TIV) was calculated by summing the gray matter, white-matter, and CSF segmentations at a threshold of 0.1. The segmentations were then registered diffeomorphically to a common group-average 1 mm isotropic template (using the "*Shoot toolbox*" of SPM, Ashburner and Friston, 2011) to produce deformation fields that were subsequently used in the segmentation step to warp the brainstem tissue priors to individual subject space.

A method for generating brainstem specific tissue probability maps (TPMs) and subsequent segmentation had been previously developed (Lambert et al., 2013). In brief, the brainstem TPMs were defined using a modified multivariate Mixture of Gaussians (mmMoG) to generate spatial TPMs from 0.8 mm isotropic MT and PD maps, which were masked to include the brainstem from the origin of the cerebral aqueduct to the level of the foramen magnum. Four tissue classes were identified, three gray matter and one white-matter. These were labeled for descriptive purposes according to the predominant tissue type present: tissue class one included the substantia nigra, locus coeruleus and raphe nuclei and hence was designated "*monoaminergic gray matter*," though regions that contained dorsal cranial nerve nuclei were also included. Tissue class two consisted mainly of the nucleus reticularis throughout its length and the pontine nuclei, and hence was designated "*reticulated gray matter*." Tissue class three was specific for the periaqueductal gray (PAG) and labeled as such. Tissue class four was the brainstem white-matter. Examples of these tissue classes projected onto the corresponding sections from Duvernoy's 9.4T MRI brainstem atlas (Naidich and Duvernoy, 2009) have been provided in **Figure 1** (reprinted from Lambert et al., 2013 with permission from Elsevier).

The previously calculated tissue priors for four brainstem and one non-brainstem tissue classes were aligned with the calculated group-average template. This was achieved by realigning a whole brain template in the same space as the TPMs with the corresponding new group-average template, re-slicing to achieve 1 mm isotropic voxel size with a 4th degree spline interpolation and ensuring each individual voxel probability value summed to one across all the maps.

## **Table 1 | Imaging parameters.**


The actual brainstem segmentation step is summarized in **Figure 2**. The realigned probability maps were warped to individual subject space and used in SPM8 "*New Segment* " on individual MT and PD maps that had been cropped with a set-bounding box to include only the brainstem region. Specifically, five tissue classes were used; four within brainstem and one for everything else. Two Gaussians were used to model each tissue class except for the PAG matter (one Gaussian) and non-brainstem (eight Gaussians). Individual level segmentations were generated from the MT and PD images. All the images from each of the four tissue classes were then cropped using a common bounding box to increase computational speed and visually checked to ensure good quality. All of the brainstem tissue classes were then re-warped back to a group-average using a diffeomorphic warping algorithm (geodesic shooting) to allow re-estimation of the Jacobian determinants. The deformation fields were then used to warp each parametric map (MT, PD, R1, R2<sup>∗</sup> ) to the common template. The resulting images were masked using the non-brainstem tissue class as an exclusion mask.

## **ANATOMICAL ANALYSIS**

This work utilized two techniques to quantitatively analyze qMRI parameter maps. The first was TBM, which is a method to characterize volumetric change *in vivo*. It is similar to the already well-described voxel-based morphometry (Ashburner and Friston, 2000), however instead of statistically analyzing warped modulated tissue segmentations, the Jacobian determinant images are used (Ashburner and Friston, 2001). The second was VBQ. This is a pipeline to allow unbiased mass univariate statistical analysis of quantitative parameter maps whilst controlling for multiple comparisons with family-wise error (Draganski et al., 2011). The original VBQ method (Draganski et al., 2011) was modified slightly for this study. Because the brainstem lacks gyrification, and due to the highly accurate warping algorithm used, smoothing of the warped quantitative data was avoided, hence negating the need to produce "*warped-weighted average*" images that were required for the published approach. This had several advantages. First, only one statistical analysis was required for each

quantitative map (i.e., the warped quantitative map) rather than one for each tissue class (i.e., four), which would have been necessary with the warped-weighted average approach. Additionally, because analysis was carried out directly on the warped quantitative maps without smoothing, a higher degree of spatial accuracy could be achieved.

## **STATISTICAL ANALYSIS**

All statistical analysis was carried out using SPM8 in MAT-LAB 2010b. A two-sample *t*-test was initially used to check for scanner-associated differences in the acquisitions controlling for age, sex, and TIV, at FWE < 0.05 correction for multiple comparisons. No voxels survived correction. A design matrix for multiple linear regression model was then constructed including age, sex, and TIV as covariates. The TIV was centered around the mean, and the remaining covariates remained uncentered. An intercept was included in the model but no normalization was used. Using this design matrix, TBM was performed by analyzing the Jacobian determinant images (Ashburner and Friston, 2001), and VBQ by analyzing each warped quantitative map. For each analysis, the non-brainstem tissue class was used as an exclusion mask, ensuring only voxels that contained brainstem tissue were included. Each image type was assessed in a single design matrix to negate potential problems associated with uneven variances across the different quantitative maps. Each map was assessed for significant positive and negative correlates with age. For each contrast, FWE < 0.05 was reported. Finally, for each VBQ analysis where results were significant at FWE < 0.05, the T-maps were binarized for visualization at *p* < 0.001 uncorrected. These were used to create three-dimensional renderings, and also to assess the overlap of significant results between different modalities to better understand the interaction between the different measures. Finally, the same binarized images were used to extract the mean quantitative value for each individual. These values were plotted against the corresponding age, and the best linear fit indicated.

## **RESULTS**

For anatomical reference, a high resolution (0.5 mm isotropic) *ex vivo* combined MT T2<sup>∗</sup> MRI with anatomical annotations has been provided in **Figure 3** (figure adapted from Lambert et al., 2013 with permission from Elsevier).

## **TENSOR-BASED MORPHOMETRY**

Highly localized volumetric decreases in tissue volume were found symmetrically within the brachium conjunctivum (superior cerebellar peduncle) bilaterally. These volumetric decreases were significant at FWE < 0.05 and are summarized in **Figures 4** and **5**. There were no significant positive TBM correlates with age at FWE < 0.05 and *p* < 0.001.

## **VOXEL-BASED QUANTIFICATION ANALYSIS Negative correlates with age**

Significant negative trends with age were only observed across the MT maps within the brachium conjunctivum at FWE < 0.05. These are summarized in **Figures 4** and **5**.

## **Positive correlates with age**

As shown in **Figures 5**–**7**, there were more widespread significant increases in qMRI values with age. **Figure 6** shows the individual parameter maps increases at six axial slices through the brainstem. **Figure 7** summarizes the regional MPM increases and examines the overlap between the different parameter maps. It also demonstrates the relative proportions of each tissue type that overlap with one another.

*R1 increases with age.* Increases in R1 intensity values with age were found throughout the brainstem, but were confined to gray matter structures. These were found in bilateral inferior olivary nuclei, pontine nuclei, dorsal raphe nucleus, and bilateral substantia nigra. To better understand this finding, a sub analysis on these significant regions was performed, looking at the remaining parameters. Whilst there may be some selection bias with this

**FIGURE 3 | High resolution ex vivo combined MTT2\* MRI with anatomical annotations for reference**. Figure adapted from Lambert et al. (2013) with permission from Elsevier. Abbreviations: A8, dopaminergic center (approximate location), CP, cerebral peduncle (anterior to posterior: consisting of frontopontine, corticonuclear, corticospinal, and parietotemporal pontine

tracts); CST, corticospinal tract; CTT, central tegmental tract; ICP, inferior cerebellar peduncle; ML, medial lemniscus; MLF, medial longitudinal fasciculus; PAG, periaqueductal gray; SCT, spinocerebellar tract; SNpc, substantia nigra pars compacta; SNpr, substantia nigra pars reticulata; TST, tectospinal tract; VTA, ventral tegmental area. \*Artifact due to fixation.

approach (Vul et al., 2009), the correlation analysis was strictly performed to compare the behavior of the different MRI modalities within the regions of R1 change, to better understand what may be contributing to this R1 increase. A significant positive correlation with R2<sup>∗</sup> was found within these regions, even when the iron rich midbrain structures were excluded (Pearson's Correlation Coefficient within pons and medulla = 0.51, *p* < 1 × 10−<sup>8</sup> ), implying that it is increasing iron deposition that significantly contributes to the observed R1 signal increases within brainstem gray matter structures.

*R2*<sup>∗</sup> *increases with age.* Increases in R2<sup>∗</sup> intensity values were localized to midbrain structures. Specifically, the substantia nigra, ventral tegmental area, and red nuclei. The pontine and medullary R1 regions described above did not show increases in R2s when analyzing the entire brainstem. This discrepancy would be due to correction for multiple comparisons at the whole brainstem level compared to a small volume region of interest.

*Proton density increases with age.* Significant increases in PD were observed in the brainstem white-matter. This included (from rostral to caudal): *Fasciculus* cerebellothalamicus, brachium conjunctivum, corticospinal tract, superior cerebellar peduncle, and medial longitudinal fasciculus. Additionally, the inferior portion of the substantia nigra also exhibited increased PD.

*MT increases with age.* Only one small region demonstrated MT increases: the corticobulbar portion of the cerebral peduncle bilaterally.

## **DISCUSSION**

In this study, we have demonstrated the age effect of MT, R1, R2<sup>∗</sup> , and PD in human brainstem using an automated, nonbiased approach that is able to resolve the internal structure of the brainstem in unprecedented detail at 3T.

## **MIDBRAIN ATROPHY IN HEALTHY AGING**

Previous studies examining age-related volumetric decline in the brainstem have found no overall volume loss (Raz et al., 2001; Lee et al., 2009) but significant midbrain atrophy (Luft et al., 1999). This had previously been attributed to shrinkage of the substantia nigra (Raz, 1996), however there is sparse histopathological evidence to support this hypothesis (McCormack et al., 2004; Collier et al., 2007). Our study agrees with previous work in that brainstem aging-associated atrophy seems to be confined to the midbrain. However, our results indicate that it is predominantly volume loss of the superior cerebellar fiber bundles (brachium conjunctivum, fasciculus cerebellothalamicus (Haroian et al., 1981)] that are responsible for this finding. Additionally the decreasing myelin content, reflected by decreasing MT and increasing PD, indicates this is a regional effect due to axonal loss rather than an artifact. What is striking is the regional specificity of these findings. As with previous studies (Draganski et al., 2011), we found a marked increase in iron concentration with age in the substantia nigra and red nucleus. Though this increase is largely sequestered in neuromelanin *in vivo*, this substance will be released by dying neurons and hence can contribute to the regional damage observed (Chiueh, 2001; Papanikolaou and Pantopoulos, 2005). Not only is iron directly toxic to axons by causing rapid lipid peroxidation, it can also induce neurotoxic microglial factors to be released locally that will potentiate the regional insult (Zecca et al., 2004). These findings correlate with previously described cerebellar atrophy (Andersen et al., 2003; Draganski et al., 2011; Walhovd et al., 2011) and the resultant cortical and subcortical disconnection (Taniwaki et al., 2007; Alalade et al., 2011) hence would better account for these observations.

## **BRAINSTEM GRAY MATTER CHANGES**

Our results agree with previous studies in the observation that the pons and medulla do not show volumetric loss in aging (Sullivan et al., 2004; Lee et al., 2009). Despite this, we identified widespread and specific changes in the qMRI maps within the gray matter structures of these regions. Specifically, increasing R1 signal was observed within the pontine nuclei and inferior olive, in addition to the substantia nigra and dorsal raphe nuclei within the mesencephalon. This increase in R1 signal was significantly correlated with R2<sup>∗</sup> , suggesting that increasing iron content within these structures can at least partly account for the observed gray matter changes. Within the pons, as with elsewhere in the brain, iron is predominately found in oligodendrocytes and to a lesser extent in the microglia and astrocytes (Ozawa et al., 1994). However the amount of iron in the former remains constant throughout aging and instead increases are found within the microglia and astrocytes (Zecca et al., 2004), which may relate to changes in vascular permeability. It has been speculated that the accumulation of immunoreactive iron in the microglia may cause or predispose to a neuroinflammatory response, as seen in Parkinson's and Alzheimer's disease (Zecca et al., 2004).

Chemical exchange brings a significant contribution to the transfer of magnetization between free water and macromolecules (Henkelman et al., 1993; Kucharczyk et al., 1994). Evidence has been presented which also suggests a significant impact of chemical exchange on the resonance frequency of water protons

(Shmueli et al., 2011; Wharton and Bowtell, 2012). Physiological factors that impact chemical exchange such as tissue oxygenation (O<sup>2</sup> and CO<sup>2</sup> exchange), temperature, and pH (Kucharczyk et al., 1994; Liepinsh and Otting, 1996; van Zijl et al., 2003) may alter the MT and R2<sup>∗</sup> estimates presented here. However only small variations of MT with pH have previously been reported within

biologically plausible values (Kucharczyk et al., 1994). T1 relaxation is primarily driven by dipolar interactions between water and macromolecular protons (Koenig, 1995). Additionally the small variations in MT and T1 values of brain tissue with temperature (Lewa and Majewska, 1980; Graham et al., 1999) are unlikely to have a significant impact on our results within a range of biologically reasonable temperatures. Decreases in cerebral blood flow (CBF) that vary from region to region have been reported in normal aging (Aanerud et al., 2012). Modulated by the physiological parameters listed above, CBF may have an impact on the MPM measurements. However, the precise interaction between these remains poorly understood (Zauner and Muizelaar, 1997) and their impact on aging sparsely characterized. Several findings in this study suggest that the impact of CBF on healthy aging is minimal in the brainstem. First, biologically plausible concordance was accomplished between different parameter maps, for example the decreasing MT with increasing PD in areas of axonal and volume loss. Second, the results demonstrate respect for known anatomical boundaries such as the brachium conjunctivum with increased PD. Finally, our results are in agreement with previous histological

observations, such as the increases in iron content in the substantia nigra (Zecca et al., 2001). Future work is expected to truly clarify and disambiguate the effects of CBF on MPM measurements. Not only will this improve the biophysical interpretation of these sequences in healthy tissue, but also allow better understanding of cerebral pathology where the vascular permeability will also change (Mooradian, 1988).

## **LIMITATIONS**

There are several limitations with our study. First, to accrue 100 normal control MPMs, we utilized scans that were acquired through previous studies (FitzGerald et al., 2012; Lambert et al., 2012; Chowdhury et al., 2013). There were few subjects between the ages of 35–65. Whilst this does not invalidate the findings, it makes it impossible to better characterize the temporal characteristics of the changes i.e., are they linear or non-linear. Additionally, as the upper age limit is 75, it is unclear how these changes extrapolate to those over that age. These features also bias the study toward those who are ambulatory, independent, and self

motivating. Whilst it could also be argued that the latter criticism ensures that the experimental group consists only of healthy normal controls, it must also be acknowledged that declining mobility is a recognized feature of normal aging that may be unrelated to cortical changes, and this would certainly represent a selection bias in this work. Finally, it is currently unclear how these changes map to individual function, so further work is required to better understand this aspect.

In this work, the impact of physiological noise on image quality and the segmentation results was not explored. However, recent work using these techniques in the cortex highlights its robustness and also sensitivity to tissue microarchitecture (Dick et al., 2012; Sereno et al., 2013). Physiological noise has mostly been addressed in the context of fMRI, where image stability is paramount (Glover et al., 2000; Hutton et al., 2011). These methods cannot be directly implemented in anatomical imaging due to the different type of image acquisition. Potentially beneficial techniques include phasenavigator correction methods (Hu and Kim, 1994; Barry et al., 2008) although they may reduce the efficiency of the FLASH acquisitions. Alternatively real-time shimming methods for correction of respiratory-induced effects (Van Gelderen et al., 2007) or optical systems for fast prospective correction of subject motion (Zaitsev et al., 2006) may yield a significant reduction of physiological effects on anatomical scans.

Regarding the scanning parameters, 1 mm isotopic volumes are still reasonably large for certain brainstem structures, so our ability to fully characterize the changes are likewise limited. However, many structures in the brainstem are well above the 1 mm<sup>3</sup> threshold such as the facial nerve nucleus (mean volume = 12.95 mm<sup>3</sup> ) and hypoglossal nerve nucleus (mean volume = 14.39 mm<sup>3</sup> ) (Sherwood et al., 2005), and hence this current study is of sufficient resolution for these nuclei. Additionally, we have demonstrated that widespread and structure specific changes can be found that also correspond to known features of brainstem aging.

## **FUTURE APPLICATIONS**

This work provided a baseline of qMRI changes with aging in the brainstem. Further work is required to characterize the exact temporal dynamics of these changes, and how they extrapolate to those beyond the age of 75. Combining these techniques with diffusion-weighted imaging would allow further characterization of the changes within the white-matter and connectivity properties, and could be used to further refine the segmentations. Future work will examine these parameters in the context of

## **REFERENCES**


the Parkinsonian disorders, such as idiopathic Parkinson's disease, progressive supranuclear palsy, and multisystem atrophy. These conditions are good models of brainstem disease with which to develop the imaging techniques, as they all show postmortem histological changes within the brainstem that, to date, are only evident on MRI at advanced stages of disease. Additionally, further longitudinal work is required to understand how individuals whose qMRI parameters lay well outside their expected normal age-matched distributions develop over time, to identify whether these parameters can serve as early biomarkers of neurodegenerative disease, an important prerequisite for any future neuroprotective therapy.

## **CONCLUSION**

In conclusion, we have characterized changes in the brainstem due to normal healthy aging using qMRI and volumetric analysis. We replicate previous findings of midbrain shrinkage, and put forward a new hypothesis as to the underlying mechanism based on statistically significant regional changes in the qMRI maps. Specifically, axonal loss in the ascending cerebellar fiber bundles is responsible for the decreased brainstem volume loss rather than atrophy of the substantia nigra, as previously speculated. Additionally, it may be that the regional increases in nigral and rubral iron content, as reflected by R2<sup>∗</sup> signal, underpin these observations. Finally we demonstrate widespread brainstem changes during aging evidenced by increasing PD in the white-matter, and increasing R1 in the gray matter, the latter being significantly driven by increasing iron deposition. This work provides a baseline from which brainstem pathology can be better explored *in vivo* using 3T MRI, and non-invasive biomarkers of different neurodegenerative conditions.

## **ACKNOWLEDGMENTS**

This work was supported by Wellcome Trust Grant 075696/Z/04/Z (R.S.J.F., Sarah Tabrizi, John Ashburner). The Wellcome Trust Centre for Neuroimaging is supported by core funding from the Wellcome Trust 091593/Z/10/Z. Christian Lambert is supported by the Academy of Medical Sciences Starter Grants for Clinical Lecturers. Bogdan Draganski is supported by the Swiss National Science Foundation (320030\_135679 and NCCR Synapsy), Foundation Parkinson Switzerland, Deutsche Forschungsgemeinschaft (KFO 247/0), Novartis Foundation for medical-biological research and the Synapsis Foundation. We thank all participants in our study and the radiographers at the Functional Imaging Laboratory for their assistance acquiring data.

Tallent, E. M., et al. (2004). Agerelated, regional, hemispheric, and medial-lateral differences in myelin integrity in vivo in the normal adult brain. *AJNR Am. J. Neuroradiol.* 25, 977–984.


*J. Neurosci.* 32, 16417–16423. doi:10. 1523/JNEUROSCI.3254-12.2012


spin-lattice relaxation time T1 in biological tissues. *Bull. Cancer* 67, 525–530.


brainstem? *Neuropathol. Appl. Neurobiol.* 35, 532–554. doi:10. 1111/j.1365-2990.2009.01038.x


*Magn. Reson. Med.* 49, 440–449. doi: 10.1002/mrm.10398


motion correction using an external optical motion tracking system. *Neuroimage* 31, 1038–1050. doi:10. 1016/j.neuroimage.2006.01.039


**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: 02 May 2013; accepted: 25 July 2013; published online: 20 August 2013. Citation: Lambert C, Chowdhury R, FitzGerald THB, Fleming SM, Lutti A, Hutton C, Draganski B, Frackowiak R and Ashburner J (2013) Characterizing aging in the human brainstem using quantitative multimodal MRI analysis. Front. Hum. Neurosci. 7:462. doi: 10.3389/fnhum.2013.00462*

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

*Copyright © 2013 Lambert, Chowdhury, FitzGerald, Fleming , Lutti, Hutton, Draganski, Frackowiak and Ashburner. 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.*

## Change in brainstem gray matter concentration following a mindfulness-based intervention is correlated with improvement in psychological well-being

#### **Omar Singleton<sup>1</sup> , Britta K. Hölzel <sup>2</sup> , Mark Vangel <sup>1</sup> , Narayan Brach<sup>3</sup> , James Carmody <sup>4</sup> and SaraW. Lazar <sup>1</sup>\***

<sup>1</sup> Massachusetts General Hospital, Harvard Medical School, Boston, MA, USA

2 Institute for Medical Psychology, Charité Universitätsmedizin, Berlin, Germany

<sup>3</sup> PGSP-Stanford Psy.D. Consortium, Palo Alto University, Palo Alto, CA, USA

<sup>4</sup> University of Massachusetts Medical School, Worcester, MA, USA

## **Edited by:**

Simon Baudrexel, Goethe University Frankfurt, Germany

## **Reviewed by:**

Marc Himmelbach, University Hospital Tuebingen, Germany Christian Lambert, St George's University of London, UK

## **\*Correspondence:**

Sara W. Lazar, Massachusetts General Hospital, Harvard Medical School, 149, 13th Street, Charlestown, Boston, MA 02129, USA e-mail: lazar@nmr.mgh.harvard.edu Individuals can improve their levels of psychological well-being (PWB) through utilization of psychological interventions, including the practice of mindfulness meditation, which is defined as the non-judgmental awareness of experiences in the present moment. We recently reported that an 8-week-mindfulness-based stress reduction (MBSR) course lead to increases in gray matter concentration in several brain areas, as detected with voxelbased morphometry of magnetization prepared rapid acquisition gradient echo MRI scans, including the pons/raphe/locus coeruleus area of the brainstem. Given the role of the pons and raphe in mood and arousal, we hypothesized that changes in this region might underlie changes in well-being. A subset of 14 healthy individuals from a previously published data set completed anatomical MRI and filled out the PWB scale before and after MBSR participation. PWB change was used as the predictive regressor for changes in gray matter density within those brain regions that had previously shown pre- to post-MBSR changes. Results showed that scores on five PWB subscales as well as the PWB total score increased significantly over the MBSR course. The change was positively correlated with gray matter concentration increases in two symmetrically bilateral clusters in the brainstem.Those clusters appeared to contain the area of the pontine tegmentum, locus coeruleus, nucleus raphe pontis, and the sensory trigeminal nucleus. No clusters were negatively correlated with the change in PWB. This preliminary study suggests a neural correlate of enhanced PWB. The identified brain areas include the sites of synthesis and release of the neurotransmitters, norepinephrine and serotonin, which are involved in the modulation of arousal and mood, and have been related to a variety of affective functions as well as associated clinical dysfunctions.

**Keywords: brain stem, mindfulness, well-being, stress, psychological, raphe nuclei**

## **INTRODUCTION**

Mindfulness meditation, a practice with origins in ancient Buddhist meditation traditions, has long been reported to produce positive effects on psychological well-being (PWB) that extend beyond the time the individual is actually meditating (Ekman and Davidson, 1994; Baer, 2003; Ekman et al., 2005). Taking advantage of these benefits, mindfulness practices have been increasingly incorporated into psychotherapeutic programs over the last three decades (Kabat-Zinn, 1990; Linehan, 1993; Roemer and Orsillo, 2002; Segal et al., 2002; Luoma et al., 2008). Mindfulness is defined as the purposeful and non-judgmental awareness of present-moment experience (Kabat-Zinn, 2003).

Research has shown that interventions incorporating mindfulness training positively affect symptoms of a variety of disorders including anxiety (Roemer et al., 2009; Hofmann et al., 2010), depression (Teasdale et al., 2000; Hofmann et al., 2010), and attention deficit hyperactivity disorder (Zylowska et al., 2008; van de

Weijer-Bergsma et al., 2012). Furthermore, preliminary evidence suggests that mindfulness-based interventions can positively influence sleep and dietary patterns in clinical populations (Baer et al., 2005; Winbush et al., 2007; Dalen et al., 2010).

Recently, neuroimaging studies have begun to explore changes in neural structure and function associated with meditation practice (Davidson, 2003; Lazar et al., 2005; Brefczynski-Lewis et al., 2007; Farb et al., 2007; Pagnoni and Cekic, 2007; Slagter et al., 2007; Hölzel et al., 2008; Lutz et al., 2008). A number of anatomical MRI studies have demonstrated that individuals who have regularly practiced meditation for several years exhibit a different gray matter morphometry in multiple brain regions when compared to demographically matched controls (Lazar et al., 2005; Pagnoni and Cekic, 2007; Hölzel et al., 2008; Luders et al., 2009; Vestergaard-Poulsen et al., 2009). Recently, we reported the first longitudinal study of gray matter changes following an 8-week-mindfulnessbased stress reduction (MBSR) course (Hölzel et al., 2011). One region with enhanced gray matter concentration following the

MBSR course was in the cerebellar vermis, reaching into a region of the brain stem that included the locus coeruleus, nucleus raphe pontis, pontine tegmentum, and the sensory trigeminal nucleus (Naidich et al., 2008). The locus coeruleus, a site of synthesis of norepinephrine, has been implicated in conditions such as depression and anxiety (Aston-Jones and Cohen, 2005). Furthermore, this region may play a role in modulating serotonin release (Plaznik and Kostowski, 1983; Grenhoff et al., 1993; Ressler and Nemeroff, 1999). The modulation of levels of serotonin, which is synthesized in the raphe nuclei, has been shown to be one of the most effective treatments for mood and anxiety disorders (Masand and Gupta, 2002). The pontine tegmentum, part of the cholinergic system, is implicated in regulating selective attention, wakefulness, learning, reward, and sleep (Kobayashi and Okada, 2007;Wang and Morales, 2009).

Given that these regions are well-known to modulate several systems, including the serotonin, dopamine, and norepinephrine systems, as well as play central roles in processes such as mood, arousal, sleep, and appetite (Ressler and Nemeroff, 1999; Aston-Jones and Cohen, 2005; Winn, 2006; Kobayashi and Okada, 2007; Wang and Morales, 2009; Bailer and Kaye, 2011), we reasoned that gray matter changes in these regions might contribute to enhanced well-being following mindfulness practice. A subset of individuals in our previous study had completed a questionnaire to assess PWB. Therefore, in order to test this hypothesis, we re-analyzed this subgroup of the previous data set and investigated correlations between changes in gray matter concentration and changes in self-report measures of PWB.

## **MATERIALS AND METHODS**

## **PARTICIPANTS**

The PWB scale was administered to a subsample of 14 participants from our previous study (Hölzel et al., 2011). As described previously (Hölzel et al., 2011), participants were recruited from MBSR courses held at the Center for Mindfulness at the University of Massachusetts Medical School. Individuals were included if they presented as physically and psychologically healthy, scored ≥1 SD above the population mean on the four-item Perceived Stress Scale (PSS; Cohen and Williamson, 1988), had no significant previous meditation experience, were between 25 and 50 years old, had no contra-indications for MRI scanning (i.e., metallic implants, claustrophobia, pregnancy), and made a verbal commitment to attend all eight classes and perform the prescribed daily meditation exercises.

The participants were healthy, right-handed individuals [five male and nine female; mean age: 37.9 years (SD: 4.3 years; age range: 29–44 years)]. Participants had an average of 17.5 years of education (SD: 1.9 years). Ethnicities were: 11 Caucasians, 1 South Asian, 1 African American, 1 multi-ethnic. Participants received a \$300 discount in the MBSR course fee (which costs between \$475 and \$630, depending on the household income) for their participation in the study. Additional analyses that included data from this sample have been reported elsewhere (Hölzel et al., 2010, 2011). The study protocol was approved by the IRBs of Massachusetts General Hospital and the University of Massachusetts Medical School and written informed consent was obtained from all participants.

## **MRI DATA COLLECTION AND ANALYSIS**

Participants were scanned at the Martinos Center for Biomedical Imaging in Charlestown, MA, USA, during the 2 weeks before (Pre) and after (Post) participation in MBSR. High-resolution MRI data were acquired with a Siemens Magnetom Avanto 1.5 T scanner with standard head coil. Data sets of the whole brain were collected using a T1 weighted, magnetization prepared rapid acquisition gradient echo (MPRAGE) sequence, consisting of 128 sagittal slices (voxel size: 1.0 mm × 1.0 mm × 1.3 mm, TI = 1000 ms; TE = 3.39 ms; TR = 2730 ms, flip angle 7°, matrix 256 mm × 256 mm). Image analysis was performed with VBM tools within the SPM5 neuroimaging statistical software (Wellcome Department of Cognitive Neurology, London, www.fil.ion. ucl.ac.uk/spm/software/spm5/) based in MATLAB 7.1, release 14 (Mathworks Inc., Natick, MA, USA), using default settings unless otherwise specified. Images were manually aligned to the anterior commissure and then segmented into gray and white matter in native space (i.e., before normalization, using the "Native Space" segmentation option implemented in SPM5). For each individual, the (unmodulated) gray matter segmentations of the Pre and Post images were spatially coregistered. Normalization parameters were calculated for the Pre scan and were applied to both time-points (trilinear interpolation, 2 mm × 2 mm × 2 mm), to make sure that regional differences between the images were not removed by scanspecific spatial normalization (Driemeyer et al., 2008; Ilg et al., 2008). Images were smoothed using an 8-mm full width at half maximum Isotropic Gaussian Kernel.

## **PSYCHOLOGICAL WELL-BEING**

Psychological well-being was assessed using the 54-item version of the PWB scale (PWB) by Ryff (1989). The PWB is based on a model comprising six factors of PWB (Ryff, 1989; Ryff and Keyes, 1995): self-acceptance (positive attitude toward oneself even while aware of one's own limitations), positive relations with others (developing and maintaining warm and trusting interpersonal relationships), environmental mastery (managing one's environment so as to meet personal needs and desires), autonomy (sense of self-determination and personal authority), purpose in life (sense of meaning in one's effort and challenges), and personal growth (view of self as growing and developing, openness to new experiences). These six factors integrate into a single second-order factor (Ryff and Keyes, 1995). The 54-item version of the PWB scale has been shown to have good psychometric properties (Sewell et al., 2004). Items were rated on a six-point continuum ranging from strongly disagree to strongly agree. The total score is derived by summing the scores on the six factors.

In a regression analysis using SPM5, Pre- to Post-intervention changes in the PWB total score were correlated with changes in gray matter concentration in regions that have previously been identified as showing an increase in gray matter concentration over the 8-week-MBSR course (Hölzel et al., 2011). In the previous study, we had identified these regions by performing a paired *t*-test within the group that had undergone the MBSR program, choosing a cluster-size threshold that was corrected for multiple comparisons across the entire brain (i.e., in order to exceed the threshold of *p* < 0.05, clusters had to exceed a size of 250 voxels) and based on statistical parametric maps with an initially, uncorrected, thresholded of *p* = 0.01. For the current study, we created a mask that contained the result of that previous study as our new region of interest. To be conservative, we included into this mask all the clusters that displayed a significant increase in gray matter concentration from Pre- to Post-intervention (cf. Table 2 and Figure 2 in Hölzel et al., 2011), namely the clusters within the brainstem/cerebellum, PCC, and left TPJ, i.e., four clusters with a total of 1537 voxels. To obtain images representing the change in gray matter concentration, the Pre-intervention scan was subtracted from the Post-intervention scan. Cluster level statistics for the current analysis are reported on an alpha level of <0.05, multiple comparisons corrected for the search region (height threshold: *p* = 0.01).

## **RESULTS**

## **IMPROVEMENTS IN PSYCHOLOGICAL WELL-BEING**

A paired-samples *t*-test revealed a significant increase in PWB from Pre- to Post-intervention (Pre mean: 224.64, SD: 28.62; Post mean: 252.75, SD: 26.89; *t* = −4.03; *p* = 0.001). Pre-Post changes for five of the six scales were also significant: self-acceptance (mean pre: 34.18, SD: 8.18, mean post: 40.93, SD: 6.67, *t* = −4.21, *p* = 0.001), environmental mastery (mean pre: 31.68, SD: 5.78, mean post: 36.68, SD: 6.47, *t* = −2.90, *p* = 0.012), autonomy (mean pre: 40.07, SD: 7.41, mean post: 44.14, SD: 6.15, *t* = −2.97, *p* = 0.011), purpose in life (mean pre: 38.14, SD: 8.02, mean post: 43.14, SD: 4.93, *t* = −2.66, *p* = 0.020), and personal growth (mean pre: 43.21, SD: 4.98, mean post: 47.21, SD: 4.93, *t* = −3.61, *p* = 0.003). The sixth scale, positive relations with others, revealed a trend toward significance (mean pre: 37.36, SD: 8.85, mean post: 40.64, SD: 7.04, *t* = −1.87, *p* = 0.084). When applying the very conservative Bonferroni multiple comparison correction for the six sub-tests, pre-post changes for the scales self-acceptance and personal growth remained significant, but all other subscales missed significance.

## **CORRELATION BETWEEN CHANGES IN PSYCHOLOGICAL WELL-BEING AND CHANGES IN GRAY MATTER CONCENTRATION**

To address the question of whether increase in gray matter concentration were related to improvements in well-being, the change in the total PWB score was regressed against changes in gray matter concentration within the regions identified in Hölzel et al. (2011). Within the chosen mask, two clusters in the brainstem were identified to be positively correlated with changes in PWB [**Figure 1**; right cluster: cluster-size *k*: 43 voxels; *p* = 0.024; MNI coordinates of peak voxel (*x*, *y*, *z*): 12, −36, −30; left cluster: cluster-size *k*: 37 voxels; *p* = 0.040; MNI coordinates of peak voxel (*x*, *y*, *z*): −14, −42, −32]. The more the participants' PWB improved over the 8-week-MBSR course, the more increase in gray matter concentration was observed in these regions. According to the atlas by Naidich et al. (2008), these clusters appear to contain the area of the pontine tegmentum, locus coeruleus, nucleus raphe pontis, and the sensory trigeminal nucleus bilaterally. For illustrative purposes, values were extracted and averaged across each cluster and plotted with the change in PWB total score (**Figure 2**). The Pearson coefficients were 0.72 (*p* = 0.004) for the correlation between PWB change and change in the left brainstem cluster, and 0.76 (*p* = 0.002) for the correlation between PWB change and

**FIGURE 1 | Correlation of improvements in psychological well-being and increase in gray matter concentration in the brainstem**. Axial slices from z = −38 to −24, with an image every 2 voxels. Significant clusters are overlaid over the group averaged normalized structural MPRAGE image. This analysis includes N = 14 participants.

change in the right brainstem cluster. Importantly, these numbers are reported only for comparative purposes, and should not be interpreted by themselves, since clusters were derived through searching for correlations with PWB scores in the first place. Using the Tukey-criterion of defining outliers as those values that are further than 1.5 times the interquartile range away from the upper or lower quartile (Tukey, 1977), we identified one single outlier, namely the individual with the highest change in PWB total score. When excluding this outlier from the analysis, the correlation coefficients dropped slightly, but remained significant for the right cluster (*r* = 0.716, *p* = 0.006) and almost significant for the left cluster (*r* = 0.553, *p* = 0.050). When additionally excluding the individual with the second highest change in PWB total score, the correlation with the left brainstem cluster was no longer significant (*r* = 0.165, *p* = 0.609), but the correlation with the cluster in the right side of the brainstem remained significant (*r* = 0.59, *p* = 0.043). No clusters or voxels were negatively correlated with the change in PWB. No clusters were negatively correlated with the change in PWB.

## **DISCUSSION**

We identified a positive correlation between improvement in PWB and increase in gray matter concentration within regions of the brainstem, suggesting that these morphological changes might be part of a mechanism underlying the changes in PWB. Regions within the brainstem were found to increase in gray matter concentration over the 8 weeks (Hölzel et al., 2011), and the increase within a sub-region of the original area of change was correlated with improvements in PWB. These regions appear to include the area of the locus coeruleus, pontine tegmentum, nucleus raphe pontis, and the sensory trigeminal nucleus (DeArmond et al., 1989; Naidich et al., 2008). Several previous cross-sectional studies have investigated the impact of meditation practice on brain morphology by comparing groups of experienced meditators to nonmeditators (Lazar et al., 2005; Pagnoni and Cekic, 2007; Hölzel et al., 2008; Luders et al., 2009; Vestergaard-Poulsen et al., 2009).

None of these previous studies assessed the participants' PWB, and all used a cross-sectional study design, which is usually not as sensitive as a longitudinal design and which suffers from well-known limitations (i.e., the possibility of pre-existing conditions or other life-style differences which may confound results).

The raphe nuclei are a major site of serotonergic neurons (Michelsen et al., 2007), which project widely throughout the brain. Serotonin is implicated in numerous functions including sleep (Monti and Monti, 2000; Monti, 2011), mood, appetite, and conditioned fear. Further, stress has been shown to downregulate serotonin receptors in the raphe nuclei (Fuchs and Flügge, 2003). Modulation of the serotonin system has been profoundly effective for the treatment of a wide range of mood and anxiety disorders (Masand and Gupta, 2002) and the serotonergic neurons of the dorsal raphe nuclei have been implicated in eating disorders (Bailer et al., 2007; Bailer and Kaye, 2011). Interestingly, mindfulness training has been shown to improve a number of conditions for which altered serotonin levels have been implicated, including anxiety and depression (Teasdale et al., 2000; Baer, 2003; Kuyken et al., 2008; Roemer et al., 2009), insomnia (Kuyken et al., 2008; Ong et al., 2009), eating disorders (Kristeller et al., 2006), as well as improvements in sleep patterns (Carlson et al., 2003; Carlson and Garland, 2005; Ong et al., 2009), and attention (Jha et al., 2007). Considering the importance of serotonin for a number of factors that contribute to PWB including mood and sleep patterns, our findings of increased gray matter concentration in the raphe following MBSR that is correlated with improved well-being is highly suggestive.

The locus coeruleus is the site of synthesis and release of the neurotransmitter norepinephrine, and is thought to optimize behavioral performance by modulating arousal, regulating the interplay between focused vs. flexible responding to environmental demands, or selective vs. scanning attention (Aston-Jones et al., 2000, 2007; Aston-Jones and Cohen, 2005). The neurons of the locus coeruleus are important in a variety of cognitive, affective, and other behavioral functions, as well as associated clinical dysfunctions (e.g., depression, anxiety, sleep, and circadian disorders; for discussion, see Aston-Jones et al., 2007). It is also one of the primary sites mediating the stress response as well as a site of action

of antidepressant drugs (Brady, 1994). This may be related to the influence of the norepinephrine system on perceptions of personal control and autonomy (Bandura et al., 1985; Ryff et al., 2006), which our results (autonomy and environmental mastery subscales) show are improved following MBSR participation. Compared to healthy controls, depressed individuals display reduced gray matter density in this region (Chan-Palay and Asan, 1989; Arango et al., 1996; Ressler and Nemeroff, 1999). Norepinephrine is thought to act as a modulatory agent, modulating serotonin and dopamine release through projections into the ventral tegmental area and the dorsal raphe nuclei (Plaznik and Kostowski, 1983; Grenhoff et al., 1993; Ressler and Nemeroff, 1999). Furthermore, there is evidence that bias toward negative memories and emotions in depression may be related to norepinephrine, and that potentiation of norepinephrine results in increased recognition of positive emotions and more positive emotional bias (Harmer et al., 2003). Several studies have found changes in serum concentration of serotonin and norepinephrine, particularly decrease in norepinephrine and increase in serotonin in meditators (Infante et al., 2001; Solberg et al., 2004; Curiati et al., 2005; Yu et al., 2011).

The pontine tegmentum and its nuclei, the pedunculopontine nucleus and the laterodorsal tegmental nucleus, are also part of the brain's cholinergic system, and have been indicated as working as a modulatory system influencing learning, reward, sleep/wakefulness, motor function, and attention (Kobayashi and Okada, 2007; Wang and Morales, 2009). The pedunculopontine nucleus and laterodorsal tegmental nucleus neurons send axons to dopamine-containing areas of the ventral tegmental area and the substantia nigra, as well as to the lateral hypothalamus, thalamus, and basal ganglia, wherein glutamine and acetylcholine may act to modulate reward and learning (Yeomans et al., 1993; Steiniger and Kretschmer, 2003). In addition, the pedunculopontine nucleus may play a role in associative learning and reward as a relay for contextual information to midbrain dopamine neurons (Pan and Hyland, 2005). Executive control processes contribute to PWB (e.g., autonomy, environmental mastery; Ryff et al., 2006) and the pontine tegmentum has also been implicated in REM sleep (Fuller et al., 2007), which also correlates positively with well-being (Ryff et al., 2006).

This study comes with several important limitations: first, the sample size is extremely small, and findings are therefore unreliable. Second, with a relatively low resolution of the acquired images and additional smoothing, spatial sensitivity is limited, which is especially relevant when looking at small nuclei. Third, it has been discussed in the literature that segmentation and normalization of the brainstem is particularly problematic (Beissner et al., 2011). The exact localization of the regions identified here therefore needs to be confirmed. Fourth, the brainstem search territory was defined by our previous between-group analysis, and as such might be considered "non-independent" from the current correlation analysis (e.g., Vul et al., 2009, but see also Lieberman et al., 2009; Poldrack and Mumford, 2009). As a consequence, we would like the results of this study to be understood as purely speculative, and hope that they might be used to generate hypotheses for future rigorous research.

An extensive body of research during the last decade has established that MBSR leads to improvements in psychological health and well-being (Grossman et al., 2004; Nyklícek and Kuijpers, 2008; Shapiro et al., 2008). Interestingly, the data presented here suggest well-being is associated with brainstem regions that are the primary production sites of several neurotransmitters and which modulate basic functions of survival (sleep, appetite) and mood/arousal. Given that chronic stress increases the likelihood of developing future psychopathology, knowledge of the neurobiological mechanisms of behavioral interventions used in reducing stress and promoting well-being will be of great clinical interest.

## **ACKNOWLEDGMENTS**

We wish to express our gratitude to all participants for their cooperation. We thank the Center for Mindfulness for conducting the Mindfulness-based stress reduction courses, and Nik Olendzki and Christina Congleton for support in data collection. This research was funded by the National Institutes of Health-NCCAM (R21- AT003425-01A2) and the British Broadcasting Company (BBC). Britta K. Hölzel was supported by a Marie Curie International Outgoing Fellowship within the Seventh European Community Framework Programme and the Kusala Foundation. Sara W. Lazar was supported by National Institutes of Health funding K01AT00694. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## **REFERENCES**


**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: 18 June 2013; accepted: 16 January 2014; published online: 18 February 2014. Citation: Singleton O, Hölzel BK, Vangel M, Brach N, Carmody J and Lazar SW (2014) Change in brainstem gray matter concentration following a mindfulness-based intervention is correlated with improvement in psychological well-being. Front. Hum. Neurosci. 8:33. doi: 10.3389/fnhum.2014.00033*

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

*Copyright © 2014 Singleton, Hölzel, Vangel, Brach, Carmody and Lazar. 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.*

## The ascending reticular activating system from pontine reticular formation to the thalamus in the human brain

#### **Sang SeokYeo<sup>1</sup> , Pyung Hun Chang<sup>2</sup> and Sung Ho Jang<sup>1</sup>\***

<sup>1</sup> Department of Physical Medicine and Rehabilitation, College of Medicine, Yeungnam University, Taegu, South Korea

<sup>2</sup> Department of Robotic Engineering, Graduate School, Daegu Gyeongbuk Institute of Science & Technology, Taegu, South Korea

## **Edited by:**

Simon Baudrexel, Goethe University Frankfurt, Germany

## **Reviewed by:**

Hubertus Axer, Jena University Hospital, Germany Sandro M. Krieg, Technical University Munich, Germany Johannes Christian Klein, Goethe University of Frankfurt, Germany

## **\*Correspondence:**

Sung Ho Jang, Department of Physical Medicine and Rehabilitation, College of Medicine, Yeungnam University, 317-1, Daemyungdong, Namku, Taegu 705-717, South Korea e-mail: strokerehab@hanmail.net

**Introduction:** Action of the ascending reticular activating system (ARAS) on the cerebral cortex is responsible for achievement of consciousness. In this study, we attempted to reconstruct the lower single component of the ARAS from the reticular formation (RF) to the thalamus in the normal human brain using diffusion tensor imaging (DTI).

**Methods:**Twenty six normal healthy subjects were recruited for this study. A 1.5-T scanner was used for scanning of diffusion tensor images, and the lower single component of the ARAS was reconstructed using FMRIB software. We utilized two ROIs for reconstruction of the lower single component of the ARAS: the seed ROI – the RF of the pons at the level of the trigeminal nerve entry zone, the target ROI – the intralaminar nuclei of the thalamus at the level of the commissural plane.

**Results:** The reconstructed ARAS originated from the pontine RF, ascended through the mesencephalic tegmentum just posterior to the red nucleus, and then terminated on the intralaminar nuclei of the thalamus. No significant differences in fractional anisotropy, mean diffusivity, and tract number were observed between hemispheres (p > 0.05).

**Conclusion:** We reconstructed the lower single component of the ARAS from the RF to the thalamus in the human brain using DTI. The results of this study might be of value for the diagnosis and prognosis of patients with impaired consciousness.

**Keywords: ascending reticular activating system, diffusion tensor imaging, reticular formation, consciousness**

## **INTRODUCTION**

Consciousness is an arousal and awareness of environment and self, which is achieved through action of the ascending reticular activating system (ARAS) on the brain stem and cerebral cortex (Daube, 1986; Paus, 2000; Zeman, 2001; Gosseries et al., 2011). The ARAS is composed of several neuronal circuits connecting the brainstem to the cortex. These neuronal connections originate mainly in the reticular formation (RF) of the brainstem and project through synaptic relays in the intralaminar nucleus of thalamus to the cerebral cortex (Daube, 1986; Paus, 2000; Zeman, 2001; Afifi and Bergman, 2005; Gosseries et al., 2011). In addition, several brainstem nuclei (locus coeruleus, dorsal raphe, median raphe, pedunculopontine nucleus, parabrachial nucleus), nonspecific thalamic nuclei, hypothalamus, and basal forebrain are also included in the ARAS system (Aston-Jones et al., 2001; Parvizi and Damasio, 2003; Fuller et al., 2011). Thorough evaluation of the ARAS is important for diagnosis and management of patients with impaired consciousness, such as patients who are in a vegetative state or those with minimal consciousness (Zeman, 2001; Gosseries et al., 2011).

Conventional brain MRI, functional neuroimaging techniques, electrophysiological methods, and MR spectroscopy have been used in studies of the ARAS in the human brain (Parvizi and Damasio, 2003; Schiff, 2006; Tshibanda et al., 2009, 2010;

Gawryluk et al., 2010). However, because the ARAS cannot be clearly discriminated from adjacent neural structures, accurate identification and estimation of the ARAS in the human brain can be problematic when using these methods. In contrast, diffusion tensor imaging (DTI) allows for evaluation of white matter because of its ability to image water diffusion characteristics (Mori et al., 1999). In normal white matter, water molecules have relative freedom of movement parallel to nerve fiber tracts. However, their movements are restricted across tracts, giving rise to diffusion anisotropy of white matter. Accordingly, diffusion anisotropy has been used for evaluation of the extent of fiber change in white matter (Chang et al., 2010; Puig et al., 2010). Several recent studies have attempted to demonstrate the usefulness of DTI for evaluation of lesions in patients with impaired consciousness and connectivity of specific ARAS nuclei in the human brain (Voss et al., 2006; Perlbarg et al., 2009; Tollard et al., 2009; Tshibanda et al., 2009; Fernandez-Espejo et al., 2010, 2011; Newcombe et al., 2010; Edlow et al., 2012). However, little is known about the whole reconstruction and estimation of the ARAS in the human brain (Edlow et al., 2012).

In the current study, using DTI, we attempted to reconstruct the lower single component of the ARAS from the pontine RF to the intralaminar nuclei of the thalamus in the normal human brain.

## **MATERIALS AND METHODS**

## **SUBJECTS**

Twenty six normal healthy subjects (14 males, 12 females; mean age, 31.85 ± 9.80 years; range, 20–50) with no history of neurologic disease were recruited for this study. All subjects participated in this study as volunteers and provided written consent before undergoing DTI scanning. The study was approved by the institutional review board of our hospital.

## **DIFFUSION TENSOR IMAGE**

DTI data were acquired using a 6-channel head coil on a 1.5-T Philips Gyroscan Intera (Philips, Best, Netherlands) with single-shot echo-planar imaging. For each of the 32 noncollinear diffusion sensitizing gradients, we acquired 67 contiguous slices parallel to the anterior commissure-posterior commissure line. Imaging parameters were as follows: acquisition matrix = 96 × 96, reconstructed to matrix = 128 × 128, field of view = 221 mm × 221 mm, TR = 10,726 ms, TE = 76 ms, parallel imaging reduction factor (SENSE factor) = 2, EPI factor = 49, and *b* = 1000 s/mm<sup>2</sup> , NEX = 1, and a slice thickness of 2.3 mm (acquired isotropic voxel size 2.3 mm × 2.3 mm × 2.3 mm).

## **PROBABILISTIC FIBER TRACKING**

Analysis of diffusion-weighted imaging data was performed using the Oxford Centre for Functional Magnetic Resonance Imaging of the Brain (FMRIB) Software Library (FSL; www.fmrib. ox.ac.uk/fsl). Affine multi-scale two-dimensional registration was used for correction of head motion effect and image distortion due to eddy current. Fiber tracking was performed using a probabilistic tractography method based on a multifiber model, and applied in the current study utilizing tractography routines implemented in FMRIB Diffusion (5000 streamline samples, 0.5 mm step lengths, curvature thresholds = 0.2) (Smith et al., 2004). Advantages of probabilistic tractography, which was used in this study, include greater robustness to noise, as well as the ability to detect pathways with sharper angles and to distinguish crossing fibers (Behrens et al., 2007; Winston et al., 2011).

The pathway of the ARAS was determined by selection of fibers passing through seed regions of interest (ROI) and target (termination) ROIs. A seed ROI was placed on the RF of the pons at the level of the trigeminal nerve entry zone (Daube, 1986; Afifi and Bergman, 2005). Analysis of the medial lemniscus and rubrospinal tract was performed in order to confirm the boundary of the RF on the pons (**Figure 1A**). For analysis of the medial lemniscus, seed ROIs were placed on the anteromedial medulla and the target ROI was placed on the somatosensory cortex (Hong et al., 2010). For analysis of the rubrospinal tract, seed ROIs were placed on the red nucleus and the target ROI was placed on the contralateral dorsolateral region of the medulla (Monakow's area) (Nathan and Smith, 1982; Kwon et al., 2011). The target ROI was given on the intralaminar nuclei of the thalamus at the level of the commissural plane (Morel, 2007). In defining the intralaminar nuclei of the thalamus, we referred to a brain atlas (Morel, 2007) (**Figure 1A**). Of 5000 samples generated from the seed voxel, results for contact were visualized at a threshold minimum of 1 streamlined through each voxel for analysis. Values of fractional anisotropy (FA), mean

diffusivity (MD), and tract number of the lower single component of ARAS were measured.

## **STATISTICAL ANALYSIS**

SPSS software (v.15.0; SPSS, Chicago, IL, USA) was used for data analysis. Paired *t*-test was used for determination of the difference in values of DTI parameters of the ARAS between the right and left hemispheres. Pearson correlation test was used for determination of correlation between DTI parameters of the ARAS and age. Results were considered significant when the *p* value was <0.05.

## **RESULTS**

We reconstructed the lower single component of the ARAS between the pontine RF and intralaminar nuclei of the thalamus. The reconstructed component of the ARAS originated from the pontine RF, ascended through the mesencephalic tegmentum just posterior to the red nucleus, and then terminated on the intralaminar nuclei of the thalamus at the level of the commissural plane in all subjects (**Figure 1B**).

Mean values for FA, MD, and fiber volume of the right ARAS were 0.43 ± 0.05, 0.96 ± 0.13, and 156.96 ± 112.06, respectively, and those of the left ARAS were 0.43 ± 0.04, 0.94 ± 0.13, and 159.73 ± 72.88, respectively. No significant differences in FA, MD, and tract number were observed between the right and left hemispheres (*p* > 0.05) (**Table 1**). In addition, results of the Pearson correlation test showed no significant age-related changes of FA, MD, and tract number (*p* > 0.05).

## **DISCUSSION**

In the current study, using DTI, we reconstructed one of the main pathways of the ARAS, the lower single component of the ARAS from the RF to the thalamus in normal subjects, although theARAS consists of additional brainstem nuclei, hypothalamus, basal forebrain, and thalamocortical projections to the cerebral cortex. We selected two ROIs for reconstruction of the lower single component of the ARAS: the seed ROI, which was the RF of the pons at the level of the trigeminal nerve entry zone (Daube, 1986; Afifi and Bergman, 2005), and the target ROI, which included the intralaminar nuclei of the thalamus (the central lateral nuclei, centromedian/parafascicular nuclei, and paracentral nuclei) at the level of the commissural plane (Morel, 2007). The rostral portion of the RF of the brainstem above the trigeminal nerve entry zone is known as the ARAS; in contrast, the caudal portion of the RF is involved in motor function and autonomic function related to cardiac and respiratory function (Daube, 1986). Therefore, we placed the seed ROI in the RF at the level of the trigeminal nerve entry zone. We placed the target ROI in the intralaminar nuclei, which are the main nuclei of the ARAS among the non-specific thalamic nuclei. Therefore, we believe that because we could not include the other thalamic nuclei concerned with the ARAS, the lower single component of the ARAS that was reconstructed in the current study is not the entire lower single component of the ARAS, but the main portion of the entire lower single component of the ARAS. Consequently, the lower single component of the ARAS originated from the pontine RF, ascended through the mesencephalic tegmentum posterior to the

**FIGURE 1 | (A)** Seed regions of interest (ROI) are given on the pontine reticular formation (Red color). The target ROI is given on the intralaminar nuclei of the thalamus at the level of the commissural plane. Boundary of the intralaminar nuclei of the thalamus was defined by reference to the text book

red nucleus, and then terminated on the intralaminar nuclei of the thalamus. In addition, values for the FA, MD, and tract numbers of the reconstructed lower single component of the ARAS did not differ significantly between the right and left hemispheres. The tract number is determined by the number of voxels contained of the brain atlas (Morel, 2007). ML, medial lemniscus; RST, rubrospinal tract; RF, reticular formation; AC, anterior commissure; PC, posterior commissure. **(B)** Pathways of the reconstructed ascending reticular activating system are shown at each level of the brain in a normal subject (26-year-old male).

within a neural tract (Kwak et al., 2010). The FA value indicates the degree of directionality and integrity of white matter microstructures such as axons, myelin, and microtubules, and the ADC value indicates the magnitude of water diffusion (Assaf and Pasternak, 2008).



Values represent mean (±SD), FA, Fractional Anisotropy; MD, Mean diffusivity; MD × 10<sup>−</sup><sup>3</sup> (mm<sup>2</sup> /s).

Several studies have demonstrated the clinical usefulness of DTI by estimating some areas of the lower single component of the ARAS from the RF to the thalamus in patients with impaired consciousness (Perlbarg et al., 2009; Tollard et al., 2009; Newcombe et al., 2010; Fernandez-Espejo et al., 2011). Tollard et al. (2009) reported on the usefulness of DTI, which was performed at the subacute stage for prediction of outcome in 45 patients with severe TBI (traumatic brain injury) (absence of response to simple orders). In their study, they measured the FA value at several supratentorial and infratentorial areas, including the anterior pons, posterior pons, and midbrain, and demonstrated that the decrease of infratentorial and supratentorial FA, except in the posterior pons, allowed for prediction of unfavorable outcomes 1 year from TBI. Perlbarg et al. (2009), who performed DTI scanning in 30 patients with an absence of response to simple orders following severe TBI, reported a definite decrease in FA

## **REFERENCES**


144–155. doi:10.1016/j.neuroimage. 2006.09.018


measured in the inferior longitudinal fasciculus, midbrain (cerebral peduncle and tegmentum), posterior limb of the internal capsule, and posterior corpus callosum in the unfavorable outcome group.Newcombe et al. (2010) used DTIfor characterization of the extent and location of white matter loss in patients who were in a vegetative state secondary to TBI (seven patients) and patients with ischemic-hypoxic injury (five patients). Abnormalities in the supratentorial areas were observed in both groups; in contrast, abnormalities of the brainstem were observed only in the TBI group. Fernandez-Espejo et al. (2011) used DTI in differentiation of the neuropathology of 25 vegetative and minimally conscious patients. They concluded that minimally conscious patients and those in a vegetative state differed in subcortical white matter and thalamic regions, but appeared not to differ in the brainstem. In a recent study using high angular resolution diffusion imaging, Edlow et al. (2012) reported on neuroanatomical connectivity of the ARAS in the human brain, both *in vivo* and *ex vivo*. They demonstrated that the connectivities of specific ARAS nuclei were implicated in arousal, and those of thalamic nuclei were implicated in modulation of arousal.

In conclusion, using DTI, we reconstructed the lower single component of the ARAS from the RF to the thalamus in the human brain. We believe that the methodology used and the results of this study may be helpful to researchers studying the ARAS in the human brain. However, one of the limitations of this study is that we were not able to fully elucidate the entire ARAS system because we did not include other thalamic and brainstem nuclei in our analysis which are also involved in the ARAS. Further studies on the clinical usefulness of our findings as well as studies on the projections of the ARAS from the thalamus to the cerebral cortex are needed.

## **ACKNOWLEDGMENTS**

This work was supported by the DGIST R&D Program of the Ministry of Education, Science, and Technology of Korea (13-BD-0401).

*Neurol.* 71, 531–546. doi:10.1097/ NEN.0b013e3182588293


of the ascending arousal system. *J. Comp. Neurol.* 519, 933–956. doi: 10.1002/cne.22559


imaging and functional implications. *J. Neurol. Neurosurg. Psychiatry* 81, 552–561. doi:10.1136/jnnp. 2009.196246


**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: 30 April 2013; accepted: 11 July 2013; published online: 25 July 2013. Citation: Yeo SS, Chang PH and Jang SH (2013) The ascending reticular activating system from pontine reticular formation to the thalamus in the human brain. Front. Hum. Neurosci. 7:416. doi: 10.3389/fnhum.2013.00416*

*Copyright © 2013 Yeo, Chang and Jang . This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## Imaging white matter in human brainstem

#### **Anastasia A. Ford1,2\*, Luis Colon-Perez <sup>3</sup> ,William T. Triplett <sup>4</sup> , Joseph M. Gullett 1,5,Thomas H. Mareci 1,4 and David B. FitzGerald1,6,7**

<sup>1</sup> Department of Veterans Affairs Rehabilitation Research and Development Brain Rehabilitation Research Center, Malcom Randall VA Medical Center, Gainesville, FL, USA


## **Edited by:**

Florian Beissner, Martinos Center for Biomedical Imaging, USA

## **Reviewed by:**

Christophe Habas, Université Pierre et Marie Curie, France Johannes Christian Klein, Goethe-University of Frankfurt, Germany

### **\*Correspondence:**

Anastasia A. Ford, Department of Veterans Affairs, Brain Rehabilitation Research Center, Malcom Randall VA Medical Center, 1601 SW Archer Road 151-A, Gainesville, FL 32608, USA e-mail: sokolova@ufl.edu

The human brainstem is critical for the control of many life-sustaining functions, such as consciousness, respiration, sleep, and transfer of sensory and motor information between the brain and the spinal cord. Most of our knowledge about structure and organization of white and gray matter within the brainstem is derived from ex vivo dissection and histology studies. However, these methods cannot be applied to study structural architecture in live human participants. Tractography from diffusion-weighted magnetic resonance imaging (MRI) may provide valuable insights about white matter organization within the brainstem in vivo. However, this method presents technical challenges in vivo due to susceptibility artifacts, functionally dense anatomy, as well as pulsatile and respiratory motion. To investigate the limits of MR tractography, we present results from high angular resolution diffusion imaging of an intact excised human brainstem performed at 11.1T using isotropic resolution of 0.333, 1, and 2 mm, with the latter reflecting resolution currently used clinically. At the highest resolution, the dense fiber architecture of the brainstem is evident, but the definition of structures degrades as resolution decreases. In particular, the inferred corticopontine/corticospinal tracts (CPT/CST), superior (SCP) and middle cerebellar peduncle (MCP), and medial lemniscus (ML) pathways are clearly discernable and follow known anatomical trajectories at the highest spatial resolution. At lower resolutions, the CST/CPT, SCP, and MCP pathways are artificially enlarged due to inclusion of collinear and crossing fibers not inherent to these three pathways. The inferred ML pathways appear smaller at lower resolutions, indicating insufficient spatial information to successfully resolve smaller fiber pathways. Our results suggest that white matter tractography maps derived from the excised brainstem can be used to guide the study of the brainstem architecture using diffusion MRI in vivo.

**Keywords: diffusion-weighted imaging, tractography, brainstem, white matter, high-resolution MRI**

## **INTRODUCTION**

The brainstem is a complex neural structure crucial for sustaining survival functions. These functions include respiration, cardiovascular regulation, sleep, consciousness, and transmission of sensory and motor information between the brain and the spinal cord (Nicholls and Paton, 2009). Knowledge of the structural organization in the brainstem originates from animal and human dissection and histology studies, and, in the more recent years, from human neuroimaging studies (Harting, 1977; Steriade et al., 1988; Stieltjes et al., 2001; Habas and Cabanis, 2007; Naidich et al., 2008; Kamali et al., 2010).

High-resolution magnetic resonance imaging (MRI) can be used to visualize the brainstem *in vivo*. Unfortunately, the brainstem is quite dense functionally with many nuclei and small cortical and spinal white matter projections. These structures can be difficult to visualize on standard clinical MRI scans as there is not sufficient tissue-specific contrast. Recently developed diffusionweighted imaging (DWI) is a non-invasive technique that allows visualization of gray and white matter architecture *in vivo*. In particular, DWI tractography can be used to infer white matter pathways based on the properties of water diffusion within underlying tissue. This method provides many valuable insights into white matter organization within the brain and also has been applied to visualize brainstem neural architecture (Nagae-Poetscher et al., 2004; Behrens et al., 2007; Habas and Cabanis, 2007; Catani and Thiebaut de Schotten, 2008; Rilling et al., 2008).

The anatomical complexity and location of the brainstem, however, presents a number of methodological challenges for *in vivo* tractography. In particular, pulsatile and respiratory motion resulting from the cardiac cycle and respiratory activity may produce significant distortions of the brainstem tissue during image acquisition. Cerebrospinal fluid (CSF) has been found to flow in both rostral and caudal directions during the cardiac cycle, with velocity through the aqueduct of Sylvius of roughly 14 mm/s in the cranial to caudal direction and −12 mm/s caudal to cranial (Sweetman and Linninger, 2011). This cardiac synchronized flow in the CSF can result in brainstem pulsations. Two different techniques have been suggested for reducing the distortions due to these pulsations. One, using a navigator-corrected approach results in reduced distortions, without complete elimination of cardiac pulsation (Jiang et al., 2002). The second is to use cardiac gating, which reduces distortion, although with longer acquisition times (Jiang et al., 2002). These techniques may not fully correct for pulsatile motion and may result in incomplete white matter modeling.

The present study investigates white matter organization within the brainstem using an excised tissue sample. This approach avoids motion-related artifacts and allows the limits of resolution to be studied by acquiring diffusion-weighted data at higher resolutions and magnetic field strengths than those currently possible *in vivo*. In addition, to examine how resolution affects the visualization of brainstem white matter pathways is altered as a function of data resolution, we acquired two additional diffusion-weighted MRI datasets of the same tissue using lower acquisition resolutions. These lower resolutions corresponded to the highest resolution currently available for whole brain *in vivo* acquisitions (1 mm isotropic), and to the most commonly used acquisition resolution in whole brain *in vivo* acquisitions (2 mm isotropic).

## **MATERIALS AND METHODS**

## **THE BRAINSTEM TISSUE**

The brain, including brainstem was obtained from the University of Florida Neuromedicine Human Brain Tissue Bank. The donor was a 78-year-old male who died of a cerebrovascular accident. The post-mortem interval prior to fixation was 18 h. The tissue had been fixed in 10% formalin after removal (**Figure 1**) and, prior to MRI measurement, was then washed in phosphate-buffered saline (PBS) to remove formalin. Written informed consent was obtained post-mortem from the patient's son in compliance with Institutional Review Board guidelines of the University of Florida and North Florida/South Georgia Malcom Randall Veteran's Affairs Medical Center.

## **IMAGE ACQUISITION**

Prior to MR data acquisition, the brainstem was washed for 5 days in PBS, with periodic solution exchanges, then placed in a 41 mm plastic cylinder filled with *Fluorinert* (3M, St. Paul, MN, USA). The brainstem was imaged with a pulsed field-gradient, spin echo sequence at 470 MHz in a 11.1 T 40 cm horizontal bore imaging spectrometer (Agilent Technologies, Santa Clara, CA, USA) using a custom-built volume transmit/receive coil. Three image datasets were acquired with isotropic resolution of (1) 0.333 mm, using a frequency by phase encode matrix of in 114 slices, (2) 1.0 mm, using a frequency by phase encode matrix in 40 slices, and (3) 2.0 mm, using a frequency by phase encode matrix of in 20 slices. The images were obtained with repetition times (TR's) of 5110.9, 4000, and 3000 ms for 0.333, 1, and 2 mm resolution, respectively, with an echo time (TE) of 31 ms, using two signal averages, gradient pulses of 4 ms width, and intervals of 18 ms

**FIGURE 1 | Ex vivo tissue sample**. **(A)** Coronal view displaying the anterior portion of the tissue sample. **(B)** Sagittal view displaying the left-hand side of the tissue sample.

between gradient pulses. The TE value was kept constant across the three acquisition resolutions to ensure that T2 effects did not differ across the acquisitions and did not result in signal-to-noise ratio (SNR) differences. The total field-of-view for each acquisition was 70 mm × 40 mm × 38 mm. Low diffusion-weighted images (100 s/mm<sup>2</sup> ) were acquired in six directions, and high-diffusionweighted images (1250 s/mm<sup>2</sup> ) were acquired in 64 directions distributed uniformly over a sphere following an electrostatic repulsion scheme (Jones et al., 1999).

The above image acquisition approach results in different SNRs for the three datasets. The lowest SNR was obtained with the 0.333 mm isotropic resolution data, where the SNR for low diffusion weighting (100 s/mm<sup>2</sup> ) is ∼32, and 25 for high-diffusionweighted data (1250 s/mm<sup>2</sup> ). As suggested in a recommendation that SNR should be above 20 to assure stable tensor estimations (Bastin et al., 1998), each dataset should contain sufficient SNR so that the contrast-to-noise ratio is suitable for stable tensor estimates and to produce consistent streamline pathways. In addition, the sampling scheme (Jones et al., 1999) used in this study and number of gradient directions also improve the tensor estimates (Batchelor et al., 2003).

Although the TE was held fixed, the TRs were shortened with decreasing resolution in order to reduce the total acquisition time for the three data sets (∼38 h). Therefore the longitudinal relaxation time (T1) weighting differed slightly between the three data sets. However, the minimum TR used in this study was still longer than brain T1 values previously reported (de Graaf et al., 2006; Oros-Peusquens et al., 2008; Wright et al., 2008). In particular, Wright and colleagues reported T1 relaxation times in human brain of 1.12 and 1.93 s for white and gray matter respectively at 7 T, and de Graaf and colleagues reported rat brain white and gray matter values of 1.75 and 2.11 s respectively at 11.7 T. Since T1 relaxation times in tissue appear to increase logarithmically with magnetic field strength, the measured human brain T1's at 7 T can

be projected to be appropriately between 1.3 and 2.2 s at 11.1 T. Also the excised fixed brain stem used in this study was immersion fixed with formalin after an unknown post-mortem interval (typically hours). Formaldehyde fixation is known to shorten T1 by as much as 21% (Shepherd et al., 2009b) and T1 is known to increase with post-mortem interval with a 20% increase after 24 h (Shepherd et al., 2009a). Therefore, the effect of fixation and postmortem on T1 may be assumed to balance. Using the estimated 11.1 T values of *in vivo* human brain T1's, the signal strength would be reduced by ∼12% between the longest TR (∼5 s for 0.333 mm resolution) and shortest TR (3 s for 2 mm resolution) TR. Thus, the effect of T1 relaxation would not appear to have significant effect on the estimated streamline tractography.

## **IMAGE PROCESSING**

MR data was processed with in-house software written in Interactive Data Language (IDL; Exelis Visual Information Systems, Boulder, CO, USA). The image intensity attenuation for each voxel was fitted as a linear decay to a rank-2 tensor dependent on the diffusion weightings, then the average diffusivity (AD) and fractional anisotropy (FA) values were calculated from the resulting tensor at each voxel (Basser and Jones, 2002). To infer fiber tract streamlines, the displacement probability of water self-diffusion in each voxel is estimated using the Mixture of Wishart method distributions (Jian and Vemuri, 2007) for an average diffusion displacement of 6µm then the maximum displacement probability within tissue in each voxel is identified. This method allows the tracking of crossing and branching fibers, and makes the current algorithm superior to traditional streamline tracking techniques based on the rank-2 tensor models (Basser et al., 1994; Jian et al., 2007). Tractography is performed by seeding each voxel in the brain with a sub-voxel grid with evenly spaced seed points. From each seed point, one streamline is launched bi-directionally for each estimated displacement probability maximum contained in that voxel using the FACT algorithm (Mori et al., 1999). Each streamline front is propagated by stepping 0.5 voxel width in the direction of the maximum that is most inline with the streamline's present direction of travel. In order to prevent streamlines from looping back, angular deviation of the track is limited to 65°. If the estimated track exceeds this threshold, the streamline is stopped.

In order to determine a seeding density that would produce consistent results across the three acquisition resolutions, each dataset was seeded with 8, 27, and 64 seed points per voxel, then the percent change in edge weight (see below) was compared across the three seeding densities for each acquisition (**Table 1**). Edge weight is a dimensional quantity that represents tract connectivity strength, independent of acquisition resolution and seeding density. Our results for the 0.333 mm isotropic dataset showed that the percent change in edge weight values for all pathways differed by less than 7% as we increased the seeding density from 8 to 27 seeds per voxel. As seeding density was increased from 27 to 64 seeds per voxel, the percent change was less than 2% for all pathways (**Table 1**). For the 1 mm isotropic dataset, the percent change in edge weights was less than 3% for the corticopontine/corticospinal, and for the middle cerebellar peduncle (MCP) pathways as seeding density was increased from 8 to 27 and also from 27 to 64 seeds per voxel. The percent change for the superior cerebellar peduncle pathways was 32.49% as seeding density was increased from 8 to 27 seeds, and 19.88% as seeding density was increased from 27 to 64 seeds per voxel (**Table 1**). Due to this large difference in edge weight values, we seeded the 1 mm isotropic dataset using 125 seeds per voxel and computed the edge weight value for the superior cerebellar peduncle pathways. Edge weight value for this seeding density differed by less than 1% from the edge weight value computed by the 64 seeds per voxel seeding density. Seeding density results for the 2 mm isotropic dataset showed that for all pathways the edge weight values differed by less than 6% as we changed the seeding density from 8 to 27 seeds and from 27 to 64 seeds per voxel (**Table 1**). Overall, these results show that seeding density of 64 seeds per voxel provide a stable result, so the tractography results presented below are pathways traced using a track seeding density of 64 seeds per voxel.

These image processing methods were applied to each of the acquired scans (333.3µm, 1, and 2 mm isotropic). To determine how resolution effects the inference of fiber pathways, we compared trajectories and diffusion characteristics for pathways of interest across the three acquisition resolutions from the best available *ex vivo* resolution (0.333 mm isotropic) to the best possible resolution *in vivo* (1 mm isotropic), and to a commonly used *in vivo* resolution (2 mm isotropic) (see **Figures 3**–**5** below).

## **REGIONS OF INTEREST**

In order to define representative white matter pathways within the excised tissue, we used Duvernoy's atlas (Naidich et al., 2008) to create six regions of interest (ROI's) to define the boundaries of four white matter fiber pathways; (1) corticopontine/corticospinal (CST/CPT) fibers, (2) fibers passing through the superior cerebellar peduncle (SCP), (3) fibers passing through the MCP, and (4) the medial lemniscus (ML) pathways. The regions were drawn on the FA map, which provides good white/gray matter contrast.

To delineate the CS/CP pathway, we created two ROI's within the pons. The first delineated the cerebral peduncle and is drawn within the superior-most axial slice where this structure can be clearly identified. We delineated the cerebral peduncle within this slice and six additional slices inferior to it, creating a sevenslice ROI (**Figure 2A**). The inferior border of the second ROI was located within the inferior-most axial slice of the pons and followed six slices superior to it, creating a seven-slice ROI (**Figure 2B**).

To create the ROI's to identify the SCP fibers (**Figure 2C**), we first located the superior-most axial slice where the SCP could be identified definitively. Next, we identified the most inferior portion of the peduncle and recorded the corresponding axial slice number. We then computed the axial slice number corresponding to a mid-way point between the superior and inferior borders of the SCP. Lastly, we outlined the peduncle within the mid-point axial slice and the three slices inferior and superior to it, creating a seven-slice ROI. In addition, we created a second, or waypoint, ROI to infer cerebellar efferent traveling through the SCP toward the thalamus and the cortex. The waypoint ROI was located contralateral to the SCP ROI, and was drawn on the same axial slices as the


## **Table 1 | Percent change in edge weights as a function of seed density.**

Values listed in the table represent percent change in edge weight values for corticopontine/corticospinal and middle cerebellar peduncle pathways. Superior cerebellar peduncle pathways edge weight percent changes are listed in parenthesis.

**FIGURE 2 | Regions of interest used to delineate white matter pathways within the brainstem**. **(A)** cerebral peduncle mask, **(B)** inferior pons mask, **(C)** superior cerebellar peduncle (SCP) mask, **(D)** middle cerebellar peduncle mask (MCP), **(E)** superior medial lemniscus (ML) mask, and **(F)** inferior ML mask.

cerebral peduncle ROI to outline the remainder of the brainstem tissue outside of the cerebral peduncle.

The MCP ROI was drawn following the steps described below. First, we identified the most superior axial slice corresponding to the superior border of the MCP. Next, we identified the inferior extent of the MCP and recorded the corresponding axial slice number. We then computed the axial slice number corresponding to the mid-point between the superior and inferior extents of the MPC. We outlined the MCP within that slice and the three slices superior and inferior to it creating a seven-slice ROI (**Figure 2D**). We used the cerebral peduncle ROI (**Figure 2A**) as a waypoint to infer pathways traveling to the cerebellum from the corticopontine nuclei and the cortex. The cerebral peduncle ROI was drawn on the contralateral side of the brainstem tissue.

To delineate the ML pathways we created two ROIs. The first ROI was drawn on an axial slice corresponding to the superior extent of the pathway determined by examining the sagittal view of the FA image, as well as on nine consecutive inferior axial slices (**Figure 2E**). The final ROI included voxels comprising the ML pathway within each of the 10 axial slices. The second ROI was drawn on an axial slice corresponding to the inferior extent of the ML pathway identified by the sagittal view of the FA map, as well as on nine superior axial slices (**Figure 2F**). The second ROI consisted of 10 axial slices and included voxels comprising ML pathways within each of these slices. To ease the tracking process we chose to include more axial slices in the ML ROIs as compared to the remaining ROIs because this pathway is smaller than the other three pathways we examined.

All ROI's were drawn on an FA image derived from the highest resolution diffusion scan (0.333 mm) and registered to 1, 2 mm scans using linear registration with nearest-neighbor interpolation using FSL FLIRT (Jenkinson and Smith, 2001). **Table 2** lists surface area values for each of these ROI in all five datasets. Surface area measurements are used in edge weight value calculations (see Eq. 1 below).

## **TRACTOGRAPHY ANALYSIS**

To infer the corticopontine/corticospinal tracts (CPT/CST) pathways, we filtered whole brainstem streamline tractography results to include only fibers that intersect the cerebral peduncle ROI then filtered these fiber pathways further with the inferior pons ROI so only pathways connecting both ROI's were kept for further analysis. In order to delineate the SCP pathways, we filtered whole brainstem streamline tractography results with the SCP ROI to infer pathways passing through this region. Next, we filtered the resulting pathways with the waypoint ROI to infer cerebellar efferent fibers traveling toward the thalamus and the cortex. Similarly, to infer MCP pathways we first filtered whole brainstem streamline tractography results with the MCP ROI. Then we applied the cerebral peduncle ROI as a waypoint to visualize cerebellar afferents traveling from the cortex and pontine nuclei toward the cerebellum. The ML pathways were inferred by first intersecting the whole brainstem tractography results with the superior ML ROI, and then further filtering these results by applying the inferior ML ROI.


**Table 2 | Surface areas of regions of interest used to infer white matter pathways.**

In order to visualize how the appearance of the streamline pathways is changed as a function of acquisition resolution, we created overlap images of the four pathways inferred at three acquisition resolutions (**Figure 7**). Specifically, after delineating each of the four streamline pathways in each dataset, we registered the resulting streamline pathways from the 1 to 2 mm datasets to the 0.333 mm dataset using linear nearest-neighbor transformation using FSL FLIRT (Jenkinson and Smith, 2001). The streamline pathways were binarized; if a voxel was occupied by a streamline pathway, it was assigned a value of 1. If a voxel was not occupied by a streamline pathway, it was assigned a value of 0. The binarized streamline pathways (0.333 mm pathways, 1 mm pathways registered to 0.333 mm, and 2 mm pathways registered to 0.333 mm) were combined into a single pathway by summing the voxel values. If a voxel belonged to a given streamline pathway in all three datasets this voxel would be assigned a value of 3, if only two datasets had this voxel in common it would be assigned a value of 2. Similarly, if a voxel was unique to one of the datasets, it would be assigned a value of 1. This approach allowed us to examine the spatial variability of the inferred streamline pathways as a function of acquisition resolution. In addition, we calculated percent overlap between the 0.333 mm streamlines and registered 1 and 2 mm streamlines (**Table 4**).

## **QUANTITATIVE TRACTOGRAPHY MEASURES**

After calculating the streamline tracts along these pathways, we quantified the strength of fiber connectivity for the resulting tracts using measures of tract volume and edge weight. Tract volume represents the number of voxels occupied by each of the streamline tracts multiplied by the volume of each voxel (measured in mm<sup>3</sup> ) and the edge weight represents strength of connectivity of each pathway and is derived as an application of graph theory (Hagmann et al., 2008; Colon-Perez et al., 2012).

Edge weight, *w*(*eij*), is a scalar, dimensionless quantity independent of the acquisition resolution and seeding density. This measure is computed for inferred fiber tracts between two ROI (Eq. 1 below). From graph theory perspective, each ROI defines a network node (e.g., *n<sup>i</sup>* or *nj*) and each pathway connecting a pair of nodes, *n<sup>i</sup>* and *n<sup>j</sup>* , represents an network edge (e.g.,*eij*). In Eq. 1, the edge weight is dependent on the inverse sum of the streamline path length, *l*(*fp,m*), between the nodes. The inner summation involves all seed points (*P*voxel) in each voxel (*p* = *1* to *P*voxel) that belong the edge. The outer summation involves all voxels (*M*) in the edge (*m* = *1* to *M*), connecting the nodes *n<sup>i</sup>* and *n<sup>j</sup>* , that originate only from voxels in the streamline path connecting the nodes (located at *xp,m*, *yp,m*, and *zp,m*).

$$\begin{aligned} \, \, \_W \left( e\_{\vec{\imath} \vec{\jmath}} \right) &= \left( \frac{V\_{\text{voxel}}}{P\_{\text{voxel}}} \right) \left( \frac{2}{A\_i + A\_j} \right) \\ &\times \sum\_{m=1}^{M} \sum\_{p=1}^{P} \frac{\delta \left( x - x\_{p,m} \right) \delta \left( y - y\_{p,m} \right) \delta \left( z - z\_{p,m} \right)}{l \left( f\_{p,m} \right)} \end{aligned} (1)$$

This double sum over the inverse path length [1/*l*(*fp,m*)] effectively provides a count of the number of unique streamlines in the pathway, which is not dependent on the length of the streamlines.

Terms to the left of the double sum represent scaling factors that result in a dimensionless quantity. Specifically, *A<sup>i</sup>* and *A<sup>j</sup>* are surface areas of the two ROI defining the nodes. Scaling by the average of these surface areas allows direct comparison of edge weights for pathways with various size ROI. Multiplying by the ratio of voxel volume, Vvoxel, and the number of seed points, Pvoxel, represents a normalization by the seed point density used during the streamline tractography analysis to ensure that the edge weight is independent of the acquisition resolution and seed density. This approach allows us to compare edge weights for data acquired using different acquisition resolutions and analyzed using different number of seeds per voxel.

## **RESULTS**

Using our tractography approach we delineated the corticopontine/corticospinal pathways, MCP pathways, superior cerebellar peduncle pathways, and ML pathways. **Table 3** lists track volumes and edge weight values for each of these pathways inferred from the three acquisition resolution datasets (0.333, 1, and 2 mm isotropic). Overall, we note that tract volumes and edge weight values of the fiber bundles increase in magnitude as acquisition resolution is decreased from 0.333 to 2 mm (with exception of the superior cerebellar pathways where tract volume and edge weight are smaller for 1 mm dataset as compared with 0.333 and 2 mm results and the ML pathways where tract volume and edge weight for the 2 mm dataset are smaller than for the 0.333 and 1 mm datasets.). **Table 2** indicates that surface areas of ROI used to infer the pathways did not vary significantly as a function of data resolution (with the exception of the ML ROI). This observation implies that larger track volumes and edge weight values observed in lower acquisition resolution datasets did not result from volume changes in the ROI. Rather, as resolution is decreased from 0.333, to 1 and 2 mm resolution, volume averaging effects artificially inflate track volumes and edge weight values as extraneous pathways are included along with the pathways of interest. Larger acquisition voxels provide less information about branching and


**Table 3 | Tract volumes and edge weight values of the brainstem white matter pathways inferred from three acquisition resolution datasets.**

<sup>a</sup>Corticopontine/corticospinal tract; <sup>b</sup>superior cerebellar peduncle tracts; <sup>c</sup>middle cerebellar peduncle tracts; <sup>d</sup>medial lemniscus tracts.

crossing fibers as well as about tissue-type boundaries. This is particularly important for brainstem anatomy as gray matter nuclei and their white matter projections are dense and relatively small in size.

We note that the surfaces areas of the ROI used to delineate the ML pathways are smaller in 1 and 2 mm datasets as compared to the 0.333 mm dataset (**Table 2**). This finding seems to suggest that due to the effects of volume averaging, these lower resolution scans do not contain enough spatial information to successfully identify voxels within small white matter pathways such as the ML. This effect greatly limits the size of the pathways that can be successfully resolved at these lower resolutions, as the ROIs necessary to delineate the pathways are difficult to identify accurately, especially in the 2 mm dataset.

## **CORTICOPONTINE/CORTICOSPINAL PATHWAYS**

**Figure 3** depicts the corticopontine/corticospinal pathways inferred from the three datasets with different acquisition resolutions and from the interpolated datasets. Corticopontine/corticospinal pathways inferred from the 0.333 mm isotropic dataset connect the cerebral peduncle with the inferior pons. The streamline pathway occupies most of the anterior-posterior extent of the pons. In addition, a small secondary streamline pathway can be seen on the sagittal view of the pathways (**Figure 3A**, top row).

Corticopontine/corticospinal pathways inferred in the 1 mm dataset also connect the cerebral peduncle and the inferior pons ROI's. In addition, and unlike the 0.333 mm dataset, an additional, secondary streamline courses within the posterior portion of the pons (**Figure 3B**, top row). This secondary streamline likely includes a small portion of the ML and/or pathways within the anterior lateral system (ALS), due to its close spatial proximity to the inferior pons ROI (Haines, 2004). The inferred streamline pathways also appear to include a small portion of the fibers that course medially and cross over to the contralateral side of the brainstem (**Figure 3B**, bottom row). These contralateral pathways are likely to be the cerebellar afferents that pass through the cerebral peduncle and cross over to the contralateral MCP. Increase prevalence of extraneous streamline pathways, observed in the 1 mm dataset, is likely the result of increased volume averaging present in this lower resolution scan.

**Figure 3C** depicts corticopontine/corticospinal streamline pathways inferred in the 2 mm dataset. Overall the majority of the pathways look consistent with those inferred from the 0.333 and 1 mm datasets. Interestingly, the extraneous pathways, likely to be comprised by the ML/ALS fibers observed in the 1 mm dataset,

**FIGURE 3 | Corticopontine/corticospinal pathways: (A) 0.333 mm isotropic acquisition resolution, (B) 1 mm resolution, (C) 2 mm resolution**. Top row: sagittal view of the pathways with a mid-sagittal slice of the FA maps serving as a background. Middle row: coronal view of the pathways with an axial slice of the FA maps serving as a background. Bottom row: axial view of the pathways with an axial slice of the FA maps serving as a background. Regions of interest used to delineate the pathways (cerebral peduncle and inferior pons) depicted in pink. The regions of interest are not included in the axial view images for ease of visualization of the pathways. Three-dimensional cubes to the right of the images represent spatial rotations of the pathways, where faces of the cubes represent orientations of the x, y, and z axes in space: R (right, x-axis), A (anterior, y-axis), S (superior, z-axis). Color gradient within the pathways represent local fiber orientation: red – medial/lateral; blue – superior/inferior; green – anterior/posterior.

are not present in the 2 mm dataset. The 2 mm resolution dataset does contain extraneous streamline pathways, but they are located within the posterior medial portion of the pons and are likely to be the medial longitudinal fasciculus (MLF) tracts (**Figure 3C** bottom row).

## **SUPERIOR CEREBELLAR PEDUNCLE**

**Table 3** lists tract volumes and edge weight values for the superior cerebellar peduncle streamline pathways inferred from the three datasets with different acquisition resolutions and from the interpolated datasets. Both the track volume and edge weight value for the 0.333 and the 2 mm datasets are larger in magnitude than those for the 1 mm dataset. **Figure 4** depicts the superior cerebellar streamline pathways inferred in 0.333, 1, and 2 mm datasets.

Superior cerebellar peduncle streamline pathways, that carry cerebellar efferents toward the thalamus and cortex and traveling through the SCP, cross over to the contralateral side of the brainstem within the cerebellar peduncle decussation (indicated by red color within the streamline pathways in **Figure 4**), and ascend further within the superior extent of our brainstem sample. We observe that these streamline pathways appear to include some extraneous tracts in the 0.333 mm dataset, where these additional pathways descend down into the contralateral side of the brainstem (**Figure 4A**). These extraneous pathways are absent in the 1 mm dataset, where a vast majority of the pathways follow known anatomical trajectory of the cerebellar efferents (Haines, 2004). The majority of the SCP streamline pathways within the 2 mm dataset originate in the medial extent of the SCP, allowing us to visualize only a fraction of the SCP pathways (**Figure 4C**).

## **MIDDLE CEREBELLAR PEDUNCLE**

**Table 3** indicates that the MCP streamline pathways, that carry cerebellar afferents from the thalamus and cortex, increase in volume and have larger edge weight values as resolution decreases. Specifically, the 0.333 mm isotropic dataset has the smallest tract volume, followed by the 2 and 1 mm datasets.

**Figure 5** depicts the MCP streamline pathways in the three datasets. Overall the streamline pathways project from the cerebral peduncle to the ipsilateral and contralateral pontine nuclei anteriorly, and to the ipsilateral and contralateral retrolenticular tegmental nuclei posteriorly before entering the MCP. The majority of the streamline pathways in the 0.333 and 1 mm datasets are located predominantly within the anterior extent of the pons consistent with the notion that the cerebellar afferents are more prominent within this region (**Figures 5A,B**) (Naidich et al., 2008). In the 2 mm resolution dataset, projections to the retrolenticular tegmental nuclei are more prominent that those in 0.333 and 1 mm datasets, as streamline pathways within the posterior extent of the pons are more robust (**Figure 5C**).

## **MEDIAL LEMNISCUS PATHWAYS**

Medial lemniscus pathways carry sensory information from the gracile and cuneate nuclei to the thalamus. **Figure 6** represents these pathways in the excised brainstem sample. We note that in the three acquired datasets, the streamline pathways travel within the posterior portion of the pons close to the midline of the brainstem (**Figures 6A–C**). In the 0.333 mm dataset, the streamlines are tightly packed into a single pathway (**Figure 6A**), while a secondary lateral streamline pathway is present in the 1 mm data (**Figure 6B**). This secondary streamline pathway is likely to be comprised of the ALS fibers located in the posterior lateral extent of the brainstem. The streamline pathways in the 2 mm dataset appear to be significantly smaller than those in the 0.333 and 1 mm datasets

**FIGURE 4 | Superior cerebellar peduncle pathways: (A) 0.333 mm isotropic acquisition resolution dataset, (B) 1 mm isotropic resolution dataset, (C) 2 mm isotropic resolution dataset**. An axial slice of the FA maps serves as a background image. Superior cerebellar peduncle (SCP) and a waypoint mask are depicted in yellow. Color gradient and three-dimensional cube colors are the same as in **Figure 3**.

(**Figure 6C**). This observation is supported by lower tract volume and edge weight for the 2 mm ML pathways (**Table 3**). The smaller number of streamline pathways observed in the 2 mm dataset is likely a result of the size of the ROIs used to delineate this pathway. Specifically, as we noted above, the superior and inferior ML ROIs had significantly smaller surface areas as compared with the 0.333 and 1 mm datasets. Insufficient spatial information resulting from volume averaging present in the 2 mm dataset results in loss of gray matter/white matter contrast in the FA image, making it difficult to identify voxels that belong to the ML pathway. Specifically, the gray matter surrounding the ML pathway contributes more as resolution at lower resolution and this results in lower FA values and a smaller number of voxels that can be definitively identified as the ML pathway. Spatial resolution is therefore a limiting factor for the size of the pathways that can be successfully resolved using tractography, especially at lower acquisition resolutions.

## **SPATIAL VARIABILITY OF THE INFERRED STREAMLINE PATHWAYS AS FUNCTION OF ACQUISITION RESOLUTION**

In order to examine spatial variability of the inferred streamline pathways due to acquisition resolution we registered streamlines inferred in the 1 and 2 mm datasets to the 0.333 mm dataset and examined the percent overlap for each of the four pathways. **Table 4** displays the resulting percent overlap values. For the CPT/CST streamline pathways, the percent overlap was 63.12% for the 1 mm dataset and 68.63% for the 2 mm dataset. **Figure 7A** depicts the overlapping streamline pathways for this pathway. We note that the core of the pathway is common to the three acquisition resolution datasets as indicated by the highest amount of spatial overlap (yellow color in **Figure 7** represents voxel values of 3 which implies that the streamline pathways from all three acquisition resolution datasets had these voxels in common). Voxels surrounding the core of the pathway were common to two acquisition resolution datasets (red color in **Figure 7** indicates that two datasets had these voxels in common). Most of the voxels

comprising the extraneous pathways, such as ML and contralateral projections, were common to a single dataset (as indicated by blue color in **Figure 7A**). These results show that although the core of the pathway is inferred consistently as resolution is decreased from 0.333 to 1 and 2 mm additional extraneous pathways are also picked up due to volume averaging.

**Figure 7B** represents the overlap map for the SCP streamline pathways in the three acquisition datasets. We note that there are no overlapping voxels common to all three datasets (as indicated by the lack of yellow colored voxels). Overall, there was a small amount of overlap between streamlines from the 0.333 mm dataset and 1 and 2 mm data (**Table 4**). We observed a 28.29% overlap between 0.333 and 1 mm datasets, and 10.21% overlap between 0.333 and 2 mm datasets. This result is consistent with the earlier finding that the SCP streamline pathways have distinct trajectories in the three acquisition datasets (**Figure 4**).

**Figure 7C** represents the overlap map for the MCP pathways. We note that similarly to the CPT/CST overlap map, the core of the MCP pathway is common to the three acquisition dataset streamlines, as indicated by the presence of yellow colored voxels (**Figure 7C**). Specifically, the core of the pathway clusters around the MCP (**Figure 7C** axial view), and the contralateral cerebral peduncle (**Figure 7C** sagittal view). In addition, a smaller portion of overlapping voxels can be seen within the contralateral pons (**Figure 7C** axial view). The small number of overlapping voxels within the contralateral pons is not surprising if we examine the MCP pathways in the three acquisition resolution datasets (**Figure 5**). This is to be expected because the 0.333 and 1 mm streamline pathways occupy predominantly the anterior portion of the pons, while the 2 mm streamline pathways travel within the posterior portion. These differences in trajectories result in a small amount of overlap within the contralateral pons. Therefore, the majority of overlapping voxels are located within close vicinity of the MCP ROI and within contralateral cerebral peduncle. Specifically, we found a 70.10% overlap between the 0.333 and 1 mm streamlines, and a 67.95% overlap between the 0.333 and 2 mm streamlines (**Table 4**).

**Figure 7D** represents the overlap map of the ML streamline pathways. We note that the majority of the voxels common to the streamline pathways in the three acquisition resolutions are clustered within the core of the pathway. However, since the 2 mm dataset inferred streamline pathway is much smaller than those in the 0.333 and 1 mm datasets, the overlapping core (depicted as yellow colored voxels in **Figure 7D**) accounts for less than a half of all voxels comprising the overlap map. In addition, the percent overlap between the 0.333 and 1 mm dataset streamlines was 47.31%, and only 30.44% between the 0.333 and 2 mm datasets (**Table 4**). The smaller amount of overlap between the 0.333 and

**Table 4 | Percent overlap between streamline pathways traced in the 0.333 mm dataset and those traced in the 1 and 2 mm datasets.**


2 mm streamline pathways results from volume averaging effects present at this lower acquisition resolution.

## **DISCUSSION**

Tractography studies have provided many valuable insights into white matter organization within the brain (Basser et al., 2000; Johansen-Berg et al., 2005; Catani and Thiebaut de Schotten, 2008; Frey et al., 2008; Rilling et al., 2008; Edlow et al., 2012). Brainstem white matter visualization, however, has proved to be challenging due to the relatively small size of the pathways, high density of their distribution throughout the structure, and image distortions associated with *in vivo* acquisitions (Stieltjes et al., 2001; Habas and Cabanis, 2007; Kamali et al., 2010). The present study aimed to address some of these technical challenges by acquiring high-resolution *ex vivo* data using high field strength MR scanner. In addition, we investigated how the visualization of white mater tracts at the highest currently available acquisition resolution (0.333 mm isotropic) compares with the best *in vivo* acquisition resolution (1 mm isotropic) and the most commonly used *in vivo* acquisition resolution (2 mm isotropic). Our results indicate that although most of white matter streamline pathways examined in this study follow known anatomical trajectories in the three acquisition resolution datasets, volume averaging effects result in inclusion of extraneous pathways in lower resolution datasets. Inclusion of these additional pathways results in enlargement of the inferred fiber bundles reflected by larger tract volumes and edge weight values observed in the present study. Visualizations of white matter pathways in 1 and 2 mm datasets contain less spatial detail and appear more smooth and contiguous as compared with the 0.333 mm pathways (**Figures 3**–**5**). Volume averaging effects observed in the present study in 1 and 2 mm datasets were previously reported by *in vivo* tractography studies using similar or lower acquisition resolutions (Stieltjes et al., 2001; Habas and Cabanis, 2007).

Another important result observed in the present study was the finding that the largest changes in overall appearance of the pathways as well as increases in track volumes and edge weight values occur as resolution was decreased from 0.333 to 1 mm (**Table 3**, with exception of the SCP and the ML pathways). Specifically, we noted a 43% increase in tract volume and a 140% increase in edge weight for the corticopontine/corticospinal pathways as resolution was decreased from 0.333 to 1 mm. Similarly, a 94% increase in track volume and a 2500% increase in edge weight was observed for the MCP pathways. In contrast, when resolution was decreased from 1 to 2 mm, we observed a 3% decrease in track volume and a 56% increase in edge weight for the corticopontine/corticospinal pathways. For the MCP pathways, we observed a less than 2% reduction in track volume and a 2% increase in edge weight values. Large changes in track volumes and edge weight values as resolution is decreased from 0.333 to 1 mm is indicative of substantial loss of spatial information and volume averaging effects. This finding implies that the relative tissue-type contrast loss due to a decrease in resolution from 0.333 to 1 mm results in significantly less accurate modeling of white matter bundles using tractography. As resolution is decreased further from 1 to 2 mm the inferred white matter pathways do not appear to undergo as significant a change as reflected by smaller increases in track volumes and edge weights. This implies that the spatial information present in a 1 mm isotropic dataset does not possess sufficient quality to successfully infer and quantify brainstem white matter pathways.

Two notable exceptions to this trend are the SCP and ML streamline pathways. In particular, for the SCP streamline pathways,a 35% reduction in track volume and a 61% reduction in edge weight was observed as resolution was decreased from 0.333 to 1 mm. These reductions are likely to be due to the reduction in the surface area of the waypoint ROI used to infer these pathways. The waypoint ROI for the 0.333 mm dataset was 581.68 mm<sup>2</sup> , while for the 1 mm dataset its surface area was reduced to 300.24 mm<sup>2</sup> . Surface area reduction for this ROI did not result from poor registration of the ROI from 0.333 to 1 mm dataset, but rather from the fact that a number of voxels along the axial perimeter of the tissue had low diffusion signal due to volume averaging. As a result, following tensor fit, these voxels were assigned low FA values (less than 0.1), which in turn degraded the axial perimeter of the tissue ROI. During the registration process the voxels were classified as background rather than tissue, which resulted in their exclusion from the ROI. This ROI was affected by the volume averaging effects because it contained tissue within the posterior extent of the pons where the MR signal is naturally lower due to presence of many gray matter nuclei.

Similarly, volume averaging effects present at lower acquisition resolutions also affected the size of the ROIs used to delineate ML pathways. In the 2 mm dataset, ML ROIs were less than half the size of the 0.333 mm ROIs (**Table 2**: 0.333 mm ROIs had surface areas of 78.22 and 113.47 mm<sup>2</sup> , while 2 mm ROIs had surface areas of 27.71 and 55.43 mm<sup>2</sup> ). The ML pathway is a fairly small white matter pathway surrounded by gray matter nuclei. As spatial resolution is decreased, surrounding gray matter results in lower FA values due to volume averaging. This effect degrades the tissuetype contrast on the FA image and makes it difficult to identify the ML pathway, reflected by the smaller ML ROIs in our study. Smaller ROIs are a likely contributor to the decreased tract volumes and edge weights observed in our study. We found that the greatest decrease in both tract volume and edge weight was present as acquisition resolution was decreased from 1 to 2 mm (**Table 3**). Specifically, in the 1 mm dataset, tract volume and edge weight were very close to those observed in the 0.333 mm dataset. However, in the 2 mm dataset we found a 48% decrease in tract volume and a 95.6% decrease in edge weight as compared with the 0.333 mm dataset. These results indicate that the spatial resolution of an acquisition scan may greatly limit the size of white matter pathways that can be successfully inferred using diffusion tractography.

In summary, the results of the present study support the approach that sub-millimeter acquisition resolutions may be necessary to successfully infer and quantify white matter pathways within the brainstem using diffusion-weighted MRI. Previous studies employing 0.6 mm isotropic resolution have found a similar trend where crossing and branching fibers are resolved more efficiently at this high resolution, as compared with 1 and 2 mm acquisitions (Edlow et al., 2012). The present study used a higher acquisition resolution to further test effects of volume averaging in brainstem white matter tractography. Our results indicate that compared with results from highest resolution 0.333 mm isotropic dataset, pathways inferred in 1 and 2 mm datasets include a substantial amount of extraneous pathways that artificially inflate

## **REFERENCES**


tractography with multiple fibre orientations: what can we gain? *Neuroimage* 34, 144–155. doi:10.1016/j. neuroimage.2006.09.018


tractography results. Future studies employing even higher spatial and angular acquisition resolutions are necessary to determine the optimal resolution required to accurately visualize brainstem white matter pathways. In addition, histological quantifications of cross-sectional areas and volumes of the pathways should be used as a reference to fully assess volume averaging effects.

## **ACKNOWLEDGMENTS**

Anastasia A. Ford conceptualized the study, analyzed the data, and wrote the manuscript. Luis Colon-Perez acquired the data. William T. Triplett and Luis Colon-Perez developed the tractography and network analysis software. Joseph M. Gullett analyzed the data and edited the manuscript. David B. FitzGerald conceptualized the study and helped to write the manuscript. Thomas H. Mareci conceptualized the study and helped to write the manuscript. Anastasia A. Ford, Thomas H. Mareci, and William T. Triplett were supported through VA RR&D grant B6793-C. Luis Colon-Perez and William T. Triplett were supported through the State of Florida Brain and Spinal Cord Injury Research Trust Fund. Joseph M. Gullett and David B. FitzGerald were supported through VA RR&D grant B6698W. Thomas H. Mareci, Luis Colon-Perez, and William T. Triplett were supported through NIH grant R01 NS063360. Our thanks to Jennifer Adamson, BS, MBA, CCRP, and Benoit Giasson, Ph.D., both of the Center for Translational Research in Neurodegenerative Disease, Department of Neuroscience, University of Florida.

Dissociating the human language pathways with high angular resolution diffusion fiber tractography. *J. Neurosci.* 28, 11435–11444. doi: 10.1523/JNEUROSCI.2388-08.2008


MRI using mixture of Wischarts and sparse deconvulation. *Inf. Process. Med. Imaging* 20, 384–395. doi:10. 1007/978-3-540-73273-0\_32


*3D Sectional Anatomy*. Wien, NY: Springer.


Aldehyde fixative solutions alter the water relaxation and diffusion properties of nervous tissue. *Magn. Reson. Med.* 62, 26–34. doi:10.1002/ mrm.21977


IR-TSE, and MPRAGE: results and optimization. *MAGMA* 21, 121–130. doi:10.1007/s10334-008-0104-8

**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: 04 May 2013; paper pending published: 17 May 2013; accepted: 08 July 2013; published online: 24 July 2013.*

*Citation: Ford AA, Colon-Perez L, Triplett WT, Gullett JM, Mareci TH and FitzGerald DB (2013) Imaging white matter in human brainstem. Front. Hum. Neurosci. 7:400. doi: 10.3389/fnhum.2013.00400*

*Copyright © 2013 Ford, Colon-Perez, Triplett, Gullett, Mareci and FitzGerald. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*