# PROCEEDINGS OF THE INTERNATIONAL SCHOOL ON MAGNETIC RESONANCE AND BRAIN FUNCTION – XII WORKSHOP

EDITED BY : Federico Giove and Itamar Ronen PUBLISHED IN : Frontiers in Neuroscience and Frontiers in Physics

#### Frontiers Copyright Statement

© Copyright 2007-2018 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-88945-554-6 DOI 10.3389/978-2-88945-554-6

#### About Frontiers

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

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### Dedication to Quality

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

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
## PROCEEDINGS OF THE INTERNATIONAL SCHOOL ON MAGNETIC RESONANCE AND BRAIN FUNCTION – XII WORKSHOP

Topic Editors:

Federico Giove, Centro Fermi – Museo Storico Della Fisica e Centro di Studi e Ricerche Enrico Fermi, Fondazione Santa Lucia IRCCS, Rome, Italy Itamar Ronen, Leiden University Medical Center, Netherlands

Image: Daniele Mascali, based on data from S. Tommasin, D. Mascali et al. Scale-invariant rearrangement of resting state networks in the human brain under sustained stimulation. Neuroimage 179:570-581 (2018) doi: 10.1016/j.neuroimage.2018.06.006

In the last thirty years, Magnetic Resonance has generated a wide revolution in biomedical research and in medical imaging in general. More recently, the "in vivo" studies of the human brain were extended by new original ways to the dynamic study of function and metabolism of the human brain. The enormous interest in expanding the investigation of the brain is emphasizing the search for new NMR methods capable of extracting information of so-far obscure aspects of the brain function. In fact, many quantitative approaches have been proposed in order to complement the information obtained by functional MRI, and several multimodal and multiparametric approaches have been developed to exploit the information, either functional or structural, made available by the flexible contrast generation typical of MRI, and to combine it with complementary information. The XII workshop of the International School on Magnetic Resonance and Brain Function, held in Erice between 17 April and 6 May, 2016, was specially devoted to novel approaches aimed at better structural characterization of brain diseases, and at investigating frontiers MRI approaches to better understand the brain function. The papers included in this eBook offer a broad overview of the subjects covered during the Workshop, including applications of multiparametric MRI to neurological diseases, multimodal combination of MRI with electrophysiology, advanced methods for the investigation of brain networks and of brain physiology, and perspectives towards brain state reading.

Citation: Giove, F., Ronen, I., eds (2018). Proceedings of the International School on Magnetic Resonance and Brain Function – XII Workshop. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-554-6

# Table of Contents

Lauri J. Lehto, Aloma A. Albors, Alejandra Sierra, Laura Tolppanen, Lynn E. Eberly, Silvia Mangia, Antti Nurmi, Shalom Michaeli and Olli Gröhn

Jeroen Van Schependom and Guy Nagels

*91 Task-Related Modulations of BOLD Low-Frequency Fluctuations Within the Default Mode Network*

Silvia Tommasin, Daniele Mascali, Tommaso Gili, Ibrahim Eid Assan, Marta Moraschi, Michela Fratini, Richard G. Wise, Emiliano Macaluso, Silvia Mangia and Federico Giove

Ian D. Driver, Richard G. Wise and Kevin Murphy

*137 Reconstruction of Arm Movement Directions From Human Motor Cortex Using fMRI*

Seungkyu Nam and Dae-Shik Kim

# Editorial: Proceedings of the International School on Magnetic Resonance and Brain Function – XII Workshop

Federico Giove1,2 \* and Itamar Ronen3

1 MARBILab, Centro Fermi – Museo Storico Della Fisica e Centro di Studi e Ricerche Enrico Fermi, Rome, Italy, 2 Fondazione Santa Lucia IRCCS, Rome, Italy, 3 C. J. Center for High Field MRI Research, Leiden University Medical Center, Leiden, Netherlands

Keywords: International School on Magnetic Resonance and Brain Function (ISMRBF), MRI and fMRI, MRI techniques, multimodal MRI, neuroscience

**Editorial on the Research Topic**

#### **Proceedings of the International School on Magnetic Resonance and Brain Function – XII Workshop**

This Special Research Topic includes the proceedings of the XII workshop of the International School on Magnetic Resonance and Brain Function (ISMRBF), held between April 30th and May, 6th 2016 in Erice (Italy) at the premises of the Ettore Majorana Foundation and Centre for Scientific Culture. ISMRBF is a biennial event that gathers scientists whose research is focused on the crucial intersection between technological advances in MRI and emerging applications in neuroscience.

Indeed, many recent improvements in MRI instrumentation and techniques have been triggered by the increasing demands from a vast community of neuroscientists, which are exploiting MRI as one the most powerful tools for non-invasive investigation of human brain structure and function. Another long lasting topic covered by ISMRBF is the use of multimodal and multiparametric approaches. MRI is intrinsically multiparametric, thanks to the flexibility of the framework for contrast generation, that can be made sensitive to multiple biophysical phenomena. MRI can be also combined with other distinct modalities, such as electrophysiology, optical imaging and PET based molecular imaging. Fusion of information within- and across modalities has shown an excellent potential to complement the strengths and weaknesses of MRI as a stand-alone technique.

The papers included in this Research Topic offer a broad overview of the subjects covered by the Workshop. Mangia et al. contributed a paper that shows the potential of multiparametric MRI in application to neurological diseases, characterizing how Parkinson's disease and idiopathic rapid eye movement sleep behavior disorder are associated with tissue properties in subcortical regions. The authors exploited a battery of conventional (diffusion tensor imaging, DTI) and non-conventional metrics (adiabatic and non-adiabatic rotating frame relaxation), and concluded that rotating frame relaxation combined with functional connectivity may provide sensitive and specific signatures of progression of neurological diseases. The paper by Lehto et al. used similar MRI approaches on rats. In this case, MRI outcomes were correlated with myelin content in pharmacologically induced demyelinating lesions. Techniques based on non-adiabatic rotating frame relaxation showed good specificity and excellent sensitivity in the detection of lesions, and good quantitative correlation with actual myelin content. In comparison, the sensitivity reached by a battery of more conventional MRI approaches was lower. Authors concluded that non-adiabatic rotating frame relaxation is a valuable technique for the investigation of demyelinating diseases.

Edited and reviewed by: Ewald Moser, Medizinische Universität Wien, Austria

> \*Correspondence: Federico Giove<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biomedical Physics, a section of the journal Frontiers in Physics

Received: 07 February 2018 Accepted: 12 February 2018 Published: 06 March 2018

#### Citation:

Giove F and Ronen I (2018) Editorial: Proceedings of the International School on Magnetic Resonance and Brain Function – XII Workshop. Front. Phys. 6:18. doi: 10.3389/fphy.2018.00018

The contribution by Havlicek et al. expanded the multimodal approach to functional MRI. By combining multiple MRI and electrophysiological datasets from humans, non-human primates and cats, they discussed the generative model for BOLD (Blood-Oxygenation Level-Dependent) they recently introduced, and showed that it can explain the dynamic relationship between haemodynamic and neural events underlying the generation of the BOLD contrast. Hawsawi et al. dealt with multimodal combination of MRI and electrophysiological measurements in humans from the complementary point of view of patient safety. In their paper, the authors explored the safety issues related to simultaneous MRI and scalp or intracranial EEG, by identifying the physical phenomena that play a role in patient safety and reviewing the relevant literature. The authors concluded that risks can be managed, but that simultaneous invasive electrophysiology and MRI is characterized by a set of potential health hazards so complex, that careful case by case evaluation is needed.

The paper by Bednarik et al. is a retrospective study on patients affected by type 1 diabetes, investigating the connection between the pathology itself, impaired hypoglycaemia awareness, and cortical atrophy. Authors showed that type 1 diabetes is associated with gray matter atrophy compared to matched controls. The effect was larger in subjects with impaired hypoglycaemia awareness, likely because of increased exposure to hypoglycaemic episodes. Authors indicated that future studies should clarify whether the reported subtle reductions of gray matter volume are associated with cognitive impairment. The review by Van Schependom and Nagels discusses the relationship between neuroimaging and cognitive scores in multiple sclerosis (MS). The authors highlighted that improving the capability of neuroimaging-derived metric to assess pathology-related features can be useful not only for the characterization of the pathophysiological mechanisms underlying MS damage, but also to improve and speed-up the development of innovative therapies.

Tommasin et al. studied how the performance of a sustained cognitive stimulus affects the features of slow BOLD oscillations underlying functional connectivity within the Default Mode Network (DMN). They found that the stimulus reduced both amplitude of oscillation and functional connectivity within the DMN, but that the coupling between the reported changes was spatially heterogeneous and likely frequency specific. They concluded that the coupling between amplitude of BOLD fluctuations and functional connectivity is a potentially useful parameter to study brain networks. Bordier et al. presented a simulation work aimed at exploring optimal procedures for sparsification in graph theory analysis of topology of functional networks. Sparsification is often needed to simplify the processing and the interpretation of networks, however it is prone to over-pruning the network, and the identification of an optimal threshold is an open issue. The authors used synthetic data to show that percolation data-driven approaches can be used to identify a threshold that maximize the information of the pruned network.

The next couple of papers deal with MR approaches to study neurophysiology and neurochemistry. Sonnay et al. reviewed the applications of 13C-MR to brain energetics, and in particular how the combination of advanced MR techniques, mathematical models of metabolism and nutrients transport, and increasing magnetic fields can help the understanding the specific metabolic response to stimulation within astrocytes and neurons, as well as the metabolic interactions between these cells. Driver et al. studied an important issue regarding the measurement of brain oxygen consumption (CMRO2 ) by MR via the socalled calibrated-BOLD method, i.e., the usual isometabolic assumption during hypercapnia. The authors showed that their new approach, using the most relaxed and physiologically sensible hypothesis of linear dependence between a graded hypercapnic challenge and CMRO2 , allows more accurate estimates of basal CMRO2 as well as activation-related CMRO2 change.

Finally, Nam and Kim tackled the issue of decoding brain states by fMRI, by showing that brain reading can be achieved non-invasively even in the motor domain. Specifically, they demonstrated that directional arm movement can be predicted in humans exploiting simple linear encoding models.

The next workshop of the ISMRBF will be held in Erice, April 22–29 2018. The focus will be, as usual, on brain functional and microstructural MRI and compatible approaches. Information about the workshop can be found at the website http://ismrbf. marbilab.eu.

## AUTHOR CONTRIBUTIONS

All authors listed have made an equal contribution to the work, and approved it for publication.

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

Copyright © 2018 Giove and Ronen. 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) and the copyright owner 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.

## Multi-modal Brain MRI in Subjects with PD and iRBD

Silvia Mangia1 \*, Alena Svatkova2, 3, Daniele Mascali 4 , Mikko J. Nissi 5 , Philip C. Burton1 , Petr Bednarik 1, 3, Edward J. Auerbach1 , Federico Giove4, 6, Lynn E. Eberly 7 , Michael J. Howell 8 , Igor Nestrasil 2 , Paul J. Tuite8 and Shalom Michaeli 1

*1 Department of Radiology, Center for Magnetic Resonance Research (CMRR), University of Minnesota, Minneapolis, MN, United States, 2 Department of Pediatrics, University of Minnesota, Minneapolis, MN, United States, 3 Central European Institute of Technology (CEITEC), Masaryk University, Brno, Czechia, 4 MARBILab, Centro Fermi - Museo Storico Della Fisica e Centro di Studi e Ricerche Enrico Fermi, Rome, Italy, 5 Department of Applied Physics, University of Eastern Finland, Kuopio, Finland, 6 Fondazione Santa Lucia IRCCS, Rome, Italy, 7 Division of Biostatistics, University of Minnesota, Minneapolis, MN, United States, 8 Department of Neurology, University of Minnesota, Minneapolis, MN, United States*

#### Edited by:

*Kevin J. Black, Washington University in St. Louis, United States*

#### Reviewed by:

*Paolo Vitali, Fondazione Istituto Neurologico Nazionale Casimiro Mondino (IRCCS), Italy Anelyssa D'Abreu, Brown University, United States*

> \*Correspondence: *Silvia Mangia<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience*

Received: *20 April 2017* Accepted: *04 December 2017* Published: *19 December 2017*

#### Citation:

*Mangia S, Svatkova A, Mascali D, Nissi MJ, Burton PC, Bednarik P, Auerbach EJ, Giove F, Eberly LE, Howell MJ, Nestrasil I, Tuite PJ and Michaeli S (2017) Multi-modal Brain MRI in Subjects with PD and iRBD. Front. Neurosci. 11:709. doi: 10.3389/fnins.2017.00709* Idiopathic rapid eye movement sleep behavior disorder (iRBD) is a condition that often evolves into Parkinson's disease (PD). Therefore, by monitoring iRBD it is possible to track the neurodegeneration of individuals who may progress to PD. Here we aimed at piloting the characterization of brain tissue properties in mid-brain subcortical regions of 10 healthy subjects, 8 iRBD, and 9 early-diagnosed PD. We used a battery of magnetic resonance imaging (MRI) contrasts at 3 T, including adiabatic and non-adiabatic rotating frame techniques developed by our group, along with diffusion tensor imaging (DTI) and resting-state fMRI. Adiabatic T1ρ and T2ρ, and non-adiabatic RAFF4 (Relaxation Along a Fictitious Field in the rotating frame of rank 4) were found to have lower coefficient of variations and higher sensitivity to detect group differences as compared to DTI parameters such as fractional anisotropy and mean diffusivity. Significantly longer T1ρ were observed in the amygdala of PD subjects vs. controls, along with a trend of lower functional connectivity as measured by regional homogeneity, thereby supporting the notion that amygdalar dysfunction occurs in PD. Significant abnormalities in reward networks occurred in iRBD subjects, who manifested lower network strength of the accumbens. In agreement with previous studies, significantly longer T1ρ occurred in the substantia nigra compacta of PD vs. controls, indicative of neuronal degeneration, while regional homogeneity was lower in the substantia nigra reticulata. Finally, other trend-level findings were observed, i.e., lower RAFF4 and T2ρ in the midbrain of iRBD subjects vs. controls, possibly indicating changes in non-motor features as opposed to motor function in the iRBD group. We conclude that rotating frame relaxation methods along with functional connectivity measures are valuable to characterize iRBD and PD subjects, and with proper validation in larger cohorts may provide pathological signatures of iRBD and PD.

Keywords: rotating frame MRI, Parkinson's disease, iRBD, functional connectivity, DTI

## INTRODUCTION

Idiopathic rapid eye movement (REM) sleep behavior disorder (iRBD) is a sleep disorder associated with abnormal pontomedullary synuclein protein. iRBD evolves into Parkinson's disease (PD) or another synucleinopathy, e.g., dementia with Lewy bodies (LBD) or multiple system atrophy (MSA), in up to 90% of cases (Schenck et al., 1986, 1996; Boeve et al., 2001; Luk et al., 2012; Boeve, 2013; Bunzeck et al., 2013). Thus, iRBD may provide an opportunity to understand the evolution of neurodegeneration as well as a window in which to intervene with disease altering therapies. It is thought that the α-synuclein pathology in more caudal brainstem regions in iRBD ascends to the substantia nigra (SN) of the midbrain as the disease evolves into PD or other synuclein disorders (Boeve, 2013). Iron, which is present in higher quantities in the substantia nigra and other structures, may also play a role in the loss of dopaminergic neurons through its pro-oxidant effects (Zecca et al., 2008; Kaur et al., 2009). In addition, several studies have demonstrated abnormalities in structural integrity and functional connectivity in patients with iRBD (Hanyu et al., 2012; Ellmore et al., 2013).

Magnetic resonance imaging (MRI) offers several noninvasive contrasts that are useful in evaluating the pathological hallmarks of patients with iRBD, and any MRI findings in iRBD patients can be compared with those in individuals with PD. Previous studies have focused on evaluating structural approaches such as diffusion tensor imaging (DTI) and voxelbased morphometry (VBM) in patients with PD. This has included work that showed changes in DTI fractional anisotropy (FA) (Menke et al., 2009, 2010; Vaillancourt et al., 2009; Peran et al., 2010; Du et al., 2011) and mean diffusivity (MD) (Kamagata et al., 2012) in the substantia nigra of patients with PD. Others have shown that DTI may be useful in differentiating PD from other degenerative parkinsonian disorders and may also have the potential to monitor for disease progression (Prodoehl et al., 2013). Additionally, VBM demonstrated nondopaminergic brainstem alterations in iRBD (Scherfler et al., 2011). VBM and DTI have also provided insights into non-motor features such as depression (Feldmann et al., 2008; Kostic et al., 2010), smell (Ibarretxe-Bilbao et al., 2010; Rolheiser et al., 2011) and cognitive changes (Ibarretxe-Bilbao et al., 2009; Song et al., 2011) in PD. Therefore, DTI and VBM may prove useful in evaluating iRBD as it evolves into a Parkinsonian condition. So far, it has only been possible to show discriminative properties of the MRI methods using fairly large numbers of patients and controls, and these methods have not yet revealed if they are able to characterize the clinical severity of PD (Gorell et al., 1995; Martin et al., 2008). Therefore, our group previously employed novel microstructural imaging methods to map the SN (Michaeli et al., 2007; Nestrasil et al., 2010), as well as other brainstem and subcortical regions, in PD (Karagulle Kendi et al., 2008; Tuite et al., 2012).

In our current work, we employed several MRI modalities that are sensitive to various brain tissue properties, such as microstructural integrity, iron loads, and functional connectivity. Namely, we utilized novel rotating frame relaxation mapping methods developed in our group including adiabatic T1ρ

(Michaeli et al., 2006), T2ρ (Michaeli et al., 2004) and RAFF4 (Relaxation Along a Fictitious Field in the rotating frame of rank 4) (Liimatainen et al., 2010, 2011, 2015), together with DTI and resting-state functional MRI (rsfMRI). Our previous studies conducted in PD (Michaeli, 2007; Michaeli et al., 2007; Nestrasil et al., 2010; Tuite et al., 2012) and other brain conditions (Sierra et al., 2008; Jokivarsi et al., 2009; Mangia et al., 2013; Satzer et al., 2015) have shown that the proposed rotating frame methods provide information on different relaxation mechanisms in the human brain. T2ρ is sensitive to diffusion and exchange of water protons in environments with different local magnetic susceptibilities and reflects iron content with higher sensitivity than conventional T2 MRI (Michaeli et al., 2005, 2007). T1ρ provides an indication of neuronal loss and can be used to assess PD nigral degeneration with a higher sensitivity than is possible with conventional T1-based MRI (Michaeli et al., 2006, 2007). RAFF4 methodology, which is derived from adiabatic T1ρ and T2ρ (Liimatainen et al., 2010, 2011, 2015), has high sensitivity to myelin content (Satzer et al., 2015; Hakkarainen et al., 2016), and thus could provide information on microstructural integrity.

Rotating frame relaxations offer a conceptual advantage as compared to free-precession relaxation methods for characterizing tissue properties thanks to their inherent sensitivity to slower motional regimes (Michaeli et al., 2008). Indeed, while water molecules undergo a large variety of motional regimes in tissue, the motions which are most sensitive to the intricate nature of tissue microstructure and composition are those in intermediate-slow regimes. Rotating frame relaxations during frequency swept pulses offer additional practical advantages as compared to continuous-wave T1ρ, including a minimized sensitivity to B1 distributions and a capability to simultaneously tune the contrast to multiple effective frequencies (Mangia et al., 2009). Although not common for clinical studies of neurological disorders, rotating frame relaxations with frequency swept pulses were proven as robust and sensitive methods for a variety of in vivo investigations conducted by our group (Michaeli et al., 2007; Sierra et al., 2008; Jokivarsi et al., 2009; Nestrasil et al., 2010; Griffith et al., 2011; Liimatainen et al., 2012; Tuite et al., 2012; Mangia et al., 2013; Satzer et al., 2015; Hakkarainen et al., 2016) and others (Andronesi et al., 2014; Casula et al., 2017; Okuaki et al., 2017).

Our previous applications of mapping adiabatic rotating frame relaxations in humans had been limited to singleslice acquisitions performed on a non-clinical 4 Tesla scanner (Michaeli et al., 2007; Nestrasil et al., 2010; Tuite et al., 2012; Mangia et al., 2013). With the present work, we piloted the extension of these methodologies to multi-slice acquisitions on a clinical 3 Tesla scanner, and we thus embarked on a small cross-sectional study of 10 healthy controls, 9 patients with iRBD, and 10 patients with PD. Rather than focusing on the endpoint of progression to a clinical PD diagnosis for iRBD subjects, the ultimate goal is to use 3 T MRI for characterizing the spectrum of abnormalities resulting from inexorable synuclein deposition and other pathological features in both PD and iRBD. Our study is a step forward for establishing biomarkers which can characterize PD and iRBD, and ultimately the development of parkinsonian syndromes from iRBD subjects.

## METHODS

## Subjects

PD and iRBD subjects were recruited from the University of Minnesota Movement Disorders Clinic. Ten healthy controls, 9 iRBD patients and 10 individuals with PD were recruited and underwent the MRI protocol. However, all data from one iRBD subject and one PD subject were entirely discarded due to major motion observed during the acquisition. All the other subjects successfully completed the MRI protocol (**Table 1**), except one control subject who did not finish the DTI acquisitions, and one PD subject who presented excessive motion during the restingstate acquisition.

Exclusion criteria included history of stroke, seizures, neurosurgical procedures, arrhythmias, Montreal cognitive assessment (MoCA) scores lower than 25, and incompatibility with MR safety criteria. This study was carried out in accordance with the recommendations of The Code of Federal Regulations, Institutional Review Board, with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Institutional Review Board: Human Subjects Committee of the University of Minnesota. Clinical assessments were performed on anti-parkinsonian medications (for subjects on such treatment) and included the National Institute of Neurological Disorders and Stroke (NINDS) Common Data Element PD medical and family history, demographic information, the Hamilton Anxiety and Depression Rating Scales (Hamilton, 1959, 1960, 1980), United States version of the Parkinson disease questionnaire (PDQ-39) that reflects upon quality of life (Bushnell and Martin, 1999), the Unified Parkinson Disease Rating Scale (UPDRS), the Schwab and England, and Hoehn and Yahr (Fahn and Elton, 1987; Goetz et al., 2007, 2008; Dibble et al., 2010). In particular, we also used the UPDRS part III score to characterize our patient populations, as UPDRS part III represents an objective evaluation of motor symptoms. We assessed for features of iRBD with the validated REM sleep behavioral disorder questionnaire-Hong Kong (RBDQ-HK, English version) (Li et al., 2010) in all subjects. For cognitive function, we included the MoCA and the CogState (Nasreddine et al., 2005; Hoops et al., 2009; Gagnon et al., 2010).

## Imaging Protocol

Studies were performed on a 3 T Siemens Prisma system using a 32-channel RF receive coil. Standard B0 shimming were performed at the beginning of the imaging session to minimize B0 inhomogeneities.

Adiabatic T1ρ, T2ρ and non-adiabatic RAFF4 measurements were collected from 30 AC-PC aligned slices covering the brainstem and basal ganglia using a segmented GRE readout with 4 segments, voxel size = 1.6 × 1.6 × 3.6 mm3 , GRAPPA= 3, TE = 3.18 ms and TR = 2 s. For adiabatic relaxation measurements, hyperbolic secant (HS) pulses were used with BW=1.6 kHz, pulse duration Tp = 6 ms, peak power ω max 1 /(2π) = 800 Hz, 5 acquisitions with number of pulses = 0, 4, 8, 12, 16 using MLEV4 phase cycling (Michaeli et al., 2004, 2005, 2006, 2008). In T2ρ acquisitions, magnetization was on -z prior to readout. For RAFF4, Tp was 4.52 ms for one P-packet. The number of P-packets was 0, 4, 8, 12, 16 with ω max 1 /(2π) = 327 Hz. RAFF4 acquisitions were separately conducted without and with a global inversion (achieved by the same HS pulse used for T1ρ) prior to the train of P-packets, to produce sampling of the recovery and decay curve, respectively, to the steady-state magnetization induced during RAFF4.

T1-weighted, T2-weighted images, DTI and rsfMRI were collected covering the whole brain. For T1-weighted acquisitions MPRAGE with TR = 2150 ms, TE = 2.47 ms, TI = 1,100 ms, voxel size = 1 × 1 × 1 mm3 , flip angle = 8 degrees and GRAPPA = 2 was used. T2-weighted images were collected using SPACE sequence, voxel size = 1 × 1 × 1 mm3 , TE = 147 ms and GRAPPA = 2. Two DTI datasets were acquired with different phase encoding (i.e., anterior-posterior and posterior-anterior) utilizing TR = 2,820 ms, TE = 72.6 ms; multi band (MB) = 4 (https://www.cmrr.umn.edu/multiband), 7 non-diffusion weighed (b0) images, and 93 diffusion weighted images with b-value of 750 s/mm2 (47 images) and 1,500 s/mm2 (46 images) and voxel size = 1.8 × 1.8 × 1.8 mm3 . Finally, rsfMRI data was obtained using gradient echo Echo Planar Imaging (EPI) sequence, TR = 900 ms, multi band (MB) = 4; TE = 30 ms; voxel

*Characteristics of subjects whose data underwent further analyses. All subjects completed the whole MRI protocol, except one control subject who did not finish the DTI acquisitions. One PD subject was later excluded from functional connectivity analyses, because the resting-state functional scans exceeded the established criteria for acceptable motion.*

size = 3 ×3 × 3 mm3 , matrix size = 64 × 64, 48 AC-PC aligned slices with interleaved slice acquisition, 502 volumes.

## Image Analysis and ROI Definition

Brain segmentation of the T1-weighted scan was carried out with FreeSurfer (FS) version 5.3.0 (http://surfer.nmr.mgh. harvard.edu) installed on the supercomputer at the Minnesota Supercomputing Institute. Automated processing has been described previously (Fischl et al., 2002). A trained operator visually inspected each subject's data to ensure accuracy of the segmentation.

The following region of interest (ROI) masks were derived from FS automatic labeling of each subject's brain anatomy: caudate nucleus, pallidum, putamen, accumbens, amygdala, thalamus, hippocampus, and brainstem. With the exception of the brainstem, separate left and right hemisphere ROIs were generated for each structure. The brainstem was then manually separated to midbrain and pons regions. Additionally, left and right ROIs encompassing the substantia nigra reticulata (SNr), and separately the substantia nigra compacta (SNc), were manually drawn on each subject's T2-weighted image using itkSNAP software v.3.2.0 (Yushkevich et al., 2006) following the methodology described in our previous paper (Nestrasil et al., 2010). Representative examples of ROI selections are shown in **Figures 1**, **2**.

Rigid-body motion correction was performed on adiabatic T1ρ, T2ρ and RAFF4 measurements using the MCFLIRT algorithm in FSL (Jenkinson et al., 2012) with default options (Jenkinson et al., 2002). For T1ρ and T2ρ measures, all acquired images were firstly co-registered to the first acquisition and then to the averaged volume. The same algorithm was also applied on the two separate decay and recovery RAFF4 datasets to minimize the effect of movement between the two acquisitions. Relaxation time constant maps were calculated with MATLAB (MathWorks, Natick, MA) using in-house written routines and the Aedes software package (http://aedes.uef.fi). In particular, we used a 2-parameters non-linear fitting for estimating M0 and relaxation time in T1ρ and T2ρ acquisitions, and a 4-parameters non-linear fitting for estimating M0, Mss (steady state value of magnetization), -Mz (initial magnetization value measured from the negative hemisphere) and relaxation time for RAFF4 acquisitions. The FreeSurfer BB-register algorithm (Greve and Fischl, 2009) was utilized to align the first acquired image of T1ρ, T2ρ, and RAFF4 measurements to the FS-preprocessed T1 weighted images. The transformation matrix was then applied to the parametric maps utilizing trilinear interpolation. The registration accuracy was visually verified. The quality of BBregister initialization was additionally quantified for the adiabatic relaxation maps by using the quality assessment value provided by the algorithm.

DTI datasets with opposite phase encoding were utilized to reduce susceptibility artifacts and field inhomogeneities using FSL TOPUP (Andersson et al., 2003). After motion and eddycurrent correction (Andersson and Sotiropoulos, 2016), data were skull-stripped (Smith, 2002) and BB-register was applied to co-register DTI data to the T1-weighted image (Greve and Fischl, 2009). The parameters utilized for BB-register matched those used for motion-correcting the relaxation measurements, except that the first volume of relaxation measures was replaced by the b0 image in order to calculate the transformation matrix for DTI. The accuracy of registration was verified visually. The DTIFIT tensor model in FSL (Jenkinson et al., 2012) was fit to the co-registered DTI data to generate fractional anisotropy (FA) and mean diffusivity (MD) maps.

T1ρ, T2ρ, RAFF4, FA, and MD values of co-registered images were averaged across voxels within each ROI for each subject. Where applicable, the masked MRI ROIs extracted from the left and right hemispheres were combined before averaging. Representative examples of relaxation and DTI maps are shown in **Figure 3**.

## Resting State Functional Connectivity Pre-processing

Resting-state data were pre-processed and analyzed using AFNI (Cox, 1996), FS, and in-house routines based on Matlab R2013b. Data underwent the following pre-processing pipeline: (1) removal of first 10 volumes to allow for signal stabilization and for subjects to adapt to the new condition, (2) despiking (AFNI, 3dDespike), (3) slice-time correction (AFNI, 3dTshift), (4) realignment to the first EPI volume (AFNI, 3dvolreg), (5) coregistration of the first EPI volume to the FS-preprocessed T1-weighted image (FS, bbregister), and (6) normalization to 2 × 2 × 2 mm3 MNI space (FS, mri\_vol2vol). To minimize data resampling, steps five and six were applied simultaneously, concatenating the two transformations (FS, mri\_vol2vol). In addition, a general linear model was constructed to regress out baseline drift, motion and spurious variance from MNI normalized data (AFNI, 3dTproject). The following regressors were included in the model: a second order Legendre polynomial, a basis of sine and cosine to model frequencies outside the band 0.01–0.1 Hz, the six estimated motion parameters and their first derivative, the first five eigenvectors from the time series inside WM (after two voxels erosion; FS mask) and the first five eigenvectors from the time series inside CSF (after one voxel erosion; FS mask), following the aCompCor approach (Behzadi et al., 2007). To further reduce the impact of motion, data were censored by removing time points with more than 0.33 mm of displacement (estimated as the Euclidian norm of motion parameter derivatives) along with each previous time point. Censoring was applied during the regression step, removing time points from both data and regressors. The censoring threshold of 0.33 mm was selected so that at least 20 degrees of freedom (DOF) remained in the data. One PD subject could not meet this criterion for any reasonable choice of the censoring threshold and was therefore excluded from any further connectivity analysis, leaving a total of 10 controls, 8 iRBD and 8 PD subjects. The remaining DOF did not differ among groups (controls: 51.5 ± 7.9, RBD: 46.4 ±11.2, PD: 48.1 ± 7.8; ANOVA p = 0.48). Finally, conventional smoothing was not applied to the data since functional outcomes were not compared voxel-by-voxel but at the ROI level. However, since intersubject variability in data smoothness, arising from different amounts of interpolation due to motion correction and/or spatial normalization, directly impacts local connectivity estimation, we

FIGURE 2 | Substantia nigra ROIs (SNr, in blue, and SNc, in red) overlaid on T2-w images and adiabatic T1ρ , and T2ρ maps from one representative control subject. Images are shown either without (top row) or with (bottom row) the overlay of the ROIs.

set the effective data smoothness to an isotropic Gaussian FWHM of 4 mm (Maximo et al., 2013). This step was accomplished with the AFNI program 3dBlurToFWHM with the non-default option to avoid local smoothness estimation (i.e., whole-brain estimation only), in order to preserve the neuronal related smoothness.

#### Movement Assessment

Framewise displacement (FD), as defined in Power et al. (2012), was computed to index the subject's movements during the resting-state scan. There was no significant difference in FD between controls and PD (controls: 0.15 ± 0.06 mm; PD: 0.18 ± 0.04; control vs. PD p = 0.10), however, iRBD subjects moved more than controls (iRBD: 0.22 ± 0.05; iRBD vs. control p = 0.004) and PD subjects (iRBD vs. PD p = 0.10). To partially take into account this issue we included FD as a nuisance covariate in all statistical tests of the functional metrics.

#### Metrics

Two functional metrics were calculated: seed-based functional connectivity (FC) and regional homogeneity (ReHo).

#### Seed-Based FC

Time series within each of the ROIs described previously were averaged, yielding a single time series that served as seed. The correlation (Pearson's r) between the seed time series and the time series at each other voxel was then computed to quantify FC between the seed and the rest of the brain. To take into account the different number of DOF across subjects (due to the censoring procedure), correlation values were transformed in standard scores (Z) dividing the Fisher transformed r-value by the standard error of the Fisher's distribution:

$$Z = \operatorname{arctanh}(r) \* \sqrt{DOF - 3}.\tag{1}$$

FIGURE 3 | Representative maps of T1ρ , T2ρ , RAFF4, FA, and MD from one each of a control, iRBD and PD subject. ROIs of interest visible on the selected slices are also shown superimposed on the T1-weighted (T1w) images, including thalamus (beige), pallidum (dark blue), putamen (light blue), caudate (red), and accumbens (yellow).

For each subject and ROI, the seed-based FC was averaged within the functional network identified in the healthy control group. Specifically, the functional network was defined by a one-sample, one-tail, t-test on the seed-based FC maps of the healthy controls (p < 0.0005, minimum cluster size 20). We refer to this aggregate value of seed-based FC as network strength.

### ReHo

Regional homogeneity was implemented as a metric to assess local connectivity (Zang et al., 2004). For each voxel in the brain, ReHo was computed as the Kendall's coefficient of concordance among the voxel time series and the time series of its 18 nearest neighbors (AFNI, 3dReHo). Resulting brain maps were then transformed into z-scores by subtracting the whole-brain mean and dividing by the whole-brain standard deviation. Aggregate ReHo values were finally extracted by averaging across voxels in each ROI previously described.

## Statistical Analyses

Subject age, clinical characteristics, and movement (FD) were summarized by group with means and standard deviations (SD) and compared using linear models that allowed for groupspecific variances and provided Holm-corrected p-values for multiple comparisons among the groups (Holm, 1979); sex distribution was compared using Fisher's exact test. The average relaxation (T1ρ, T2ρ, and RAFF4), diffusion (FA, MD) and functional connectivity (network strength and ReHo) parameters were analyzed per ROI. For each of the 11 ROIs separately, each MR parameter was summarized by group (PD, iRBD, controls) using means and standard deviations, and compared across groups using linear models which allowed group-specific variances; we fit models with and without correction for age. Per ROI, type I error was controlled with Holm's adjustment for the 2 comparisons (PD vs. control and iRBD vs. control) and false discovery rate correction (Benjamini and Hochberg, 1995) was used for multiple testing of the 7 MRI outcomes. Results are presented both without and with the multiple-comparisons and the multiple-testing corrections (**Figure 4**, Supplementary Table 1).

## RESULTS

Individuals enrolled in this study spanned from middle-age to older adults (**Table 1**). The age of PD and iRBD subjects was 66 years on average for both groups, while the healthy controls were on average 9 years younger (p < 0.02, corrected, for each of PD

and iRBD vs. controls); the sex distribution was similar across all 3 groups. Disease duration spanned from 1 month to 4.5 years in the PD subjects, and from 1 to 15 years in the iRBD patients. As expected, PD subjects had significantly higher (worse) UPDRS scores and UPDRS-III scores than iRBD subjects (p = 0.00003 and p = 0.00009, respectively), but their MoCA values were comparable.

The adiabatic relaxation maps were well-registered to the anatomical images, with an average value of quality assessment of 0.21 ± 0.04 (normal range is between 0 and 1.1; 0 represents perfect registration initialization). The range of adiabatic T1ρ (∼130–180 ms), T2ρ (∼55–90 ms) and RAFF4 (∼ 240–410 ms) measured at 3 T from the gray matter structures of interest were in good agreement with previous studies conducted at 4 T (Michaeli et al., 2007; Nestrasil et al., 2010; Liimatainen et al., 2015). In addition, rotating frame relaxation maps were highly reproducible among subjects. Indeed, between-subject variations among healthy controls were 3% for adiabatic T1ρ, T2ρ, and RAFF4 from the regions of interest. Diffusion parameters were more variable than relaxation parameters, with between-subjects variations of 9 and 4% for FA and MD, respectively. Betweensubjects variations were considerably higher for functional connectivity parameters, namely 22 and 35% for network strength and ReHo, respectively.

Rotating frame relaxation parameters tended to be higher in patients than in controls for most of the ROIs (**Figure 4**, Supplementary Table 1). When age was not used as a covariate, significant Holm-corrected group-differences between PD patients and controls were observed for T1ρ in the amygdala (p = 0.0001), hippocampus (p = 0.03), pons (p = 0.011), putamen (p = 0.014), substantia nigra compacta (p = 0.018), and thalamus (p = 0.013), while between iRBD and controls T1ρ differences occurred in the amygdala (p = 0.012), hippocampus (p = 0.038), midbrain (p = 0.016), pons (p = 0.011), and thalamus (p = 0.013). Group differences between PD and controls were also observed for T2ρ in the amygdala (p = 0.033) and thalamus (p = 0.008), while none of the RAFF4 findings were significant. Holm-corrected group differences between iRBD and controls occurred for T2ρ and RAFF4 in the midbrain (p = 0.024 and p = 0.005, respectively). On the other hand, diffusion parameters exhibited only one Holm-corrected group difference, namely MD in SNr of iRB vs. controls (p = 0.001). Network strengths were smaller in patients than in controls for most regions. Holm-corrected group differences were significant for the network of the accumbens of both patient groups vs. controls (p < 0.013), and for the networks of the amygdala and SNr of PD vs. controls (p = 0.040 and p = 0.028, respectively). Finally, differences in ReHo were observed only between PD and controls in the SNc (p = 0.021) and SNr (p = 0.007).

Including age as a covariate demonstrated that most of the above mentioned differences were largely accounted for by age (**Figure 4**, Supplementary Table 1), especially for T1ρ outcomes. A lengthening of T1ρ in the amygdala of PD patients as compared to controls was however still robustly detected with Holmcorrection after age-adjustment (p = 0.005). PD patients also maintained a significantly longer T1ρ in the SNc (p = 0.046) as compared to controls. Between iRBD and controls, only trendlevel group differences were observed after age-adjustments, namely for T1ρ in the amygdala (p = 0.075), and for T2ρ and RAFF4 in the midbrain (p = 0.08 and p = 0.064, respectively). Network strength was still lower in the accumbens of iRBD patients as compared to controls (p = 0.003), but only a trend was observed between PD and controls (p = 0.085). Moreover, age-adjustments and Holm-correction maintained a significantly lower ReHo in the SNr of PD vs. control (p = 0.035), while only trend level differences were reached for the amygdala (p = 0.067).

Finally, after additionally correcting FDR for multiple testing (7 modalities), the longer T1ρ in the amygdala of PD vs. controls remained significant (p = 0.038), along with the smaller network strength of the accumbens of iRBD vs. controls (p = 0.018).

## DISCUSSION

In the current study we recruited an iRBD group in addition to a mild-moderately affected PD group with the ultimate goal of identifying similarities in MRI signatures of both diseases. Both patient groups were selected to have no diagnosis of dementia, with MoCA scores not less than 25 (**Table 1**, and Supplementary Figure 1). In addition, while iRBD subjects had findings on their UPDRS-III evaluations (**Table 1**, and Supplementary Figure 1), they were not thought to meet a clinical diagnosis of PD (Hughes et al., 1992) at enrollment. Motivated by the group-differences previously observed between PD subjects and controls when using adiabatic T1ρ and T2ρ (Michaeli et al., 2007; Nestrasil et al., 2010; Tuite et al., 2012), here we piloted a more extended imaging protocol which included T1ρ and T2ρ, RAFF4, DTI, and rsfMRI, and performed multi-slice acquisitions from the midbrain to characterize multiple subcortical brain regions rather than single slices of the brainstem as we used previously. Notably, mutli-slice mapping of T1ρ, T2ρ, and RAFF4 had never been implemented before, especially on a 3 Tesla clinical scanner. Here we demonstrate that 30-slices relaxation mapping is feasible with 5-min of acquisition time for each of T1ρ and T2ρ and 10-min for RAFF4. Thus the total acquisition time is similar to those employed for clinical MRI studies, and therefore such developments represent a critical step toward a future translation of the rotating frame methodologies to clinical trials. A multislice approach is needed for the purpose of obtaining more data that may allow an examination of regions involved in both motor and non-motor aspects of disease and that could be potentially useful in developing a prognostic biomarker (to predict whether iRBD patients are likely to develop PD) as well as a PD staging biomarker. With regard to the prognostic biomarker, based on clinical data, e.g., patients with greater cognitive impairment are more likely to have a different course than those without. If MRI methods can demonstrate pathology in sub-cortical regions involved in non-motor function, such as the amygdala or hippocampus, this information could be potentially useful in prognosis.

In principle, other MRI contrasts such as T1-Magnetization Transfer (MT) signal (so-called neuromelanin-sensitive scan Ogisu et al., 2013) and T1 relaxometry (Menke et al., 2009) could have been used, as they had been shown to reveal SN degeneration in PD. Such methods had been proven useful for visualizing and estimating the volumes of the SN structures. However, the neuromelanin-sensitive scan (Ogisu et al., 2013) is a qualitative approach based on signal intensity and it thus precludes a quantitative evaluation of SN tissue characteristics. On the other hand, T1-relaxometry is a quantitative method that allows quantification of tissue relaxation properties. However, the T1 measured in the SN did not reveal a difference between PD and controls (Menke et al., 2009), consistent with a lower sensitivity of free-precession relaxations to detect tissue degeneration as compared to rotating frame relaxations.

An important finding of the current study is that rotating frame relaxation parameters provided more robust characterizations of brain tissue as compared to DTI metrics, both in terms of observed coefficient of variations and in terms of their sensitivity to detect age-related and diseaserelated differences. The most prominent observation was that PD and iRBD groups share a lengthening of T1ρ in a region that is non-motor related and appears abnormal during cognitive impairment and depression, i.e., the amygdala. However, group differences vs. controls were clearly more pronounced in PD patients than in iRBD patients, as they remained robustly significant for PD patients also after ageadjustments and corrections for multiple testing, whereas they reached only trend level in comparisons of iRBD subjects vs. controls. Dysfunction in olfactory, cardiovascular, sleep, sensory and cognitive functions is often present during PD, and even can precede the motor symptoms and the disease onset (for example, see Jellinger, 2015). Therefore, the differences between PD and iRBD vs. controls in nonmotor areas detected by rotating frame T1ρ MRI further support the critical importance of considering non-motor abnormalities in PD and iRBD. The T1ρ results obtained in the amygdala of PD patients were overall consistent with previous findings of structural amygadalar deterioration in these patients (Bouchard et al., 2008; van Mierlo et al., 2015; Vriend et al., 2016). Since significant differences of amygadalar volumes were not observed in the present study (data not shown), and our PD subjects did not have significant mood or cognitive impairments, one may speculate that T1ρ MRI may detect subtle tissue changes in the amygadala which occur prior to structural changes and/or symptom manifestations. Interestingly, the lengthening of T1ρ in the amygdala as well as in SN regions of PD vs. controls, which is suggestive of neuronal degeneration in these areas, was accompanied by a deterioration of local functional connectivity as revealed by regional homogeneity.

Group differences were seen in rotating frame MRI parameters for other regions such as the thalamus, although such differences were largely explained by age. It is, however, still interesting to note that motor dysfunction in PD has been related to the change in activity of thalamic neurons in the motor circuits (Halliday, 2009). The slight lengthening of T1ρ in the SNc of PD patients as compared to controls is in agreement with previous studies (Michaeli et al., 2007; Nestrasil et al., 2010), although the finding was not robust enough to survive after correcting FDR for multiple testing, most likely due to the limited number of subjects and/or due to lack of on-to-one age matching of PD and controls (although age was taken into account in the statistical analyses). Importantly, adiabatic T2ρ was shown to pick up the difference in the SN of PD vs. control subjects in our previous studies conducted at 4 T (Michaeli et al., 2007; Nestrasil et al., 2010), but not in the current work. The differences in the outcomes among the current and previous studies could be attributed to the hardware used, i.e., higher magnetic field strength 4 T vs. 3 T (the lower field decreases sensitivity). Also, the multislice acquisition used in the current study may have impacted the efficiency of the contrast preparation module thus "diluting" the contrast as compared to the single slice used for the 4 T study. Importantly, in order to allow extended brain coverage, the spatial resolution used in this study was lower than previously used, thus leading to higher partial volume effects. Moreover, although the manual ROI segmentation of the SN was carefully performed, inaccuracies in placement are possible as there are unclear borders in T2-weighted images (Bolding et al., 2013). Differences between current and previous studies may also be due to different subject characteristics, with the most notable being a shorter disease duration of the current cohort (on average 2.5 years) as compared to our previous studies (on average 5 years in Michaeli et al., 2007; and 7 years in Nestrasil et al., 2010). Similar arguments apply for the pons, where T1ρ differences were previously observed (Tuite et al., 2012). In fact, in the present cohort we observed lengthening of T1ρ in the pons of PD vs. controls only when comparisons were not adjusted for age.

Finally, the long-distance connectivity as measured by network strength was significantly lower for the accumbens of iRBD subjects relative to controls, while such difference was not noted for the PD subjects. On the other hand, rotating frame relaxation and DTI parameters did not show evident abnormalities in the accumbens. Another putative signature of iRBD subjects was a longer RAFF4 and T2ρ in the midbrain, that could again be related to changes in non-motor features as opposed to motor function in the iRBD group.

In future studies, it would be important to follow the subjects for at least 3 years to identify iRBD patients who remain stable, those who progressed to PD, and those who progressed to LBD. Changes in nigrostriatal circuits are expected mainly in iRBD converted to PD, and to a lessser extent in iRBD converted to LBD. At the same time, it will be critical to fully characterize the iRBD and PD subjects with extensive behavioral and cognitive tests in order to fully understand the implications of the MRI findings. In addition, the lack in the present study of a complete age matching among controls and patients was a nuisance that had to be resolved with the inclusion of age-adjustments in the statistical analyses, and will have to be avoided in future investigations. Yet it allowed appreciating for the first time the interesting property that rotating frame relaxations may be exquisitely sensitive to age, a feature that may even be exploited in investigations that target aging itself. From an acquisition perspective, the rotating frame relaxation protocols also allow great deal of flexibility. Depending on the hypothesis to be tested, for instance if a focus on only hippocampus and amygdala are desired, one may design a specific rotating frame multi-slice acquisition protocol that achieves high spatial resolution with a reduced number of slices, while guaranteeing full coverage of the selected regions of interest. If acquisition time is limited, one may choose to collect only adiabatic T1ρ and resting-state fMRI, and skip entirely DTI which did not reveal sufficient sensitivity to detect group differences in this pilot investigation. On the other hand, other developments are currently underway to further optimize the acquisition times along with the intrinsic contrasts of the readout scheme. One promising direction is to implement magnetization preparation modules within a 3D readout sequence with virtually zero-echo times sequence, such as SWeep Imaging with Fourier Transformation (SWIFT) (Idiyatullin et al., 2006). Initial steps in this directions of research have already been proven successful in animal studies (Zhang et al., 2016).

Assessment of PD and iRBD subjects through combined microstructural and functional approaches may be a platform for comprehensive analysis of ongoing neurodegeneration. In essence, this approach may allow for staging of degeneration. In the present study, several outcomes showed trend-level differences in multiple regions, indicating the need of larger cohorts in future studies. Therefore, further studies on large cohorts of patients including longitudinal investigations of how iRBD advance to PD are necessary to characterize the ability of MRI metrics to be used as biomarkers of progression of iRBD into PD or other neurodegenerative diseases.

## CONCLUSION

We piloted for the first time the multi-slice mapping of rotating frame adiabatic T1ρ and T2ρ and non-adiabatic RAFF4 on a 3 Tesla scanner for human brain studies, which here involved control subjects along with non-demented iRBD and PD subjects. Rotating frame relaxations provided greater sensitivity to detect tissue abnormalities compared to conventional diffusion modalities that are commonly used to characterize tissue microstructure. The most prominent abnormality observed was a lengthening of T1ρ in the amygdala of PD subjects, indicating neuronal degeneration in a non-motor area. This finding was even more prominent than what was detected by T1ρ in the substantia nigra, i.e., the area related to the typical motor symptoms of PD. Significantly smaller network strength of the accumbens was also observed for iRBD subjects as compared to controls.

## AUTHOR CONTRIBUTIONS

SiM participated in design of the work, acquisition, analysis, interpretation of the data and preparing the manuscript. AS participated in analysis and interpretation of the data, and preparing the manuscript. DM participated in analysis and interpretation of the data, and preparing the manuscript. MN participated in acquisition and analysis of the data, and editing the manuscript. PCB participated in analysis and interpretation of the data, and editing the manuscript. PB participated in analysis and interpretation of the data, and editing the manuscript. EA participated in acquisition of the data, and editing the manuscript. FG participated in analysis and interpretation of the data, and editing the manuscript. LE participated in design of the work, analysis and interpretation of the data, and editing the manuscript. MH participated in design of the work, interpretation of the data, and editing the manuscript. IN participated in analysis and interpretation of the data, and editing the manuscript. PT participated in design of the work, interpretation of the data, and preparing the manuscript. ShM participated in design of the work, acquisition, interpretation of the data and preparing the manuscript.

## ACKNOWLEDGMENTS

We are grateful to the study coordinators and volunteers who participated in the study and thereby made this research possible. The authors acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for providing

## REFERENCES

resources that contributed to the research results reported within this paper. URL: http://www.msi.umn.edu. In addition, the authors would like to thank Karin Shmueli and Alejandra Sierra-Lopez for helpful discussions and feedback on the manuscript. Research reported in this publication was supported by the University of Minnesota Foundation, by the National Institutes of Health (Funding P41 EB015894, P30 NS076408, 1S10OD017974, UL1TR000114, and KL2TR000113), by the Academy of Finland (grant #285909) and by the Ministry of Education, Youth and Sports of the Czech Republic under the National Sustainability Programme II (project CEITEC 2020, LQ1601). This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant, agreement #691110 (MICROBRADAM). The content is solely the responsibility of the authors and does not necessarily represent the official views of the funding bodies.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins. 2017.00709/full#supplementary-material

with idiopathic rapid eye movement sleep behaviour disorder. Neurology 46, 388–393. doi: 10.1212/WNL.46.2.388

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

Copyright © 2017 Mangia, Svatkova, Mascali, Nissi, Burton, Bednarik, Auerbach, Giove, Eberly, Howell, Nestrasil, Tuite and Michaeli. 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.

# Lysophosphatidyl Choline Induced Demyelination in Rat Probed by Relaxation along a Fictitious Field in High Rank Rotating Frame

Lauri J. Lehto1, 2, Aloma A. Albors 1 , Alejandra Sierra1 , Laura Tolppanen3 , Lynn E. Eberly 4 , Silvia Mangia2 , Antti Nurmi 3 , Shalom Michaeli 2 and Olli Gröhn1, 2 \*

1 Department of Neurobiology, A.I. Virtanen Institute for Molecular Sciences, University of Eastern Finland, Kuopio, Finland, 2 Department of Neurobiology, Center of Magnetic Resonance Research, University of Minnesota, Minneapolis, MN, United States, 3 Charles River Discovery Services, Kuopio, Finland, 4 Division of Biostatistics, University of Minnesota, Minneapolis, MN, United States

#### Edited by:

Federico Giove, Centro Fermi, Italy

#### Reviewed by: Kai-Hsiang Chuang,

The University of Queensland, Australia Li-Wei Kuo, National Health Research Institutes, Taiwan

> \*Correspondence: Olli Gröhn<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 24 March 2017 Accepted: 13 July 2017 Published: 03 August 2017

#### Citation:

Lehto LJ, Albors AA, Sierra A, Tolppanen L, Eberly LE, Mangia S, Nurmi A, Michaeli S and Gröhn O (2017) Lysophosphatidyl Choline Induced Demyelination in Rat Probed by Relaxation along a Fictitious Field in High Rank Rotating Frame. Front. Neurosci. 11:433. doi: 10.3389/fnins.2017.00433 In this work a new MRI modality entitled Relaxation Along a Fictitious Field in the rotating frame of rank 4 (RAFF4) was evaluated in its ability to detect lower myelin content in lysophosphatidyl choline (LPC)-induced demyelinating lesions. The lesions were induced in two areas of the rat brain with either uniform or complex fiber orientations, i.e., in the corpus callosum (cc) and dorsal tegmental tract (dtg), respectively. RAFF4 showed excellent ability to detect demyelinated lesions and good correlation with myelin content in both brain areas. In comparison, diffusion tensor imaging metrices, fractional anisotropy, mean diffusivity and axonal and radial diffusivity, and magnetization transfer (MT) metrices, longitudinal relaxation during off-resonance irradiation and MT ratio, either failed to detect demyelination in dtg or showed lower correlation with myelin density quantified from gold chloride stained histological sections. Good specifity of RAFF4 to myelin was confirmed by its low correlation with cell density assesed from Nissl stained sections as well as its lack of sensitivity to pH changes in the physiological range as tested in heat denaturated bovine serum albumin phantoms. The excellent ability of RAFF4 to detect myelin content and its insensitivity to fiber orientation distribution, gliosis and pH, together with low specific absorption rate, demonstrates the promise of rotating frame of rank n (RAFFn) as a valuable MRI technique for non-invasive imaging of demyelinating lesions.

Keywords: RAFF, DTI, relaxation, myelin, demyelination, white matter damage, lysophosphatidyl choline

## INTRODUCTION

The quantitative assessment of myelin in the brain is highly important for both diagnostic and monitoring purposes of a variety of disorders, including multiple sclerosis (MS), traumatic brain injury (TBI), and stroke. However, conventional magnetic resonance imaging (MRI) approaches currently used to detect demyelination are suboptimal in their ability to detect myelin. To date, the identification of specific pathological processes involving demyelination by conventional MRI methods is not possible. This is because many different pathological processes result in similarlyappearing lesions in conventional MR, e.g., inflammatory demyelination, infection, stroke or tumor may all appear as bright lesions on T2-weighted images.

Advanced MRI modalities such as diffusion tensor imaging (DTI) (Schmithorst et al., 2002), magnetization transfer (MT) (Does et al., 1998), T1/T2 ratio (Glasser and Van Essen, 2011), multiexponential T2 (McCreary et al., 2009), and direct detection using ultra-short echo time (Wilhelm et al., 2012) have potential in assessing tissue organization and myelin. However, these methodologies just partially address the problem of myelin detection (Glasser and Van Essen, 2011; Nossin-Manor et al., 2013). For instance, although DTI is sensitive to tissue microstructure, it is not specific to myelin (Le Bihan and Johansen-Berg, 2012), because metrics derived from DTI data are influenced also by the underlying macroscopic level organization of myelinated fibers. On the other hand, MT has contributions from T1 in different water pools, and the results are dependent on experimental parameters such as off-resonance frequency and saturating B1 field strength. For example, the MT ratio (MTR) showed equally strong correlations with both the degree of myelin loss and the extent of axonal loss in postmortem samples (Schmierer et al., 2004). T1/T2 ratio has been shown to provide good correlation with myelin in normal gray matter (Glasser and Van Essen, 2011), but just like T1 and T2, it is influenced by the change of multiple different factors such as edema and paramagnetic ion concentration in the pathological brain (Pirko and Johnson, 2008; Beer et al., 2016). Furthermore, in a recent study by Dula and coworkers (Dula et al., 2010), myelin water fraction, an often cited multiexponential T2 relaxation metric, varied by almost a factor of 2 between two regions in the spinal cord while the myelin volume fractions differed by only approximately 12%. The proposed explanation for such observation relied on variation in microanatomy and intercompartmental water exchange. Overall, no currently available MRI technique or other noninvasive methods have proven capable of specifically signaling the myelin loss.

Recently, we introduced a novel MRI relaxation method entitled Relaxation Along a Fictitious Field (RAFF) (Liimatainen et al., 2010, 2011) in the rotating frame of rank n (RAFFn) (Liimatainen et al., 2015). With RAFFn, the fictitious fields are created by a non-adiabatic modulation of the RF amplitude and frequency. The particular case of the fictitious fields **H**n of rank n which are stationary in the rotating frame of rank (n-1) was detailed in (Liimatainen et al., 2015). We have shown that frequency modulated pulses operating in the nonadiabatic regime in the 1st rotating frame, and thus producing a large fictitious magnetic field in the 2nd rotating frame, can be used for generating novel MRI contrast in living samples. RAFFn is conceptually different from more conventional rotating frame MRI methods. With continuous wave spin lock and adiabatic T1ρ and T2ρ methods, relaxations are governed solely by longitudinal or transverse relaxations in the rotating frame. On the other hand, relaxations during RAFF2 comprise contributions from both T1ρ and T2ρ relaxation pathways. Similarly, RAFFn comprises both T1ρ(n) and T2ρ(n) relaxation channels. Moreover, RAFFn is also conceptually different from chemical exchange saturation transfer (CEST) and/or MT since RAFFn is sensitive to all contributing relaxation mechanisms including dipole-dipole interactions, diffusion in local field gradients and exchange between spins with different chemical shifts and the cross-relaxations. On the other hand, CEST and MT provide direct measure of the saturation transfer or exchange between off-resonance spins and water. We have shown that for adiabatic T1ρ the contribution of cross relaxations (origin of MT) between spins within the bandwidth of the adiabatic pulse is minor (Michaeli et al., 2008). The contribution of off-resonance saturation to RAFFn relaxation may originate also from side bands at frequencies defined by the duration of the RAFFn pulses (Liimatainen et al., 2015), a mechanism which may enhance sensitivity of RAFFn to multiple spin pools within the myelin component. These frequencies are determined by the assembling of the RAFFn pulses into P-packets, which refocuses the magnetization in the form of rotary echo (Solomon, 1959). Formation of the rotary echoes in the rotating frame during RAFFn pulses is an additional distinction from MT and CEST techniques.

RAFFn generates flexible MRI contrast with low specific absorption rate (SAR) and is thus particularly safe for human applications. We have demonstrated that by changing the rank of the rotating frame and the orientation of the fictitious field in the rotating frame of reference, RAFFn MRI contrast can be sensitized to various motional regimes. In particular, RAFF4 and RAFF5 are exquisitely sensitive to slow/ultra-slow motion characterized by the correlation times (τc) in the millisecond time scale. This is important for sensitive and specific myelin detection as the highly organized structure of myelin contains multicompartment water environments: myelinic, intra-axonal, and interaxonal. These pools are characterized by different molecular mobilities and T2 relaxation components which correspond to ultra-short (50µs–1 ms), short- (1–50 ms), intermediate (50 ms– 0.5 s), and long-lived T2 components (>0.5 s) (Does et al., 1998; Bonilla and Snyder, 2007). The short T2 component is attributed mainly to myelin water, whereas an ultra-short component is thought to arise from carbon-bound methylene protons. Extensive experimental and theoretical studies by Mefed and coworkers have demonstrated that relaxation measurements in high–rank rotating frames of T1ρρ and T1ρρρ (second and third frame longitudinal relaxations, respectively) allow the probing of slow-ultra slow motion with the characteristic correlation times up to 10−1–10−3 s (Mefed, 1999, 2001). Our previous analysis, in agreement with those by Mefed, also shows that the sensitivity of RAFF4 and RAFF5 is shifted to the motional regime in the millisecond time scale as compared to adiabatic T1ρ (Michaeli et al., 2006) which has a maximal sensitivity to the correlation times in the microsecond time scale (Satzer et al., 2015).

Utilizing the sensitivity of RAFF4 and 5 to slow molecular motion, we have recently embarked on measuring myelin content in the normal brain and in complex mucopolysaccharidosis type I (MPS-I) pathology including demyelination (Satzer et al., 2015; Hakkarainen et al., 2016). Importantly, the highest correlation between relaxation time constants and myelin content as assessed by quantitative histology was achieved with RAFF4 and RAFF5 techniques, as compared to other relaxation based contrasts (Satzer et al., 2015; Hakkarainen et al., 2016).

Lehto et al. LPC Induced Demyelination Probed by RAFF

In the present work, we further focused on intracranial lysophosphatidyl choline (LPC) injections in the rat brain that have been used previously to model focal demyelination lesions in MS (Waxman et al., 1979; Ransohoff, 2012). LPC injection causes demyelination in the white matter within 3-5 days with only mild inflammatory response, thus allowing a clean model to study demyelination. Another important feature of the model is that location of the white matter lesion can be freely chosen, thus allowing to study the influence of underlying macroscopic structure of the myelinated fibers on MRI contrasts. The optimal MRI contrast for detection of demyelination would be sensitive to the amount of myelin regardless of macroscopic level anisotropy.

The aim of this study was to characterize the ability of RAFF4 to detect demyelination processes in presence of different underlying tissue macrostructure and to compare that with some of the approaches currently used to detect white matter lesions in advanced human studies. To achieve these goals, we induced demyelination by injecting LPC in the rat corpus callosum (cc) and dorsal tegmental tract (dtg), measured RAFF4, DTI and MT parameters, and performed histology validation to assess the ability of the various MRI modalities to detect demyelination. Furthermore, influence of pH on RAFF4 contrast was evaluated using bovine serum albumin containing phantoms.

## MATERIALS AND METHODS

## Animal Model

Male Spraque Dawley rats (n = 21, Charles River, Germany; 300–350 g) were used in the experiments. Rats were group housed with a preserved 12 h light/12 h dark cycle and ad libitum access to food and water. All animal procedures were approved by the Animal Ethics Committee of the Provincial Government of Southern Finland, and conducted in accordance with the guidelines set by the European Commission Directive 2010/63/EEC.

All surgical procedures were done under inhalation anesthesia using 1.8–2.2% isoflurane in 30%/70% O2/N2O. To induce demyelinated lesions, stereotaxic injections of the LPC solution (volume of 1.5µl; concentration: 10 mg/ml; L-α-lysophosphatidylcholine from egg yolk; L-4129 Sigma-Aldrich, St. Louis, USA) were performed in selected areas of the rat brain, chosen based on their myeloarchitecture. The injections were placed either in the cc (n = 6 and n = 4, LPC and vehicle, respectively) at the stereotaxic coordinates: 0.4 mm posterior from bregma, 1.4 mm lateral from the midline, and 2.6 mm from the skull. For the dtg (n = 6 and n = 4, LPC and vehicle, respectively) the injection coordinates were at 6.3 mm posterior from bregma, 0.8 mm lateral from the midline, and 4.3 mm from the skull. The cc represents a major white matter bundle with the parallel organization of myelin bundles, while dtg is a typical example of more complex myeloarchitecture with a high density of heterogeneously oriented axons. Control animals underwent identical protocol but were injected with 1.5µl of vehicle solution [0.1 M sodium phosphate buffer (PBS)] instead of LPC.

## MRI

All the animals were imaged 3 days after injection, when there is already significant demyelination in this animal model, while inflammatory reaction and remyelination typically develops later (Waxman et al., 1979). All MR experiments were performed using a horizontal 7 T magnet system (Bruker Pharmascan, Entlingen, Germany).

The sites of injections were localized using T2-weighted fast spin-echo (FSE) images with the following parameters: TR = 4 s, echo spacing 12 ms, TEeff = 48.0 ms, necho = 8, FOV = 25.6 × 25.6 mm2 , matrix size = 256 × 256, number of slices = 20 and slice thickness = 0.75 mm. The imaging slice for RAFFn, MT and DTI (middle slice) was centered to the FSE slice next to the injection site to avoid including mechanical damage induced by the injection.

## RAFFn Technique

A detailed description of the RAFFn technique was presented elsewhere (Liimatainen et al., 2015). We chose to use RAFF4, as RAFF4 and 5 showed highest correlation to myelin content in normal rat brain in our previous study (Hakkarainen et al., 2016), and as RAFF4 relaxation is faster than RAFF5 thus requiring shorter pulse train duration than RAFF5. With RAFFn, the sine/cosine pulses were used for the modulation of the amplitude and frequency, respectively. The orientation of the effective field of RAFFn in the rotating frame of rank (n-1) relative to the Z (n−1) axis of quantization is defined by the angle α (n), as was previously detailed (Liimatainen et al., 2010, 2011, 2015). Here, the α (n) was set to 45◦ in each rotating frame of rank n-1, and the amplitude of the effective field Hn was maintained at the same level. During RAFFn pulses, the rotary echoes are generated when using four pulse elements assembled into a Ppacket according to the scheme PP−1PπP −1 π (Liimatainen et al., 2010). An instantaneous flip of the effective field Hn is indeed performed during each P-packet to ensure refocusing of **M** on the Z (n−1) axis and to form rotating frame rotary spin echo (Solomon, 1959). The peak RF amplitude of RAFF4 pulses was γB1 = 323 Hz. The time duration of each PP−1Pπ P −1 π packet defined as Tp = 4π/( √ 2ω max 1 ) was set to 4.525 ms. The signal intensity decay was measured by incremental pulse trains of Ppacket, with an inversion pulse to account for steady state. RAFF4 pulse train durations were 0, 27, 54, 81, and 108 ms.

Fast spin-echo (FSE) was used as a readout pulse sequence with TR = 4 s, TEeff = 8.3 ms, necho = 8, FOV = 32.0 × 32.0 mm2 , matrix size = 256 × 256, number of slices = 1 and slice thickness = 0.5 mm, leading to total acquisition time 21 min 20 s for one relaxation time constant map.

## Comparator Techniques: DTI and MT

MT measurements were conducted using the same FSE readout sequence as described above for RAFF4. The modified inversion MT protocol with two consecutive acquisitions was used (Mangia et al., 2011). First, the signal decay during off-resonance irradiation with **M** oriented along +Z was acquired, and, second, the signal recovery was measured when **M** was inverted along −Z. A square saturation pulse with γB1 = 200 Hz was placed 8 kHz off-resonance with an incrementing pulse duration of 0, 0.3, 0.6, 0.9, and 1.2 s. Total acquisition time was 21 min 20 s. T1sat, steady state magnetization MSS and magnetization in fully relaxed state, M0, were solved using pixel-by-pixel analysis with monoexponential decay and rise functions to the same steady state value, as described in Mangia et al. (2011). MTR was calculated as MTR = 1–MSS/M0.

For DTI, segmented spin-echo EPI was used with TR = 2 s, TE = 30.0 ms, nshots = 6, number of averages = 12, FOV = 26.5 × 18.0 mm2 , matrix size = 212 × 144, number of slices = 9, slice thickness = 0.5 mm, b = 1,000 s/mm2 , diffusion directions = 42 leading to total acquisition time of 1 h 55 min. Mean diffusivity (MD), fractional anisotropy (FA), radial and axial diffusivity (RD, AD) maps were calculated from DTI data.

## ROI Analysis with MRI

Relaxation time constant maps and parametric maps from MT and DTI data were calculated in MATLAB (MathWorks, Natick, MA). The regions of interest (ROIs) were hand-drawn using the Aedes software package (http://aedes.uef.fi). Six ROIs were drawn in the cc, three contralateral (1, 2, and 3) and three ipsilateral (4, 5, and 6) to the injection site (**Figure 3**). Mean values of ROI 4 further averaged over LPC or vehicle animals were used to evaluate the relative contrast (RC) of each of the different MRI metrics between LPC and vehicle injected animals. For metrics with lower intensity in vehicles, RC = [LPC–vehicle]/LPC ∗ 100%. For metrics with higher intensity in vehicles, RC = [vehicle–LPC]/vehicle ∗ 100%. This approach ensures that RC is comparable between techniques regardless of the direction of change between LPC and vehicle animals. Further, by dividing the difference with the higher value makes RC a conservative approach as it has a maximum of 100%. In addition, relative contrast-to-noise ratio (RCNR) was calculated so that RCNR = [higher value–lower value]/SD(mean values of ROI4 of the vehicles). To have a broader dynamic range of myelinated tissue, all 6 ROIs were used to calculate correlations between MRI and myelin content based on histological stainings as described below. For the dtg, an ROI covering the lesion was drawn and copied to the contralateral side to the corresponding anatomical location. For the vehicle injected animals, a 3-by-3 voxel ROI was drawn at the injection site and copied to the same anatomical location on the contralateral side. RC between LPC and vehicle injected animals were then calculated as for the cc, using the contralateral ROIs.

## Histological Procedures

Immediately after the MRI scans, all the animals were transcardially perfused with 0.9% NaCl (30 ml/min, 2 min, 4 ◦C), followed by 4% paraformaldehyde (PFA) solution in 0.1 M phosphate buffer (pH 7.4, 30 ml/min, 25 min, 4◦C). Fixed brains were removed from the skull, and post-fixed for 4 h in 4% PFA. Then, the brains were cryoprotected in 20% glycerol in 0.02 M potassium phosphate-buffered saline (pH 7.4) for 36 h, and frozen in dry ice. The frozen brains were stored at −70◦C until sectioning.

The brains were cut in five series of 30-µm thick coronal sections using a sliding microtome. The first series of sections was stored in 10% formalin, at room temperature. The series from second to fifth were stored in a cryoprotectant tissue-collecting solution (30% ethylene glycol, 25% glycerol in 0.05 M PBS) at −20◦C until further processing.

The first series of sections was stained with Nissl (thionin) to study cytoarchitecture, cell death and gliosis. The second series of sections was stained with gold chloride to assess the myeloarchitecture. For myelin staining, sections were mounted on gelatin-coated slides and dried at 37◦C. They were then incubated in a 0.2% gold chloride solution (HAuCl4•3H2O, G-4022 Sigma-Aldrich, St. Louis, USA) in 0.02 M PBS (pH 7.4) containing 0.09% NaCl for 7–8 h at room temperature in the dark. Then, the slides were washed twice for 4 min in 0.02 M PBS in 0.09% NaCl and placed in a 2.5% sodium thiosulfate solution for 5 min. After three 10 min washes in the buffer solution, sections were dehydrated through an ascending series of ethanol, cleared in xylene and coverslipped with DePeX (BDH, Laboratory Supplies, Dorset, UK).

## Histological Analyses

High-resolution photomicrographs of both myelin- and Nisslstained sections of the cc and the dtg were obtained using a light microscope (Leica DMRB, Wetzlar, Germany) equipped with a digital camera (DXM1200F, Nikon Instruments Inc., Japan). Three consecutive sections of these areas were analyzed covering a volume of 450µm that corresponded to the slice thickness in MRI.

The optical density on myelin- and Nissl-stained sections was quantified using ImageJ software (version 1.41 http:// rsb.info.nih.gov/ij/). For cc and dtg, ROIs were drawn corresponding in location, size and shape to the MRI ROIs. Optical density (OD) was averaged for each ROI over three consecutive sections to cover the volume analyzed in MRI. After conversion of the images to gray-scale, OD values were obtained from the ROIs. For correction of possible staining differences between sections and brains, OD values from healthy reference areas were obtained. For the cc lesions, white matter in the striatum was used as a reference and for dtg lesions, cc (in the sections of dtg) was used as a reference.

## Transmission Electron Microscopy (TEM)

One extra rat, injected with LPC as described above, was used to study the ultrastructure of the myelin sheaths in the corpus callosum 3 days after the LPC injection using TEM. T2-weighted images were acquired to verify the presence and consistency of the demyelinated lesion. After MRI, the rat was perfused using 0.9% NaCl (30 ml/min) for 2 min followed by 4% PFA (30 ml/min) in 4◦C for 50 min. The brain was removed from the skull, post-fixed in 4% PFA/1% glutaraldehyde overnight at 4◦C, and sectioned in 1-mm thick coronal sections in a brain matrix. One 1-mm thick section matching the level of the MRI analysis was selected based on T2-weighted images. The corpus callosum of this section was dissected into eight small portions from the level of the injection site to the same level on the healthy contralateral side. These corpus callosum samples were incubated with 4% PFA in 0.1 M cacodylate

buffer (pH 7.4) at 4◦C for an overnight. Then, the samples were rinsed in 0.1 M cacodylate buffer 5 min three times, followed by a post-fixation in 1% osmium tetroxide (OsO4) in 0.1 M cacodylate buffer for 2 h. After that, the samples were again rinsed in 0.1 M cacodylate buffer for 5 min three times. Then, the samples were dehydrated in ascending percentage of alcohols for 10 min, and the last incubation was propylene oxide twice for 10 min. The samples were then infiltrated with a 1:1 mixture of propylene oxide and LX-112 resin (Ladd Research Industries Inc., USA) for 1 h, followed by incubation in LX-112 resin overnight. The embedding was done with fresh LX-112 resin in molds. The polymerization of the samples was done in an oven at 37◦C for 24 h, and then at 60◦C for 48 h.

Once the samples were embedded in the resin, four semi-thin sections of 1µm were sectioned and stained with Toluidine blue. Toluidine blue stained sections were used to study the cyto- and myeloarchitecture of the samples, and to guide TEM imaging. Six ultrathin sections were cut for TEM and mounted in three different grips. Imaging was done in a transmission electron microscope JEOL JEM-2100F using magnifications between 250 and 10,000x. Photomicrographs were taken using a digital camera connected to the microscope.

## Bovine Serum Albumin Phantoms

Heat denaturated bovine serum albumin samples (BSA) were used to investigate influence of pH on RAFF4. Heat denaturated BSA has been shown to provide a good model to study influence of pH modulated proton exchanges as relaxation mechanism in tissue, while glutaraldehyde crosslinking blocks a significant portion of exchangeable NH-groups in the protein (Mäkelä et al., 2001). BSA (fraction V, Sigma Chemicals, St. Louis, MO) was dissolved in 0.1 M Tris (pH 7.4) to yield a 10% solution. The pH was adjusted to be 6.2, 6.6, 7.0, 7.4, and 7.8 in five BSA samples, respectively, using HCl or NaOH. BSA phantoms were exposed to 65◦C for 3 min for heat denaturation.

The pH phantoms were imaged using a horizontal 9.4 T magnet (Magnex Scientific Ltd., Abington, UK) interfaced to a Varian (Agilent) DirectDrive console (Agilent Technologies, Santa Clara, CA). RF transmission and signal reception were carried out using a single loop coil with 20 mm diameter. T2 weighted FSE MR imaging was used as a readout with the parameters as follows: TR = 4 s, TEeff = 55 ms, matrix size 256 × 256, FOV = 25.6 × 25.6 mm2 , 8 echoes with 8 ms echo spacing, initial TE = 10 ms, one slice with slice thickness 1 mm. The settings of RAFF4 were same as for in vivo studies described above. T2-maps were measured using double spin-echo with two

3 ms adiabatic full passage pulses and echo times of 5, 7, 15, 23, 31, 39, and 63 ms.

## Statistical Analyses

The Mann-Whitney U-test was used to compare LPC injected animals to vehicles. The contribution of myelinated axons and cell density to the MRI metrics was assessed using Pearson's linear correlation (two-tailed) between pooled ROI analysis results from MRI and histology of myelin and Nissl stained sections, respectively. ROIs from both LPC and vehicle injected animals were included in the analysis including contralateral and ipsilateral ROIs. Correlation was calculated separately for all data, cc data and dtg data. False discovery rate (FDR) correction was done to adjust p-values for inflation of error due to multiple testing.

FIGURE 2 | ROI analysis of MRI parameters in the corpus callosum, ipsilateral side. Mean and SD using RAFF4 (A), T1sat (B), magnetization transfer ratio, MTR (C), fractional anisotropy, FA (D), mean diffusivity MD (E), and axial and radial diffusivity, AD and RD (F,G) values obtained from LPC-induced lesion ROI 4 in the corpus callosum of LPC injected (n = 6) rats and from the corresponding ROI in the vehicle injected (n = 4) rats. Statistical significance; \* < 0.05 (FDR adj. for testing 7 MRI parameters), Mann-Whitney U-Test.

## RESULTS

## Corpus Callosum–Demyelination in a Major white Matter Track with Parallel Axons MRI

Representative parametric maps from relaxation, MT and diffusion data obtained from a rat with LPC injection in the cc are shown in **Figure 1**. The imaging slice was chosen to be 0.75 mm posterior from the injection site to avoid the influence of mechanical damage caused by the injection needle. The lesion was clearly visible in all quantitative MRI maps as well as in myelin and Nissl stained histological slices on the ipsilateral site close to the injection site. In all animals, the lesion was located within the cc, extending ipsilaterally 1.95 ± 0.24 mm (mean ± SD; based on T2-weighted images) from the midline and contralaterally 0.05 ± 0.15 mm from the midline in mediolateral direction, due to diffusion of the LPC along the tract.

In the ROI analysis, all MRI metrics showed statistically significant differences (p < 0.05) between lesions induced by LPC in the ipsilateral cc (ROI 4, in **Figure 3**) and the corresponding area in vehicle animals (**Figure 2**). The RC observed with RAFF4 was 19.7%. The RC of MTR and T1sat were 7.4 and 4.7%, respectively, thus exhibiting smaller RC than RAFF4. On the other hand, MD and FA showed the largest RC: 29.8 and 22.5%, respectively. Reduced MD and FA were explained by decrease in axial diffusivity (AD) by 41.3%, while radial diffusivity (RD) did not show a statistically significant difference, consistent with axonal damage and contradicting the typical pattern of diffusion changes often associated with demyelination (Song et al., 2005). While RC of DTI metrics were higher than that of RAFF4, RCNR

of RAFF4 was higher at 5.7 in comparison to 4.4, 1.8, and 2.6 of MD, FA and AD, respectively. Correspondingly, RCNR for MTR and T1sat was 2.3 and 1.8 respectively.

## Histology and Electron Microscopy

An example of myelin- and Nissl-stained section from the same animal as was presented for MRI in **Figure 1** is shown in **Figure 3**. Three days after LPC injection in the cc, a demyelinated area over the whole ipsilateral cc was observed (ROI 4 to 6) while the myelin content on the contralateral side appeared normal (ROI 1 to 3, **Figure 3A**). Also ROI 4, closest to the midline in the ipsilateral side still presented signs of demyelination (**Figure 3C**), which gradually decreased toward the contralateral side and ROI 3 (**Figure 3B**).

On Nissl-stained sections, a slightly increased cell density was observed on the demyelinated area (**Figure 3D**), which can be attributed to mild gliosis. The cell density was higher in ROI 6 and decreased toward the midline (ROI 4) and the contralateral side (ROI 3) (**Figure 3D**). Mild gliosis overlapped with the demyelinated area observed in myelin staining (**Figure 3A**). Gliosis was diffusely distributed; however, a high number of cells were found surrounding the blood vessels (**Figure 3F**). On the contralateral side to the injection site, the density and distribution of cells appeared normal (**Figure 3E**).

In order to understand the unexpected diffusion results in cc, samples from one additional animal underwent TEM. T2 weighted imaging confirmed that the location and extent of the LPC lesion was comparable to lesions in other cc injected animals (**Figure 4A**). Toludine blue staining showed clearly decreased staining intensity around the axons in the ipsilateral cc when compared to the contralateral cc, corroborating demyelination (**Figures 4B,C**). In TEM, axonal myelin sheaths appeared thick and tightly packed on the contralateral cc. Intra-axonal space appeared normal with only a few membranes or organelles, such as mitochondria, that could restrict diffusion along axons (**Figure 4D**). In the demyelinating ipsilateral cc, axons were either completely demyelinated or surrounded by thin and disorganized myelin sheaths (**Figure 4E**). In pockets between the myelin sheaths, vacuoles and myelin debris were observed around all the demyelinated axons, forming more boundaries along the axonal length.

## Dorsal Tegmental Tract–Demyelination in an Area with a Complex Myeloarchitecture MRI

Lesions with a width of 0.59 ± 0.18 mm were clearly visible in RAFF4, T1sat and MTR maps in the dtg (**Figure 5**). In ROI analysis of the dtg, RAFF4 and the MT MRI metrics showed a statistically significant difference (p < 0.05) between lesions

axial and radial diffusivity, AD and RD (F,G) of LPC-induced lesion in the dorsal tegmental tract. Representative examples of ROIs for analyzing lesions in the dorsal tegmental tract are shown on a grayscale RAFF4 map in (H). Arrowheads indicate the lesion site.

induced by LPC and the corresponding area in vehicle injected animals, while the diffusion metrics did not show a significant difference (**Figure 6**). The RC of RAFF4 was 26.2% while the RCs of MTR and T1sat were lower at 16.4 and 12.5%, respectively. Corresponding RCNRs were 6.8, 4.8, and 5.6, respectively.

### Histology

**Figure 7** shows myelin- and Nissl-stained sections of a rat 3 days after LPC injection in the dtg. We observed a demyelinated area ipsilaterally to the injection (**Figure 7C**) while the myelin content on the contralateral side appeared normal (**Figure 7B**).

On Nissl-stained sections, we observed an increased cell density attributed to gliosis, which overlapped with the demyelinating area (**Figure 7D**). On the contralateral side to the injection site, the amount and distribution of cells were normal (**Figure 7E**). Similarly as in the cc lesion, gliosis was diffusely distributed, and a high TABLE 1 | Pearson correlation of MRI parameters with optical density in myelin and Nissl stained sections including ipsilateral and contralateral ROIs and both LPC and vehicle injected animals.

Statistically significant R-values (p < 0.05, FDR adj. for testing 7 MRI parameters) are shown in bold.

Cc, corpus callosum; dtg, dorsal tegmental tract.

number of cells were found around blood vessels (**Figure 7F**).

## Correlation between MRI Parameters and Histology

RAFFn and MT parameters provided good correlation (R > 0.70) with myelin staining optical density when data from both brain areas and all ROIs were pooled (**Table 1**, **Figure 8**), while the correlation with myelin density was clearly lower for FA (R = 0.47) and for MD (R = 0.59). When only ROIs in the cc were included in the correlation analysis, all MRI parameters, except RD, markedly correlated (R > 0.66) with myelin density. However, the correlation with diffusion parameters was lower than that of RAFFn or MT despite larger relative changes in diffusion parameters than RAFFn and MT (**Table 1**). Importantly, in the dtg, RAFF4 and MT parameters had excellent correlation with myelin staining (R > 0.70), while there was no statistically significant correlation between diffusion parameters and myelin density.

There was a mild correlation between RAFF4 and Nissl staining intensity in the dtg (R < 0.59), however no statistically significant correlation in the cc and a weak correlation overall (R < 0.35). This indicates that cellularity is not the determining factor to contrast, but may have a small contribution for example in the case of profound gliosis (**Table 1**).

## Influence of pH in Protein Phantoms

Relaxation time maps were measured with different techniques from a BSA phantom with varying pH values. While T2 showed clear dependence on pH, due to altered proton exchange, as expected, RAFF4 was virtually unaffected by pH at values around a physiological pH of 7.4, in the range of 6.2–7.8 (**Figure 9**).

## DISCUSSION

In the present work, the value of the RAFF4 technique for detection of demyelination in the brain was investigated in a well-controlled LPC-injection rat model. Our major finding was that RAFF4 was able to robustly detect demyelination in both the cc and dtg, which have a different organization of myelinated axons and had a different amount of gliosis after LPC injection. Furthermore, RAFF4 provided excellent correlation with myelin content in both cases, while DTI was unable to reliably assess demyelinated area in the dtg, likely due to its sensitivity to underlying cytoarchitecture and gliosis. Interestingly, MT parameters measured with a recently introduced modified inversion technique (Mangia et al., 2011) were able to detect demyelination in both the cc and dtg and showed good correlation with myelin content, though the overall ability of MT at detecting demyelination was slightly lower than that of RAFF4.

Excellent correlation of RAFF4 to myelin content can be attributed to sensitivity of RAFF:n to correlation time regime in the ms-range (Satzer et al., 2015; Hakkarainen et al., 2016), which likely corresponds to exchange and dipolar interaction

of myelin water and dipolar interaction with carbon-bound methylene protons. In our previous study, there was no significant difference between RAFF4 and RAFF5 in their ability to detect myelin in normal brain (Hakkarainen et al., 2016) indicating that contribution of myelin to spectral density function is relatively flat in this correlation time range. As RAFF4 has shorter relaxation time constant than RAFF5, and therefore requires shorter magnetization preparation pulse train, RAFF4 was chosen over RAFF5 for sensitizing MRI for myelin. In our previous studies, RAFF1-3 showed sensitivity to shorter correlation time regime than RAFF4-5, with reduced gray/white matter contrast and reduced sensitivity to myelin (Liimatainen et al., 2010; Satzer et al., 2015). The different sensitivities of RAFFn techniques indicate how these techniques can be tuned to detect different motional regimes corresponding to different cellular components and pathological processes.

A conventional DTI approach was used as a comparative method, as it is one of the most commonly used advanced MRI approaches in clinical settings for assessment of white matter abnormalities. In the cc, DTI was able to detect LPC lesion even better than RAFF4 based on contrast between LPC and vehicle injected animals. However, large variability of the DTI metrics led to worse contrast-to-noise ratio between the two groups in comparison to RAFF4. In addition, the pattern of diffusion changes was unexpected: decreased MD, which was associated with decreased AD, while RD was unchanged. Demyelination is classically associated with increased RD (Song et al., 2005), as water molecules can pass axonal membranes more freely in the absence of myelin. In our study, demyelination, 3 days after LPC injection, was confirmed by quantitative analysis of myelin stained histological sections. Our results indicate substantial but not complete demyelination at this subacute time point. A previous electron microscopy study in LPC injected rats showed significant thinning of the myelin sheaths around axons, but not complete dissolution of myelin especially in the early time points (Foster et al., 1980). Our TEM data showed similar changes 3 days after LPC injection. In spite of clear demyelination, some disorganized myelin sheaths were still observed, which preserved RD close to normal. Furthermore, disorganized myelin sheaths with pockets, myelin debris, and vacuoles increased the number of diffusion restricting boundaries in the intra axonal space explaining the decreased AD, FA, and MD. These findings emphasize the complexity of using simple diffusion metrics as markers for a specific pathological feature such as demyelination. On the other hand, DTI, and especially more advanced high angular resolution diffusion-weighted imagingtype diffusion imaging approaches, clearly have potential to provide information that goes far beyond the myelin content and warrants further studies. The optimal MRI protocol to assess white matter pathology should contain both a technique with high sensitivity and specificity to myelin and capability to characterize microstructural changes.

Importantly, RAFF4 and RAFF5 have been shown to correlate with myelin density to a greater extent than MT in normal brain (Hakkarainen et al., 2016). Our results here show that using advanced MTR mapping (Mangia et al., 2011), MTR can also provide excellent correlation with myelin in LPC induced demyelination although lower than that of RAFF4. In addition to anisochronous and isochronous rotating frame relaxation pathways contributing to RAFF4 and 5, RAFF4 and 5 may share with MT cross-relaxation pathways. Therefore, these two techniques could provide complementary information for characterizing tissue integrity which could be utilized by proper modeling (Underhill et al., 2011).

A benefit of the LPC model is that it causes demyelination without significant edema formation and only mild to moderate gliosis. Therefore the LPC model allows assessment of demyelination without multiple simultaneous contributing factors. In practically all human pathological conditions, demyelination is co-localized with a number of other pathological processes including inflammation, edema and acidification. There was more gliosis in the dtg than in cc after LPC injections, which likely explains why there was statistically significant correlation with Nissl staining in dtg but not in cc. However, in both cases correlation with myelin content remained good, and in general RAFF4 correlation with myelin density was higher than with other MRI methods studied. Importantly, the correlation of RAFF4 with Nissl stain proves that mechanisms other than myelin must contribute to the RAFF4 contrast. Therefore, RAFF4 is not entirely specific to myelin, as expected for any relaxation parameters that unavoidably depends on multiple biological processes. Yet, the insensitivity of RAFF4 to at least some processes such as pH changes, along with the distinct feature of RAFF4 to detect slow and ultra-slow motional regimes typical in myelin, overall suggest also a better specificity of RAFF4 to myelin as compared to other methods.

A unique feature of RAFF4 is that sensitivity to slow molecular motion can be achieved within the SAR limits of human studies, making the RAFF4 technique readily applicable to clinical settings. Indeed, excellent gray/white matter contrast with high rank RAFF4 has already been demonstrated in human brain with reasonable acquisition times (∼10 min) (Liimatainen et al., 2015). While sensitivity to ms-time range can be also achieved with conventional spin-lock techniques, such as continues-wave T1ρ with spin-lock field in kHz range (Sepponen et al., 1985; Gröhn et al., 2000), human studies are compromised because of high SAR of on-resonance spin-lock pulses.

## CONCLUSION

The excellent ability of RAFF4 to detect myelin content in healthy and pathological tissues and its insensitivity to fiber orientation distribution, gliosis and pH together with low SAR, promises RAFFn to become a useful technique for visualizing demyelinating lesions. Future work will aim to investigate later time points after LPC-induced demyelination in order to follow the remyelination of the LPC injected animals and to test specificity of RAFF4 for myelin in other disease models with more complex pathology which for instance include edema.

## REFERENCES

## AUTHOR CONTRIBUTIONS

LL, AA, and AS participated in design of the work, acquisition, analysis, interpretation of data and preparing the manuscript. LT participated in design of the work, acquisition and preparing the manuscript. LE participated in design of the work, analysis and preparing the manuscript. SMa participated in design of the work, interpretation of the data and preparing the manuscript. AN participated in design of the work and preparing the manuscript. SMi and OG participated in design of the work, analysis, interpretation of data and preparing the manuscript.

## FUNDING

The work was supported by EU H2020 Marie Skłodowska RISE project #691110 (MICROBRADAM), Academy of Finland grants #275453 and #298007, EU FP7-HEALTH project #602102 (EPITARGET), NIH BTRC grants, P41 EB015894, P30 NS057091. OG was supported by Professor Pool grant from Svenska Kulturfonden, SMi was supported by Fulbright-Saastamoinen Foundation Grant in Health and Environmental Science. LL was supported by MnDRIVE and the Emil Aaltonen Foundation.

## ACKNOWLEDGMENTS

Electron microscopy was performed utilizing Sib-labs infrastructure at the University of Eastern Finland.

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

Copyright © 2017 Lehto, Albors, Sierra, Tolppanen, Eberly, Mangia, Nurmi, Michaeli and Gröhn. 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.

# Determining Excitatory and Inhibitory Neuronal Activity from Multimodal fMRI Data Using a Generative Hemodynamic Model

#### Martin Havlicek\*, Dimo Ivanov, Alard Roebroeck and Kamil Uludag˘

Department of Cognitive Neuroscience, Faculty of Psychology and Neuroscience, Maastricht University, Maastricht, Netherlands

#### Edited by:

Itamar Ronen, Leiden University, Netherlands

#### Reviewed by:

Nelson Jesús Trujillo-Barreto, University of Manchester, United Kingdom Fahmeed Hyder, Yale University, United States

\*Correspondence:

Martin Havlicek<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 16 May 2017 Accepted: 23 October 2017 Published: 10 November 2017

#### Citation:

Havlicek M, Ivanov D, Roebroeck A and Uludag K (2017) Determining ˘ Excitatory and Inhibitory Neuronal Activity from Multimodal fMRI Data Using a Generative Hemodynamic Model. Front. Neurosci. 11:616. doi: 10.3389/fnins.2017.00616 Hemodynamic responses, in general, and the blood oxygenation level-dependent (BOLD) fMRI signal, in particular, provide an indirect measure of neuronal activity. There is strong evidence that the BOLD response correlates well with post-synaptic changes, induced by changes in the excitatory and inhibitory (E-I) balance between active neuronal populations. Typical BOLD responses exhibit transients, such as the early-overshoot and post-stimulus undershoot, that can be linked to transients in neuronal activity, but they can also result from vascular uncoupling between cerebral blood flow (CBF) and venous cerebral blood volume (venous CBV). Recently, we have proposed a novel generative hemodynamic model of the BOLD signal within the dynamic causal modeling framework, inspired by physiological observations, called P-DCM (Havlicek et al., 2015). We demonstrated the generative model's ability to more accurately model commonly observed neuronal and vascular transients in single regions but also effective connectivity between multiple brain areas (Havlicek et al., 2017b). In this paper, we additionally demonstrate the versatility of the generative model to jointly explain dynamic relationships between neuronal and hemodynamic physiological variables underlying the BOLD signal using multi-modal data. For this purpose, we utilized three distinct data-sets of experimentally induced responses in the primary visual areas measured in human, cat, and monkey brain, respectively: (1) CBF and BOLD responses; (2) CBF, total CBV, and BOLD responses (Jin and Kim, 2008); and (3) positive and negative neuronal and BOLD responses (Shmuel et al., 2006). By fitting the generative model to the three multi-modal experimental data-sets, we showed that the presence or absence of dynamic features in the BOLD signal is not an unambiguous indication of presence or absence of those features on the neuronal level. Nevertheless, the generative model that takes into account the dynamics of the physiological mechanisms underlying the BOLD response allowed dissociating neuronal from vascular transients and deducing excitatory and inhibitory neuronal activity time-courses from BOLD data alone and from multi-modal data.

Keywords: excitatory-inhibitory, multi-modal data, fMRI signal modeling, response transients, neuronal adaptation, post-stimulus BOLD undershoot, hemodynamic uncoupling, DCM

## INTRODUCTION

Functional magnetic resonance imaging (fMRI) is a widely used non-invasive technique to assess brain function. The fMRI signal reflects neuronal activity only indirectly through the measurements of accompanying hemodynamic processes at temporal resolution typically on the order of seconds and spatial resolution typically on the order of tens of cubic millimeters. In general, neuronal activation causes a series of physiological events, including localized changes in cerebral blood flow (CBF), cerebral metabolic rate of oxygen (CMRO2), cerebral blood volume (CBV), and deoxyhemoglobin content. These physiological variables form the basis of the blood oxygenation level-dependent (BOLD) signal (Ogawa et al., 1990), the most commonly used fMRI approach. However, CBF and CBV can also be directly measured with various fMRI techniques (e.g., see Wong et al., 1997; Lu et al., 2003; Liu and Brown, 2007; Huber et al., 2014b).

In many studies, it has been found that there is high correspondence between response properties measured in fMRI and other hemodynamic techniques and those measured from invasive electrical recordings, mostly acquired in non-human primates, cats and rodents (Logothetis et al., 2001, 2010; Kim et al., 2004; Niessing et al., 2005; Devor et al., 2007, 2013; Logothetis, 2008; Muckli, 2010; Boynton, 2011; Hillman, 2014). In particular, CBF and BOLD responses show better correlation with post-synaptic local field potentials (LFPs) than with spiking activity (multi-unit activity, MUA), suggesting that the hemodynamic response reflects stronger the input to a neuronal population in a brain area and intrinsic processing (Lauritzen, 2005) rather than the output of that area (Goense and Logothetis, 2008; Logothetis et al., 2010; Magri et al., 2012). Positive CBF and BOLD responses during stimulation are associated with an increase in neuronal activity and decrease in deoxyhemoglobin content, whereas negative CBF and BOLD responses are associated with a decrease in neuronal activity below baseline and increase in deoxyhemoglobin content (e.g., Shmuel et al., 2006 and references therein).

Typical positive BOLD responses to sustained stimulation display transients, such as response adaptation (also referred to as early-overshoot) and post-stimulus undershoot (Frahm et al., 1996; Krüger et al., 1996; Hoge et al., 1999). Similarly, neuronal responses to stimulation exhibit rapid rise followed by a decay (or adaptation) to a steady-state level and are often followed by a brief decrease below baseline after stimulus cessation (e.g., Logothetis et al., 2001). These neuronal transients (Boynton et al., 1996; Hoge et al., 1999; Bandettini and Ungerleider, 2001; Birn et al., 2001) are the result of changes in excitatory and inhibitory (E-I) balance between active neuronal populations, controlled by local micro-circuitry but also by long-range connections with other brain areas (Logothetis, 2002; Logothetis and Wandell, 2004; Shmuel et al., 2006; Hyder et al., 2010; Havlicek et al., 2017b). CBF reflects these neuronal transients in a temporally smoothed fashion (Hoge et al., 1999; Attwell and Iadecola, 2002; Uludag˘ et al., 2004; Sadaghiani et al., 2009; Attwell et al., 2010; Cauli and Hamel, 2010; Mayhew et al., 2014). Because of the complex underlying physiological processes, the BOLD response can exhibit transients not only from neuronal sources, but also due to the properties of blood vessels: the BOLD response is dominated by signal contributions originating from the venous compartments, and venous CBV can be dynamically uncoupled from CBF (i.e., venous CBV lags behind CBF), influencing the amplitude of the early-overshoot and post-stimulus undershoot of the BOLD-response (Buxton et al., 1998b; Mandeville et al., 1999; Yacoub et al., 2006; Chen and Pike, 2009). Alternatively, dynamic uncoupling between CBF and CMRO2 (i.e., CMRO2 lags behind CBF) could result in the same BOLD transients as well (Lu et al., 2004; Frahm et al., 2008; Donahue et al., 2009; Hua et al., 2011; Poser et al., 2011; van Zijl et al., 2012). However, we have recently argued, supported by modeling of experimental data, that the contribution of CBF-CMRO2 uncoupling to BOLD signal transients is much lower than that of CBF-venous CBV uncoupling (Havlicek et al., 2017a). Thus, neuronal and hemodynamic responses in different areas (or voxels) and subjects exhibit dynamic features, which can be both due to changes in E-I balance or due to biomechanical properties of the vasculature.

In standard analysis of fMRI data, linear convolution is applied between a stick or box-car functions (representing the stimulation paradigm) and assumed canonical hemodynamic response function (representing the combined transform from stimulus time-course to neuronal signal and the measured BOLD response) (Friston et al., 1995). However, as indicated above, excitatory and inhibitory neuronal responses may be nonlinearly related to stimulation (e.g., see Boynton et al., 1996; Bandettini and Ungerleider, 2001; Birn et al., 2001; Grill-Spector and Malach, 2001; Kida and Yamamoto, 2008; Mullinger et al., 2013, 2014; Pérez-González and Malmierca, 2014; Havlicek et al., 2017a; Keller et al., 2017). Furthermore, vascular transients resulting from dynamic uncoupling induce additional nonlinearities between the input function and subsequent hemodynamic variables. Therefore, the linearity assumption in BOLD data analysis may not be sufficient in many experiments and, consequently, the inferred information about the neuronal activity changes obtained from hemodynamic signals using linear (de)convolution analysis might be confounded with vascular effects. In other words, to more accurately estimate neuronal responses from the BOLD signal, a nonlinear model that accounts for dynamic relationships between neuronal and hemodynamic physiological variables underlying the BOLD response is needed.

Recently, we have introduced a physiologically-informed generative model of the BOLD signal within the framework of dynamic causal modeling (Friston et al., 2003) (called P-DCM)1 , linking excitatory and inhibitory neuronal activity to the BOLD response (Havlicek et al., 2015). In P-DCM, we employ: (i) an adaptive excitation-inhibition neuronal model that accounts for a wide range of neuronal time-course both during stimulation and post-stimulation periods; (ii) a neurovascular coupling (NVC) model that links neuronal activity to blood flow in a strictly

1As we have developed our novel generative model within the framework of the dynamic causal model (DCM) and applied it both to single ROI and to a network of ROIs, in this paper to avoid confusion, we use DCM and the generative model embedded within P-DCM interchangeably.

feedforward fashion; (iii) a balloon model (Buxton et al., 1998b) that can account for a vascular uncoupling between CBF and venous CBV; and (iv) field strength and sequence dependent parameterization of the BOLD signal equation. We compared P-DCM with other DCM models (Friston et al., 2003; Marreiros et al., 2008) and demonstrated significant improvements in the ability to model commonly observed neuronal and vascular response transients in single regions (Havlicek et al., 2015) and also within a network of several regions with task-driven activity changes (Havlicek et al., 2017b). In the latter case, we also showed a high fidelity of P-DCM to jointly explain CBF and BOLD responses simultaneously measured with the arterial spin labeling (ASL) fMRI technique, demonstrating the benefit of additional information provided by a CBF measurement for model inversion.

In general, multi-modal imaging is a powerful approach to study the relationship between neural activity and the BOLD fMRI signal. Measurements of different physiological variables can increase the ability to disambiguate neuronal and vascular effects present in the BOLD signal and potentially unravel limitations of the hemodynamic models. In the current work, we aim to explore the versatility of P-DCM to explain dynamic relationships between various combinations of measured physiological variables and to deduce the excitatory and inhibitory neuronal dynamics from hemodynamic data. This is done under the constraints of assumed physiological mechanisms and experimental manipulations. In particular, we use: (1) newly acquired CBF and BOLD responses to static and flickering stimuli in human subjects; (2) CBF, total CBV and BOLD response to square-wave grating stimulus acquired in the cat brain from the study of Jin and Kim (2008); and (3) positive and negative neuronal and BOLD responses induced by rotating visual stimuli measured in the monkey brain from the study of Shmuel et al. (2006). In our modeling, we emphasize stimulus-type-dependent modulation of response transients that can be linked to a dynamic interplay between excitatory and inhibitory activity. In addition, we allow for differences between stimulation and post-stimulation response periods and account for vascular-, magnetic field strength-, and MRI sequence-dependent properties. The current approach can also be generalized to other invasive and non-invasive multimodal data, such as EEG-fMRI, provided generative models exist for both modalities.

## METHODS

## General Description of P-DCM

The generative model in P-DCM consists of four causally-linked components that define how the neuronal signal is transformed to the measured BOLD response (see **Figure 1** for an illustration and the summary of model equations). For a more detailed description of the model and its comparison with previous DCM models, please see (Havlicek et al., 2015, 2017b).

## Neuronal Model of E-I Balance

In this model, an exogenous input u (t) (e.g., sensory stimulus) drives the change in excitatory activity, ne (t), which is directly coupled with a change in inhibitory activity, ni (t). The strength of this input, expressed with the parameter c, scales the amplitude of the neuronal response. The shape of the neuronal response

is tuned by the transient imbalance between excitatory and inhibitory activity. In particular, the typical overshoot at the onset can be produced by gradually increasing inhibitory activity that modulates the excitatory activity via negative feedback. Next, persistence of the inhibitory activity following stimulus cessation can produce post-stimulus deactivation. This temporal evolution of excitatory and inhibitory activity, including their dynamic mismatch, is controlled by the parameters σ and λ, respectively, and the strength of the inhibitory activity modulating the excitatory activity is encoded by the parameter µ. Optimization of these neuronal parameters allows modeling a broad repertoire of neuronal response adaptation profiles and of possible poststimulus deactivations separately for the stimulation period (SP) and the post-stimulation period (PSP)2 .

#### Neurovascular Coupling (NVC):

The output of the neuronal model, i.e., the excitatory activity modulated by inhibitory activity, is transformed to CBF, f(t), in a strictly feedforward fashion, via vasoactive signal a(t). Neuronal excitation/inhibition leads to arterial vasodilatation/vasoconstriction associated with increased/decreased CBF (Devor et al., 2007). Thus, the modeled neuronal response transients are conveyed to a CBF response, albeit in a smooth version. Decay and delay of the CBF response with respect to the neuronal response is regulated with three constants, ϕ, φ and χ, with only χ being optimized during model inversion.

## Hemodynamic Model:

The hemodynamic model is represented by the balloon model (Buxton et al., 1998a,b). It models the mass balance of normalized changes in CBV, v(t), and deoxyhemoglobin content, q(t), as they pass through the venous compartment. Their changes are driven by changes in the inflowing CBF, f(t), and CMRO2, m(t), respectively. It is assumed that CBF and oxygen extraction faction, E(t), are dynamically coupled, thus m (t) = f(t) · E(t) E0 , where E0 is the value of oxygen extraction faction at rest (please see discussion of this assumption in (Havlicek et al., 2015), and the Discussion section for more details). Furthermore, during steady-state, the blood leaving the venous compartment (i.e., the outflow, fout(t)) and the venous CBV are coupled via a power law relationship (Grubb et al., 1974), with exponent α, whereas during the transient periods venous CBV and CBF can be uncoupled; i.e., venous CBV can evolve more slowly than CBF. This is due to the vessel's resistance to changes in venous CBV, described by the viscoelastic time constant, τ . Theoretically, the viscoelastic time constant can have different values during SP and PSP, τSP and τPSP. The time dimension of changes in v(t) and q(t) is scaled by the mean transit time of the blood through the venous compartment at rest, t0. Note that t0 is linked to the resting venous blood volume fraction V0 via the central volume principle, t0 = V0 F0 , where F0 is the blood flow at rest.

## BOLD Signal Equation:

The BOLD signal reflects changes in the deoxyhemoglobin content, q(t), together with changes in the deoxyhemoglobin concentration, q(t)/v(t), and venous CBV, v(t). Their relative contribution is weighted by parameters that are magnetic field strength-, TE- and MRI sequence-dependent (Uludag et al., 2009; ˘ Havlicek et al., 2015, 2017a).

In summary, as illustrated in **Figure 1**, P-DCM and its parameters allows the BOLD response to exhibit transients, such as response adaptation during stimulation and post-stimulus undershoot that can have both neuronal and vascular origins. As we show below, physiological origins of these transients can be tested under the constraints of concurrent multi-modal physiological data and experimental manipulations.

## Data Description

To demonstrate the utility and versatility of P-DCM, below we describe three different data-sets acquired from human, cat, and monkey brains. Each data-set consists of a different combination of physiological measurements: I. CBF and BOLD response; II. CBF, total volume CBV, and BOLD responses, published in Jin and Kim (2008); and III. neuronal and BOLD responses, published in Shmuel et al. (2006). Any additional physiological data, next to the BOLD data (and/or experimental manipulation), provide physiologically-informed constraints on the underlying mechanisms of the BOLD response. This can result in a more accurate inference on the changes in E-I balance forming the neuronal responses (Havlicek et al., 2017b) and potentially inform about the limitations of model structure and parameter assumptions.

### I. CBF and BOLD Responses