Multi-echo EPI of human fear conditioning reveals improved BOLD detection in ventromedial prefrontal cortex

Author’s Accepted Manuscript Multi-echo EPI of Human Fear Conditioning Reveals Improved BOLD Detection in Ventromedial Prefrontal Cortex Brice Fernandez, Laura Leuchs, Philipp G. Sämann, Michael Czisch, Victor I. Spoormaker www.elsevier.com

PII: DOI: Reference:

S1053-8119(17)30395-6 http://dx.doi.org/10.1016/j.neuroimage.2017.05.005 YNIMG14017

To appear in: NeuroImage Received date: 16 November 2016 Revised date: 24 April 2017 Accepted date: 4 May 2017 Cite this article as: Brice Fernandez, Laura Leuchs, Philipp G. Sämann, Michael Czisch and Victor I. Spoormaker, Multi-echo EPI of Human Fear Conditioning Reveals Improved BOLD Detection in Ventromedial Prefrontal Cortex, NeuroImage, http://dx.doi.org/10.1016/j.neuroimage.2017.05.005 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Multi-echo EPI of Human Fear Conditioning Reveals Improved BOLD Detection in Ventromedial Prefrontal Cortex a12*

b3

b2

b2

b2

Brice Fernandez , Laura Leuchs , Philipp G. Sämann , Michael Czisch , Victor I. Spoormaker a Applications & Workflow, GE Healthcare, Munich, Germany b Neuroimaging Unit, Max Planck Institute of Psychiatry, Munich, Germany *Correspondence to: GE Healthcare, Service Hospitalier Frédéric Joliot – CEA, 4 place du Général Leclerc, 91401 Orsay Cedex, France. [email protected]

ABSTRACT Standard

weighted functional magnetic resonance imaging (fMRI) performed with echo-planar imaging

(EPI) suffers from signal loss in the ventromedial prefrontal cortex (vmPFC) due to macroscopic field inhomogeneity. However, this region is of special interest to affective neuroscience and psychiatry. The Multi-echo EPI (MEPI) approach has several advantages over EPI but its performance against EPI in the vmPFC has not yet been examined in a study with sufficient statistical power using a task specifically eliciting activity in this region. We used a fear conditioning task with MEPI to compare the performance of MEPI and EPI in vmPFC and control regions in 32 healthy young subjects. We analyzed activity associated with short (12 ms), standard (29 ms) and long (46 ms) echo times, and a voxel-wise combination of these three echo times. Behavioral data revealed successful differentiation of the conditioned versus safety stimulus; activity in the vmPFC was shown by the contrast “safety stimulus > conditioned stimulus” as in previous research and proved significantly stronger with the combined MEPI than standard single-echo EPI. Then, we aimed to demonstrate that the additional cluster extent (ventral extension) detected in the vmPFC with MEPI reflects activation in a relevant cluster (i.e., not just nonneuronal noise). To do this, we used resting state data from the same subjects to show that the timecourse of this region was both connected to bilateral amygdala and the default mode network. Overall, we demonstrate that MEPI (by means of the weighted sum combination approach) outperforms standard EPI in vmPFC; MEPI performs always at least as good as the best echo time for a given brain region but

1

Applications & Workflow, GE Healthcare, Oskar-Schlemmer-Str. 11, 80807 Munich, Germany Current address: GE Healthcare, Service Hospitalier Frédéric Joliot – CEA, 4 place du Général Leclerc, 91401 Orsay Cedex, France 3 Max Planck Institute of Psychiatry, Kraepelinstr. 2-10, 80804 Munich, Germany. 2

1

provides all necessary echo times for an optimal BOLD sensitivity for the whole brain. This is relevant for affective neuroscience and psychiatry given the critical role of the vmPFC in emotion regulation. ABBREVIATION BOLD, Blood-Oxygen-Level Dependent; dACC, dorsal Anterior Cingulate Cortex; DMN, Default Mode Network; EPI, Echo Planar Imaging; FWE, Family-wise Error; MEPI, Multi-echo EPI; PET, Positron Emission Tomography; PTSD, Post-Traumatic Stress Disorder; MDD, Major Depressive Disorder; mPFC, medial Prefrontal Cortex; vmPFC, ventromedial Prefrontal Cortex

KEYWORDS Echo Planar Imaging (EPI), Multi-echo EPI (MEPI), Fear conditioning, Ventromedial Prefrontal Cortex, Default Mode Network (DMN), Functional Connectivity Orientation: S superior I

Inferior

A

anterior

P

posterior

L

Left

R

Right

INTRODUCTION In the last two decades, functional magnetic resonance imaging (fMRI) has become one of the most popular imaging tools to study brain function in a variety of applications like surgical planning (Vlieger et al., 2004), basic neuroscience (Poldrack, 2012) and clinical neuropsychiatric research (Zhan & Yu, 2015; Mitterschiffthaler et al., 2006; Etkin & Wager, 2007; Groenewold et al., 2013). Gradient-echo echo-planar imaging (EPI) is the most commonly used and accepted pulse sequence for fMRI because of its relative robustness to motion, its high temporal resolution and its good sensitivity to blood-oxygen-level dependent (BOLD) effect (Mansfield, 1977; Ogawa et al., 1990). However, there are considerable issues with this pulse sequence mainly due to the susceptibility-induced field gradients (Koch et al., 2009; Farzaneh et al., 2

1990), which lead to severe signal and BOLD sensitivity loss (Deichmann et al., 2002) especially pronounced in the ventromedial prefrontal cortex (vmPFC) at 3.0 T. In particular, Deichmann et al. (2002) showed that in susceptibility affected brain regions such as the most ventral parts of the vmPFC, the local echo time (i.e. local to a given voxel or brain region) is shifted to lower or higher values, possibly leading to a signal void when the echo occurs outside of the EPI readout window. This implies that a whole range of echo times (

) are needed to properly sample the BOLD response in the vmPFC. Additionally, still in

the vmPFC, the more ventral the considered voxel is, the more shifted is the local echo time, meaning that the vmPFC is not uniformly affected by susceptibility induced gradient fields. In the most ventral parts of the vmPFC, the BOLD signal is expected to be better sampled by a short echo time ( the upper part of the vmPFC it is expected to be better sampled by a longer

ms), while in ms).

Unfortunately, the vmPFC is strongly affected by susceptibility-induced field gradients while being a key region for many affective processes. The medial part of this region has been repeatedly associated with fear conditioning and extinction (Milad & Quirk, 2012; Groenewold et al., 2013; Hayes et al., 2012; Milad & Rauch, 2007). Fear extinction is impaired in anxiety disorders and post-traumatic stress disorder (PTSD) (Duits et al., 2015) and hypoactivity in vmPFC during fear extinction learning and recall has been observed in patients with PTSD as well (Milad et al., 2009). The most ventral parts of the mPFC (overlapping with the medial orbitofrontal cortex) have direct projections to amygdala, ventral striatum, lateral hypothalamus, and hippocampus (Milad & Rauch, 2007), which are also involved in affective processing and stress response integration, making the vmPFC ideally posed as a mediator of circuitry subserving in emotional responses. The role of the vmPFC can extend to more cognitive processes as well, as demonstrated by one study showing that cognitive emotion regulation (i.e., subjects trying to modulate their own emotional responses to stimuli) is also mediated by the vmPFC (Delgado et al., 2008). Interestingly, studies on fear learning using standard EPI generally do not show involvement of the most ventral parts of mPFC in fear learning (Fullana et al., 2015), whereas an anatomical study (that uses acquisitions less prone to signal loss and distortion) focusing on thickness of vmPFC in relation to extinction learning reported the strongest correlation in the most ventral part of the mPFC (Milad et al., 2005).

3

Furthermore, a general emotion and stress regulation function of the vmPFC might also help understand why this region, at least the more subgenual parts, showed differences in functional connectivity during the resting-state in major depressive disorder (MDD) patients (Greicius et al., 2007). Particularly its contribution to the default mode network (DMN) was heightened in patients, which was found generally plausible as the DMN is associated with an internal mode of information processing including episodic memory processing (Buckner et al., 2008). Additionally, positron emission tomography (PET) imaging has highlighted a role of the subgenual mPFC in MDD (Drevets et al., 1997). There are over 2,000 PubMed listed fMRI publications on the DMN in the last decade, and interestingly, PET studies reported integration of the most ventral parts of mPFC in the DMN (Buckner et al., 2008) whereas this is less clear for fMRI studies (Andrews-Hanna et al., 2014), even though the medial but not the lateral parts of the ventral parts of the PFC have been found to correlate with the DMN (Zald et al., 2014). We propose that inconsistencies in the literature on the involvement of the vmPFC in task activity or resting state functional connectivity may be due to technical limitations of standard EPI pulse sequences even though other potential influences (e.g., study design, statistical analysis choices, sample sizes) cannot be excluded. In the literature, many solutions have been proposed to mitigate signal dropouts in fMRI. These can be classified in three categories: (i) shimming techniques including z-shimming (Constable & Spencer, 1999; Deichmann et al., 2002; Weiskopf et al., 2006), 2D/3D x/y/z-shimming (Glover, 1999; Weiskopf et al., 2007), passive shimming (Cusack et al., 2005; Koch et al., 2006) and active shimming (Hsu & Glover, 2005; Juchem et al., 2014); (ii) tailored RF pulses (Cho & Ro, 1992; Stenger et al., 2000; Wastling & Barker, 2015; Yip et al., 2006) and through-plane phase precompensated RF pulses (Chen & Wyrwicz, 1999; Yang et al., 2010; Yang et al., 2011; Yang et al., 2012; Yip et al., 2009) and; (iii) other pulse sequence and reconstruction techniques including parallel imaging (Preibisch et al., 2003; Preibisch et al., 2008), improved k-space sampling strategies like spiral in/out (Glover & Law, 2001) combined gradient and spin-echo EPI (Schwarzbauer et al., 2010) and Multi-echo EPI (MEPI) (Poser et al., 2006; Posse, 2012; Posse et al., 1999) which can easily take advantage of parallel imaging (Schmiedeskamp et al., 2010). Active and passive shimming techniques usually require additional hardware that is not always compatible with a full fMRI setup (e.g. 32 channel coil, video screen, goggles, electroencephalography, 4

electrooculography, eye-tracker and so on) and can lead to subject discomfort. Slice shimming techniques have the disadvantage of significantly lengthening the repetition time (

). Tailored and phase

precompensated RF pulses have unclear effects on BOLD sensitivity (Yip et al., 2009) and temporal signal-to-noise ratio (Wastling & Barker, 2015). The combined gradient and spin-echo EPI approach (Schwarzbauer et al., 2010) might possibly be difficult to use in resting state as it implies mixing two different types of signals with different origin and sensitivity which might not be correlated (i.e. gradientecho and spin-echo). The use of parallel imaging such as SENSE and GRAPPA (Griswold et al., 2002; Pruessmann et al., 1999) is now well established in fMRI (Lütcke et al., 2006; Schmidt et al., 2005) and fits particularly well with MEPI as previously reported (Schmiedeskamp et al., 2010). Moreover, recent advances in analysis strategies based on MEPI showed the possibility of separating BOLD from nonBOLD components and improving specificity (Kundu et al., 2012; Kundu et al., 2013; Kundu et al., 2015). It is worth noting that a

of 30 ms is typically used at 3T which is lower than the optimum for gray matter

but represents a good compromise between BOLD sensitivity, image quality and temporal resolution (Norris, 2006). All these previous results on MEPI and on the BOLD sensitivity gave additional motivation to further develop MEPI and its associated analysis pipeline (Posse et al., 1999). Another motivation for MEPI is that

varies across the brain (Peters et al., 2007) and consequently, the optimal

varies

accordingly (Gati et al., 1997). In the present study, we are seeking a generic methodology that allows for whole brain acquisition with reasonable spatial and temporal resolution (i.e. whole brain coverage at 3-4 mm isotropic spatial resolution in 2-2.5 s acquisition time per volume) and a better BOLD sensitivity in the vmPFC without compromising the BOLD sensitivity in B0-homogeneous brain regions. Since the methodology is intended to be used routinely in fMRI, it should be not more complicated to use than the standard EPI. Consequently, we chose to develop and evaluate MEPI associated with parallel imaging and the weighted sum combination approach (Poser et al., 2006). This approach has the advantages of including the standard EPI and of being compatible with the latest advanced development like multiband EPI with blipped CAIPI (Setsompop et al., 2012). Several papers showed that MEPI had a tendency to perform better than the standard EPI approach (Poser et al., 2006; Posse et al., 1999; Weiskopf et al., 2005; Posse, 2012; Schmiedeskamp et al., 2010), and a recent study indicated an enhanced sensitivity of MEPI 5

in the orbitofrontal cortex using an emotional learning and a reward based learning task (Kirilina et al., 2016). However, according to several meta-analyses, the emotional learning and reward based learning tasks employed do not seem to specifically elicit activity in the ventral part of the mPFC (Garrison et al., 2013; Hayes et al., 2014; Sabatinelli et al., 2011; Zhang et al., 2013). Consequently, what is needed to remove uncertainty about the superior performance of MEPI in vmPFC is a whole brain, statistically robust comparison between MEPI data and standard EPI data with appropriate statistical power using a task specifically eliciting activity in vmPFC. We aimed to do this by using a fear conditioning task which is known to elicit activity in the vmPFC, the target region, the dorsal anterior cingulate cortex (dACC) and the bilateral insula as control regions (Fullana et al., 2015). We hypothesized that MEPI with the weighted sum combination approach would yield a statistically stronger effect in the vmPFC for the contrast safety stimulus > conditioned stimulus (

) than standard EPI, particularly in the most ventral parts. As a

second hypothesis, we expect the anticipated differences in vmPFC to be anatomically and neurophysiologically relevant, meaning that the tentative vmPFC cluster should be connected to the bilateral amygdala in a separate data-set of resting-state data, among other regions. THEORY This section will introduce the notation, convention and definition which will be used throughout the paper. The pulse sequence developed in this work is a single shot multi-echo echo-planar-imaging (MEPI) similar to the one used by Poser et al. (2006) and is schematically described in Figure 1. Here, each EPI read-out has a given echo-time

and leads to one image named here an “echo”. Consecutive echoes are

separated by a tunable parameter

. For simplicity, the different echoes are named

the echo number. The signal intensity of

is governed by the

) with

with

indicating

decay:

)

. Deichmann et al. (2002) derived the BOLD sensitivity (

) of

as:

, and in the same paper demonstrated theoretically and experimentally the effects of through-plane and inplane susceptibility induced inhomogeneity. It is shown that, in susceptibility affected brain regions, the local echo time (i.e. local to a given voxel or brain region) is shifted to lower or higher values, possibly 6

leading to a signal void when the echo occurs outside of the EPI readout window. It is worth noting that a of 30 ms is usually used at 3T which is somewhat lower than the optimum for gray matter but represents a good compromise between BOLD sensitivity, image quality and temporal resolution (Norris, 2006). These results gave additional motivation to further develop MEPI and its associated analysis pipeline as suggested earlier (Posse et al., 1999). Another motivation for MEPI is that brain (Peters et al., 2007) and consequently, the optimal

varies across the

varies accordingly (Gati et al., 1997). Hence,

MEPI should allow a wide range of local echo times to be optimally sampled in the presence of brain induced susceptibility gradients. To analyze these data in an optimal way, one approach is to create a composite dataset (

) using a weighted sum of the different echoes (Poser et al., 2006) on a voxel by

voxel basis as: ∑

with

(1)

the number of echoes and ∑

The weights can be computed based on their relative expected BOLD sensitivity as proposed by Poser et al. (2006): ̂)

( ∑

(

with ̂ being an estimate of

̂)

[2]

for a given voxel. ̂ can be computed directly from the MEPI data using

nonlinear fitting.

7

Figure 1: Description of the pulse sequence used in this work. Additional echo-navigators were acquired prior to each EPI read-out but only the data acquired prior to the first echo were used to update the static correction maps over time (see text for details).

METHODS Subjects In this study, 36 healthy subjects (17 females), aged 25.0 ± 3.2 years (mean ± standard deviation), were enrolled after being checked for standard MRI exclusion criteria (pacemaker, cochlear implants, insulinpumps, in-body metal splinters, claustrophobia etc.) and general medical, neurological and psychiatric history. A short structural MRI screening protocol was performed to exclude brain anomalies or normal variants incompatible with spatial fMRI processing. All subjects provided written informed consent after a full explanation of the protocol and were reimbursed for their participation. The study protocol was in accordance with the last version of the Helsinki Declaration and approved by the local ethics committee (Medical Faculty of the Ludwig-Maximilian University, Munich, Germany). An a priori power analysis indicated that an effective sample of 27 subjects was needed to achieve 80% power ( detecting a medium effect (Cohen’s

) for

) in a one-tailed paired samples t-test as computed by

G*Power 3.1 (Faul et al., 2007) and took typical

dropout rates into account.

Paradigm In the fear conditioning task (see Figure 2), three visual stimuli (three squares in different colors) were presented 10 times each in a pseudorandomized order using Presentation Software (Neurobehavioral 8

Systems Inc., Berkeley, CA, USA) with inter-trial-intervals varying between 12-14 s. Two of these stimuli were paired with mild electric shocks, which were administered to the right wrist through MR compatible gold electrodes (Digitimer Stimulator, Digitimer Ltd., Hertfordshire, UK). Shock intensity was individually titrated to a level that was reported by the subjects to be uncomfortable but not painful. One of the conditioned stimuli was always paired with a shock (100% reinforcement – labelled

), one was

paired with shocks in 60% of the trials (pseudo-randomized, 60% reinforcement – labelled two stimuli are summarized as labelled

). Stimulus duration of

). These

. The third stimulus was never paired with a shock (safety stimulus – and

was 4 s, the shock to the

occurred at 3.5 s and

lasted for 2 ms.

Figure 2: Description of the fear conditioning task. The three visual stimuli (colored squares) were either never paired with a shock ( summarized as

), or paired with a shock in 60% (

) or 100% of the trials (

).

and

are

. The three basic conditioned stimuli were presented 10 times each in a pseudorandom order as

depicted above. Three behavioral ratings were performed before, in the middle of and after the fear conditioning task.

Before the fear conditioning session started, subjects underwent a brief habituation session in which the three geometric forms were presented two times each. After this habituation session, the fear conditioning paradigm began with a behavioral rating, in which the three stimuli were presented successively and subjects had to indicate the shock expectancy per stimulus on a scale from 0% to 100% (ten steps in between, ten seconds duration per rating). These behavioral ratings were repeated in the middle of the 9

fear conditioning task (after each stimulus had been presented 5 times) and after the fear conditioning task, which served as an index for whether subjects learned the stimulus contingencies. In total, the fear conditioning paradigm lasted around 10 minutes (236 volumes). The statistical contrast

was expected to show increased activity in the vmPFC (Fullana et al.,

2015), which is the main region of interest for our analysis as it is strongly affected by susceptibility gradients. The reverse contrast

commonly shows an increase of activity in dACC and bilateral

insula (Fullana et al., 2015). These regions are known to be unaffected by signal loss (i.e.

-homogenous

brain regions) and were used as control regions. MR acquisition All data were acquired on a 3.0T MR scanner (MR750, GE Healthcare, Waukesha, WI, USA) using a 32channel head coil (MR Instruments Inc., Minneapolis, MN, USA). A MEPI pulse sequence with SENSElike parallel imaging capabilities (ASSET, GE Healthcare, Waukesha, WI, USA) was developed as shown in Figure 1. The pulse sequence included the static phase correction for ghost removal (Maier et al., 1992) and real-time self-navigated phase correction (Hinks & Xu, 2013) derived from the EPI pulse sequence provided by the manufacturer. The static phase correction maps were derived from non-phase encoded acquisitions (part of the pre-scan) using only the data of the first echo, and the same correction maps were applied to all echoes. For the real-time phase correction, additional echo-navigators were acquired prior to each EPI read-out but only the data acquired prior to the first echo were used to update the static correction maps over time. This was to avoid introducing noise between subsequent echoes due to possible phase estimation errors when correction maps are derived from data acquired with a longer and consequently with a lower SNR. All the images were reconstructed offline using the Orchestra Software Development Kit (GE Healthcare, Waukesha, WI, USA). For each subject, after standard 3-plane localizer and calibration acquisitions, a sagittal high resolution 3D weighted (3D-T1w) volume was acquired using a 3D inversion-recovery fast spoiled gradient echo (3DIR-FSPGR). The following parameters were used: FOV ,

ms,

ms, inversion time

acceleration factor of 2 leading to a resolution of

x

ms, bandwidth x

3

x

mm , matrix

x

x

kHz, parallel imaging (ASSET) x

3

mm and a scan time of 4:48 min. 10

This was followed by a 10-minute resting state acquisition using MEPI, after which the fear conditioning paradigm was performed, also using MEPI. The fMRI acquisitions were performed using MEPI with echoes (

) for both the resting state and the fear conditioning paradigm. Other common parameters

were: field-of-view (FOV)

cm,

mm gap leading to

x

factor 2, bandwidth

kHz,

x x

matrix, fat saturation pre-pulse, slice thickness of 3

mm resolution, 36 slices, parallel imaging (ASSET) acceleration s, flip angle

ms, respectively. The second echo the gold standard for EPI. Finally, the Parameters were: scan time of 2:16 min. A

s,

mm and

,

ms and with

ms leading to ms is considered to be

map was acquired with the same spatial parameters as MEPI. ms,

ms, bandwidth

kHz, flip angle

for a

map allows to measure a voxel-wise frequency shift (in Hz) from the Larmor

frequency. This provides a measurement of how homogenous/inhomogeneous the

field is for a given

voxel, due to the underlying brain induced susceptibility field gradients, and hence it gives a measurement of how much a voxel is affected by signal loss and local echo time shift. Preprocessing Data processing and analyses were done with Statistical Parametric Maps (SPM) version 8 (SPM8 http://www.fil.ion.ucl.ac.uk/spm/software/spm8) and Matlab (Mathworks, Natick, MA, USA). The 3D-T1w and the MEPI data were first converted into the NIFTI-format. Here, we assumed that the different echoes are intrinsically registered at a given time point. This assumption is quite reasonable as compared to

is negligible

. Additionally, spatial transformation was estimated based on one echo (i.e.

) and applied

to all echoes. This way, no noise due to performance difference of the registration algorithms on the different echoes would be introduced to the time courses. The different echoes of MEPI data (i.e.

,

,

) were individually corrected for slice acquisition timing

difference (slice timing correction). For each echo, the first 4 volumes were excluded from the slice timing correction and directly copied for later use (to prevent spill-over effects from slice contrast differences in the initial volumes due to interleaved acquisition). After that, motion correction parameters were estimated using only

data for each time point. This was done by aligning each time point of

volume using a rigid-body model alignment procedure (6 parameters). All echoes

,

to the mean and

were then 11

realigned using the motion parameters estimated from x

and interpolated to the native resolution (

3

mm ). At this point, the mean volumes of the realigned

(excluding the 4 first volumes) and used to estimate the ̂ (

,

and

were computed

̂ ) maps by mono-exponential least

square fitting with the Levenberg–Marquardt algorithm. The weighted sum combined echo computed using equations [(1] and [2]. Next,

x

was then

and the 3D-T1w were co-registered using an affine

transformation (12 parameters) and the same spatial transformations were applied to

,

and

.

Finally, the 3D-T1w was normalized to a T1 weighted template in the MNI space using an affine transformation and non-linear deformations; again the same spatial transformations were applied to ,

and

3

. The 3D-T1w was then interpolated to 1 x 1 x 1 mm and

,

,

,

,

3

to 2 x 2 x 2 mm . For

display purposes, the normalized 3D-T1w of all subjects were averaged to create a study specific T1w template. Moreover, the B0 maps were also co-registered to the 3D-T1w and normalized to MNI space using the spatial transformations estimated for the 3D-T1w. Nuisance removal The nuisance removal steps were performed differently on the fear conditioning data and the resting state data. For both the fear conditioning data and the resting state data, the nuisance removal step consisted of running a general linear model (GLM) with the nuisance regressors only. Then, after the nuisance removal GLM fitting, the residualized images (i.e. the error term ε of the GLM) were used for the statistical analysis described at the next section. In each case, the first 4 volumes were excluded from the nuisance removal GLM. For the fear conditioning data, nuisance regressors included the 6 motion parameters derived from the realignment step, their derivatives, a linear term and a constant term. This led to a total of 14 nuisance regressors. After a high-pass pre-filtering (cut-off 256 s) of the time courses to account for potential slow signal drifts, the nuisance regressors were transformed into Z-scores and included in the nuisance removal GLM. With this preprocessing and nuisance removal strategy, the exact same nuisance regressors were used for the normalized

,

,

and

for the fear conditioning data.

In this fear conditioning task analysis, we only used regressors that were the same across echoes for optimal comparison. Regressors extracted from the white matter (WM) and Cerebrospinal fluid (CSF) 12

compartments are also often incorporated as nuisance regressors. However, as these regressors would naturally differ for

,

,

and

, the “denoising” effects and thus the structure of the residuals (the

error term ε of the GLM) could vary and confound statistical comparisons. Moreover, in robust task fMRI contrasts, the effects of interest are typically large enough to be detected without such additional WM and CSF regressors. Consequently, we did not include these regressors into the nuisance removal GLM. For resting state analyses, however, it is critical to have optimally corrected BOLD signal fluctuations in the region-of-interest and the most specific correlations to other voxels in the brain. This means that for the resting state data, we included the echo-specific WM and CSF nuisance regressors. This was performed with the CompCor method (Behzadi et al., 2007). For this purpose, a CSF mask, a WM mask and a mask containing the voxels with the upper 2% standard deviation over time (2pSTD) were created from the normalized data. The CSF mask was computed by thresholding the ̂ map generated during the preprocessing according to Kundu et al. (2015) with

( Here, we exploit the fact that the

)

of the CSF is significantly longer than that of other tissue types

(Kundu et al., 2015). The resulting mask is then slightly eroded using a small structural element (i.e. 2D kernel defined as [0 1 0; 1 1 1; 0 1 0] in Matlab notation) to be sure that our final mask corresponds to more centrally located CSF voxels with no partial volume effect from neighboring gray matter (GM) voxels. The WM mask was computed by the intersection of two different primary masks. The first mask was obtained by segmentation of the first time point of

(the worst case for signal dropouts in

inhomogeneous brain regions but with a good white/grey matter contrast). The second mask was a brain mask computed by means of mathematical morphology on the mean (over time) of

. Then, the

intersection of the two masks was computed and the resulting mask was slightly eroded as well. With this approach, the resulting mask had a high likelihood of containing only pure WM voxels common to and, by construction, 2pSTD mask of , and

and

,

and

,

,

. The 2pSTD mask was simply computed by taking the intersection of the . In this way, the three masks (WM, CSF, 2pSTD) masks can be used on

. Next, these masks were used for extracting the corresponding time courses in

,

,

,

. From each a separate principal component analysis (PCA) was performed on these 13

compartments, then the three first principal components (PCs) were computed and added to the nuisance removal GLM. In summary, the nuisance regressors consisted of 14 generic regressors (6 motion parameters, 6 motion derivatives, 1 linear, 1 constant) and 9 echo-specific regressors (3 PCs from the CSF mask, 3 PCs from the WM mask and 3 PCs from the 2pSTD mask) resulting in total of 23 regressors. As performed for the fear conditioning data, a high-pass pre-filtering (cut-off 256 s) of the time courses was included, the nuisance regressors were transformed into Z-scores and included in the nuisance removal GLM (the first 4 volumes were excluded). After the nuisance removal steps, the times courses were variance normalized. Then the seed (defined from fear conditioning task analysis as clarified later) time courses were extracted. The resting-state time-courses were then filtered with a bandpass filter (cutoff frequencies [0.01 - 0.1] Hz).

Finally,

,

,

and

were smoothed after nuisance removal using a Gaussian kernel of 6 x 6 x 6 mm

3

at full width at half maximum for both the resting state and the fear conditioning data. Note that all GLMs were computed using an explicit whole brain mask in MNI space as distributed with SPM8 (i.e. for nuisance removal, first and second levels). Brain regions were identified and labelled according to the Automated Anatomical Labeling (AAL) atlas and toolbox (Tzourio-Mazoyer et al., 2002).

Statistical analysis The first level models (i.e. GLM) for fear conditioning were applied after nuisance removal and comprised the three stimulus regressors (the safety stimulus and

which includes

) convolved with the canonical hemodynamic response function. After hyperparameter

estimation, the contrasts vector:

and the conditioned stimuli

(contrast vector:

) and its reverse

(contrast

) were computed for each subject.

Second level random effects analyses were performed for statistical inference on the group level. First, a one-sample t-test was computed for and

and its reverse

independently. To perform a comparison of

,

,

to check the results of and

,

,

, a full factorial analysis of variance

(ANOVA) design was used with one factor “method” with 4 dependent levels, representing . The main effects contrasts were computed to map the brain regions where

,

,

, and

,

and differ. 14

Finally, a paired t-test was computed to map

for the contrast

and its reverse

. To control for multiple testing, familywise error (FWE) correction was applied. Statistical maps were sampled at

uncorrected and clusters surviving

cluster-wise (denoted

) were deemed significant. For the sake of completeness, the paired t-test also computed for the contrast

and its reverse

and

were

.

The behavioral ratings were entered into a repeated measures ANOVA with the two within-subject factors stimulus (3 levels) and time (3 levels). Successful learning of stimulus contingencies is indicated by a significant stimulus × time interaction, with post-hoc (paired) t-tests as illustrative of the direction of the effect. When the assumption of sphericity was not met, Greenhouse-Geisser corrections were performed (but note that original df-values are provided for readability purposes). Again, significance was set to a level of

.

Seed-based analysis Based on the fear conditioning results, a seed region (6mm radius sphere) was derived from the vmPFC cluster from the paired t-test

. The sphere was centered at the centre of mass of the vmPFC

cluster and the same vmPFC seed was used for all subjects. Using this seed in the vmPFC, the mean signal time course was extracted from the

,

,

and

data (after nuisance removal but before spatial

smoothing) and included in the first level model (GLM) together with a constant regressor. The contrast for the regressor of interest (i.e. the seed) was then calculated. On the second level, a one sample t-test design was used to compute the results of of voxel-wise

(denoted

,

,

and

independently. Clusters surviving a threshold

) and a threshold on cluster extent of

were

deemed significant.

15

RESULTS Four datasets had to be discarded due to incomplete data, leaving 32 subjects for the analysis. )

The behavioral data showed successful fear learning. There was a significant effect of stimulus, ,

, a significant effect of time, )

stimulus × time interaction, revealed that

)

,

, all t-values

,

, and critically, a significant

. After fear learning, post-hoc paired t-tests , all p-values

.

The results of the second level analysis of the fear conditioning task (contrasts ) for

,

,

and

and

computed independently are displayed in Figure 3 (corresponding brain regions

and peak voxels are provided in Supp. Table 1). It is worth noting that in Figure 3, the vmPFC cluster in the contrast voxel of

appeared more ventral in

was

than in

and of course

. As expected, the peak

mm more ventral in the MNI z-direction than the peak voxel of

. This was expected,

as reminded in the theory section, in the vmPFC, the local echo time (e.g. local to a voxel or brain region) is shifted to lower values due to brain induced susceptibility gradient as demonstrated by Deichmann et al. (2002). Additionally, still in the vmPFC, the more ventral the considered voxel is, the more shifted the local echo time. In the most ventral parts of the vmPFC, the BOLD signal was better sampled by the shortest echo (

with a short

ms), while in the upper part of the vmPFC it was better sampled by

ms). Also, as expected,

appeared to yield a stronger response (i.e. higher T-values) in the

vmPFC cluster with a larger extent (the vmPFC cluster size was voxels for

and

voxels for

(with

voxels for

). In the reverse contrast

survive the cluster-wise FWE correction for

but only for

,

,

voxels for

,

, the cluster in dACC did not and

.

16

Figure 3: Whole brain second level analysis of the fear conditioning task for The contrast

is shown in hot colors and the reverse contrast (

are displayed in MNI space for the four methods

,

,

and

,

,

and

(one sample t-tests).

) is shown in cold colors. T-values

. All clusters shown survived

, the

crosshair is positioned such as the vmPFC, dACC and insula clusters are visible (S = superior, A = anterior, P = posterior, L = Left).

A full factorial ANOVA addressed the question for which brain regions the methods

,

,

and

differed statistically and as expected, they strongly differed in the vmPFC (see Figure 4 and Supp. Table 2 for a detailed brain region list). The main effects (similar to Figure 4) are also provided with voxel-wise FWE correction in Supp. Figure 1 with corresponding brain regions detailed in Supp. Table 3. The peakvoxel’s contrast estimates (Figure 4-b) revealed a large contribution of

in the vmPFC cluster (Figure 4-a

and Figure 4-b). Interestingly, the four methods also differed strongly in the control regions (dACC and bilateral insula). Here, however, the contrast estimates of the control regions (located in homogenous brain regions) revealed the largest contribution from that

compared with

and

(Figure 4-b). This shows

revealed the strongest differences in both vmPFC and the control regions; taking advantage of

early echoes (i.e.

) in susceptibility affected regions (vmPFC) and of later echoes (i.e.

) in control

regions (dACC and insula) not affected by susceptibility. To illustrate this finding, the mean weights in these brain regions are provided in Table 1 (see also Supp. Figure 2). There was a robust contribution of

17

(

) in the vmPFC and of

approximately the same for

in the control regions (

in all brain regions (

Figure 4: Main effect of the factor "methods" (i.e.

,

); however, the contribution was

).

,

and

) and contrasts estimates for the peak voxels of

various clusters. (a) Main effect of the factor "methods",

, in MNI space. (b) Contrast estimates and 90%

confidence intervals of the peak voxel in vmPFC (contrast

), dACC and right insula (contrast

),

with MNI coordinates provided below (S = superior, A = anterior, P = posterior, L = Left).

Table 1: Contribution (mean weights generate

standard deviation across subjects) of each echo (

,

,

) used to

in the targeted and control regions visible in the main effects analysis (see Figure 4-a).

Cluster

vmPFC

dACC

Insula (R)

Insula (L)

18

To relate these results to physical properties, the clusters shown on Figure 4-a (vmPFC, dACC and bilateral insula regions) were used to extract the mean frequency shifts (in Hz) from the subjects. The mean frequency shift was cluster,

Hz (mean

Hz in the dACC cluster,

maps for all

standard deviation) in the vmPFC

Hz in right insula cluster and

Hz in the

left insula cluster. As expected, the vmPFC cluster showed the highest frequency shift and the other clusters showed relatively small frequency shifts in comparison. This confirmed that despite the high frequency shift in susceptibility affected regions, Next, we tested on the

and

vmPFC compared to the gold standard stronger task response than

had optimal contrast estimates. contrasts whether

showed increased activity in

with a paired t-test. As anticipated,

showed a significantly

(see Figure 5-a and Supp. Table 4 for detailed brain regions).

Furthermore, in the reverse comparison,

also showed a significantly stronger (negative) task response

in the bilateral insula (see Figure 5-b). However, no significant differences were found in the dACC (using the same threshold of paired t-test

for the contrast

as in other analyses). The cluster identified in the vmPFC with the (Figure 5) represents an additional cluster extent (ventral

extension), which could not have been identified with

(the gold standard) as there is no signal with

(see Supp. Figure 3 and Supp. Figure 4 for an illustration of this additional cluster extent detected with ).

19

Figure 5: Paired t-test

for the contrast

(hot colors) and the reverse contrast

colors) displayed in MNI space. All clusters shown were significant at

(cold

. (a) Paired T-test

centered at the center of mass of the vmPFC cluster. (b) same as in (a) but displayed in a multislice view to show the significant difference in the vmPFC and the bilateral insula. No significant differences were found in the dACC (S = superior, A = anterior, P = posterior, L = Left).

For the sake of completeness, the comparisons comparisons

,

and

Figure 6. As expected, the comparison

and

were also performed. The three

, with a focus on the same regions of interest, are given in showed very few significant differences in the vmPFC

(only in the extreme upper part of the region) but showed significant differences in other brain regions (bilateral insula and dACC – see Figure 6). The comparison

highlighted the difference in the

vmPFC, whereas no difference was observed in the other regions (Figure 6). This is understandable as provides the best

to

to sample the BOLD signal in the vmPFC (hence

in the vmPFC), while with found for

most of vmPFC completely dropped out (hence the large vmPFC cluster

). For the bilateral insula and dACC,

signal in these regions while

shows little difference

provides the best

to

to sample the BOLD

has very little BOLD sensitivity in these regions (hence the large cluster 20

found for In other words,

). Notably, the comparison

,

,

did not yield any significant clusters.

derived from MEPI always performs at least as good as the best

region but provides the best

for a given brain

for the whole brain.

21

Figure 6: Paired t-test

,

and

for the contrast

(hot colors) and the reverse contrast

(cold colors) displayed in MNI space for the regions of interest (vmPFC, dACC and bilateral insula). All clusters shown were significant at vmPFC region. (b) Paired t-tests

. (a) Paired t-tests and

,

and

with a focus on the

with a focus on the bilateral insula region. (c) Paired t-tests

22

and

with a focus on the dACC region. The comparisons

are not shown for the bilateral insula

and the dACC as there were no significant differences (S = superior, A = anterior, P = posterior, L = Left).

To ensure that the vmPFC cluster identified by the paired t-test

(Figure 5-a) represents relevant

fMRI signal (rather than containing noise/artifacts or other unspecific signal fluctuation unrelated to neural activity), it should be functionally connected to the bilateral amygdala in line with previously published anatomical findings. For the resting state analysis, the seed was derived as a 6-mm sphere centered at the center-of-mass (

in MNI coordinates) of the vmPFC cluster observed in the paired t-test

(Figure 7).

Figure 7: Seed of the vmPFC cluster derived from the paired t-test results (see Figure 5) and used for the resting state analysis. This is a 6-mm sphere centered at [2 41 -20] in MNI coordinates (S = superior, A = anterior, P = posterior, L = Left).

Figure 8 depicts the results of the vmPFC seed based resting state analysis for but not the single echoes

,

,

,

,

and

. Only

,

, showed significant connectivity between the vmPFC seed and the left

and right amygdala, as expected (see blue arrows on Figure 8, lowest panels). Moreover, when considering the whole brain pattern, the vmPFC seed analysis revealed typical nodes of the DMN, particularly the posterior cingulate cortex and adjacent precuneus, and the bilateral inferior parietal lobules. While already partially present in

and

, this pattern was most pronounced in

, providing

indirect evidence at the group level that the vmPFC participates in the DMN.

23

Figure 8: Results of the vmPFC seed based analysis (one sample t-tests) on resting state data computed for and

. All clusters shown survived

. Note that

,

,

shows significant connectivity between the vmPFC

seed and the left and right amygdala (blue arrows).

24

Figure 9: Mean contrast estimates and standard error of the left (L) and right (R) amygdala ROI from the AAL atlas in the vmPFC seed based analysis of Figure 8. (a) Mean contrast estimates and standard error for Mean contrast estimates and standard error for

and

,

,

and

. (b)

together with the results of a paired t-test.

Finally, an atlas based region-of-interest (ROI) based comparison was performed to examine whether our vmPFC seed was more strongly connected to amygdala in

than

as suggested by the normative

maps of the seed analysis (Figure 8). The mean contrast in left and right amygdala (as defined by the AAL atlas) was extracted for all subjects and all methods. The values of

,

,

and

correlation between the vmPFC seed and the bilateral amygdala was strongest for

showed that the (see Figure 9-a). A

repeated measures ANOVA on these mean values showed a significant main effect of “method” (i.e. of ,

and

) for the left amygdala (

)

). The paired t-test comparison

) and the right amygdala (

,

)

(see Figure 9-b) showed that the contrast estimates of

were significantly stronger than the current gold standard

(

for the left amygdala,

for

the right amygdala). In addition to the voxel-wise results, these regional results corroborate the fact that the vmPFC cluster found in the paired t-tests the right amygdala, and more strongly in

(Figure 5-a) is functionally connected to the left and

than in

.

DISCUSSION In this study, a MEPI pulse sequence was developed and compared with the currently accepted gold standard (as represented by the second echo

with

ms) in a sufficiently powered study using a 25

fear conditioning task known to elicit activity in the vmPFC (our main targeted region) and the dACC and bilateral insula (as control regions). In summary, our observations support the superiority of MEPI in detecting activity in vmPFC, a region strongly affected by susceptibility field gradients. More specifically, a direct comparison was first performed and a statistically stronger task response was found for MEPI (i.e. ) in the vmPFC. Secondly, a seed-based resting state functional connectivity analysis in a separate data-set revealed the expected connectivity of the vmPFC to the bilateral amygdala. Finally, this susceptibility affected region (i.e. vmPFC), showed markedly high frequency shifts and were heavily weighted by

, whereas the homogeneous brain regions (i.e. dACC and bilateral insula) showed small

frequency shifts and was heavily weighted by A priori, the sample size (

.

) of our study was based on achieving sufficient power (

medium effect size in the vmPFC with a standard alpha of

) to detect a

in a one-tailed paired samples

comparison. Because the dropout in our study was lower than expected, and the observed effect size for the comparison of

of the vmPFC cluster was larger (Cohen’s

revealed that our study had sufficient power ( This also applies to the ). Only for the comparison resulting in a post-hoc power of

), post-hoc power analyses

) to detect these differences in the same comparison.

comparison, for which we also observed a large effect size (power we observed a slightly lower effect size (Cohen’s

),

. This analysis revealed that significance in our study is unlikely to be

due to an overestimated effect sizes due to random sampling variability (Button et al., 2013). To ensure a fair comparison of the BOLD information contained in the different datasets (

,

,

and

), several study design choices had to be made. On the acquisition side, achieving whole brain coverage with reasonable spatial and temporal resolution (3-4 mm isotropic in 2-2.5s) led us to a isotropic resolution with a temporal resolution of

mm

s. This is slightly slower than what can be achieved

with a conventional single echo EPI pulse sequence (e.g,

mm isotropic in

s) but it does not deviate

unacceptably from what many fMRI studies have employed. However, the state of the art is to use multiband excitation which allows to reduce the

and hence to get an increased temporal resolution

(Moeller et al., 2010; Setsompop et al., 2012; Todd et al., 2016; Xu et al., 2013). Since the standard EPI (i.e.

) was embedded in MEPI, it allowed a direct comparison of the two methods. Here, a direct

comparison meant that the comparison was performed on the same signals, sampled at approximately the 26

same time, for the same subjects in the same brain state, the time difference between the subsequent echoes being negligible compared to the hemodynamic response and the

. For the analyses, a

conservative GLM design was used for both the task and the seed-based resting state analysis. Yet to perform a fair comparison, the utmost attention was paid to the preprocessing steps in order to avoid introducing any bias or undesired variability, as different registration and normalization algorithms could induce variability between the echoes due to the different contrast (Gonzalez-Castillo et al., 2013). For this reason, all registration and normalization steps were performed on the first echo and the resulting spatial parameters were applied to the subsequent echoes. In this way, all echoes underwent the exact same spatial transformation. For the fear conditioning task, the exact same nuisance regressors and the same design matrix were used for all echoes for a better comparability; the data were also variance normalized for the same reason. The task-induced increases (and decreases) in activity were robust enough to be detected without any echo-specific nuisance regressors which might have introduced a confounding factor for comparisons among echoes. Still, the situation was different for the resting state analysis as, in this case, it is critical to have spontaneous BOLD signal fluctuations with optimal correction in the region-ofinterest and the most specific correlations to other voxels in the brain. Consequently, additional echo specific nuisance regressors were used for resting state functional connectivity analysis. One issue that might have created a bias is that the signal-to-noise ratio (SNR) of

and

were not matched.

Effectively, an inherent advantage of the normalized weighted sum approach (i.e.

) used in this study for

MEPI is an SNR enhancement (Poser et al., 2006). However, this is almost impossible to get an SNR matched EPI (

) and MEPI ( ) without introducing uncontrollable run-to-run variability; particularly

regarding a fear conditioning task with strong habituation and retest effects. Recently, Kirilina et al. (2016) compared several advanced pulse sequences (including MEPI) for fMRI; they confirmed that (i) advanced pulse sequences increased the sensitivity at the single subject level and found that (ii) the sensitivity of these advanced pulse sequences at the group level was mainly limited by the inter-subject variability – with the exception of MEPI, for which better performance in the orbitofrontal cortex was reported, although this was not statistically tested at the group level. Our results are in accordance with this work (Kirilina et al., 2016) regarding the MEPI part (regarding inter-subject variability see also Supp. Figure 5); but due to the different objectives of these two studies, several differences in the 27

study design should be highlighted. The first was that their data were acquired in 4 different runs while our data were acquired in a single run. This implies we have the same intra and inter-subject variability between the vs

,

,

and

, which allows us to draw our conclusion on the different methodologies (

) without introducing uncontrollable factors due to multiple runs. The second difference concerns the

task; Kirilina et al. (2016) did not specifically focus on the vmPFC and used emotional learning and reward based learning tasks which may not elicit activity in the ventral parts of the mPFC according to several meta-analyses (Garrison et al., 2013; Hayes et al., 2014; Sabatinelli et al., 2011; Zhang et al., 2013). We did focus on the vmPFC specifically and the fear conditioning task used in this work does elicit the vmPFC in a robust way (Fullana et al., 2015) which makes it especially suitable for a comparison of different pulse sequences when focusing on the vmPFC specifically. A third relevant difference was the MEPI acquisition parameters. Kirilina et al. (2016) used five echoes with higher acceleration factors (i.e. 3) and partial Fourier (i.e. 6/8) leading to

-

of (

ms,

ms,

ms,

ms,

ms). For this study, we used a

moderate acceleration factor of 2 without partial Fourier, leading to a narrower range of =

/

/

(i.e.

ms). The choice of MEPI acquisition parameters we made was mainly constrained by the

desired whole brain acquisition with reasonable temporal resolution (3-4mm isotropic in 2-2.5s) and embedding the standard EPI (

) in a single run. Of course, further optimization of MEPI to get a shorter

is possible and may have beneficial effects for signal detection in the vmPFC. The fourth difference relates to the echo combination technique: whereas we have based the weighing to generate

on the

relative, expected BOLD contrast contribution of each echo (Poser et al., 2006), Kirilina et al. (2016) made the weighting according to the tSNR (temporal signal-to-noise ratio) of each echo. The different weighting techniques should not make a large difference, as previously reported (Poser et al., 2006; Kettinger et al., 2016). Kettinger et al. (2016) compared four different weighting options using MEPI (with only two echoes) and a reward based task. They found that the optimal combination (normalized weighting sum similar to the approach we have used) was better than simple averaging at subject level but that these advantages do not carry over at the group level even not in vmPFC. They concluded that at the group level the simple averaging approach performed as well as the optimal combination technics and that this is likely because the main limiting factor at group level remains the inter-subject variability according to Kirilina et al. (2016). This tends to show that the different weighting techniques should not make a large difference. We have directly tested for the difference between

and

and our results are in agreement with previous studies 28

(Kirilina et al., 2016; Poser et al., 2006); MEPI shows improved performance in susceptibility-affected brain regions such as vmPFC. It is worth mentioning that our results also differed for brain regions not affected by susceptibility (dACC and insula regions), for which MEPI ( performance due to the contribution of

) still shows slightly better

.

In the present study, we chose to use a rather conservative analysis pipeline to compare the different echoes. Consequently, we had to exclude more advanced analysis techniques that have shown improved sensitivity like multi-echo independent component analysis (ME-ICA) (Kundu et al., 2012) and multi-echo independent component regression (ME-ICR) (Kundu et al., 2013; Kundu et al., 2015). Indeed, ME-ICA and ME-ICR can only be applied to the MEPI data and cannot be applied to the single echo data, which would decrease the fairness of the comparison for our research purposes. However, ME-ICA and ME-ICR are expected to further improve the sensitivity as these advanced analysis pipelines have been shown to give superior performance compared to the weighted sum approaches (like

) on MEPI resting state data

(Kundu et al., 2012; Kundu et al., 2013; Kundu et al., 2015) and also on MEPI task data (Lombardo et al., 2016). For this study, we aimed to compare the information contained in the data but not to optimize the MEPI analysis pipelines. For the fear conditioning data, there are several options to enhance and optimize the analyses. One option could have been to use the respiration waveform with RETROICOR (Glover et al., 2000). The methodology would have provided additional echo-independent regressors for the nuisance removal step, which would have allowed a better denoising of the data. Yet, the effect of respiration could be assumed to be approximately the same on

,

,

and

and should not affect the

main effect (Figure 4) and the paired t-test comparison

(Figure 5 and Figure 6) to a large extent.

Another option to improve the nuisance removal step would have been to make a “multi-echo CompCor” (as suggested by one of the reviewer) by performing a PCA across all echoes as well as voxel locations for the WM, CSF and 2pSTD mask. However, such a novel analysis strategy is not yet validated and the influence of WM and CSF fluctuations on the different echoes (as a function of TE) is still unclear. This could lead to privileging one echo versus one other. This topic requires further study to examine and validate a generic solution applicable to MEPI on all tasks. A third option to optimize the analysis pipeline would be to use the

field maps to correct the EPI distortion of the functional data. However, with the

specific EPI acquisition parameters and specific scanner hardware (in particular, the gradient system) 29

employed in this study, the EPI geometric distortion was fairly limited (i.e.,

voxel). So the spatial

smoothing of 6 mm at FWHM (~1.8 times the voxel size) sufficed to reduce the risk of significant mislocalization in the vmPFC cluster. We did not apply the field map based distortion correction since the precise localization of our clusters was not so relevant for our main objective: to compare standard EPI). Given that i) and

, iii)

,

and

,

and

were acquired almost simultaneously, ii)

with

was derived from

with

,

are affected by the EPI distortion in the exact same way (but with different

levels of signal loss), and iv) the fact that we have applied the exact same spatial transformation to and

(the

, any minor mislocalization will be the same for

,

,

,

,

. Consequently, the comparison of

and our results and conclusions remain valid, even if some voxels would be slightly mislocalized in

the vmPFC. However, all the distortion correction methods developed for EPI can be directly applied to different echoes of MEPI if necessary, for instance with higher fields, or at higher resolution or with a less powerful gradient system. MEPI can be further improved to increase spatial and/or temporal resolution, to reach shorter

and to

acquire more echoes. As already mentioned, MEPI can be used with high parallel imaging acceleration factor and partial Fourier, which would allow for the reduction of the length of the individual echo-train and hence reaching shorter

and to acquire more echoes (Kirilina et al., 2016; Poser et al., 2006; Posse et

al., 1999). MEPI is compatible with simultaneous multi-slice (also called multiband) acquisition techniques like blipped CAIPI (Setsompop et al., 2012). Depending on the chosen multiband factor, this will allow a reduction of the

by approximately the same factor as shown at 3T (Olafsson et al., 2015) and at 7T

(Boyacioglu et al., 2015). Simultaneous multi-slice multi-echo EPI can be combined with ME-ICA and MEICR for improved sensitivity and specificity with the same spatial resolution (i.e. temporal resolution (e.g.

-

mm) and a high

s) (Olafsson et al., 2015; Boyacioglu et al., 2015). The advantage of

simultaneous multi-slice acquisition in the context of MEPI lies in improving the trade-off between spatial and temporal resolution. Simultaneous multi-slice MEPI can be combined with thin slice (i.e. 1-2 mm) acquisition technique to further reduce the signal drop in susceptibility-affected brain regions (Kim et al., 2016). The performance of simultaneous multi-slice and parallel imaging, in terms of image quality and achievable acceleration factor, will strongly depend on the design and number of channels of the phased-

30

array coil. Currently available 32-channel coils already yield promising results, but optimized coils will likely allow pushing the limit further. When using a MEPI pulse sequence (with or without simultaneous multi-slice), the key point is to sample an appropriate range of vmPFC, the

s. To be able to recover the susceptibility affected brain regions like the

should be sufficiently short (Deichmann et al., 2002). In this work, we detected robust

activation at the group level using a with a shorter

of

the theoretically optimal

of

ms, but other authors have also reported successful results

ms (Kirilina et al., 2016). For the longest (

). At 3T, the

, it may not be worth sampling beyond

of grey matter is usually considered to be approximately

ms (Norris, 2006), although some authors have reported it to be around this study, the longest

was set to

ms (Peters et al., 2007). In

ms, which is close to the optimal

made a similar choice (Kirilina et al., 2016). Using a longer

. Other authors have

may not add much value and would

compromise the spatio-temporal resolution. Considering the duration of the EPI read-out (

ms in this

work) and the local echo-time as defined by Deichmann et al. (2002), the parameters used in this study allow covering of a wide range of echo-times. Regarding the number of echoes, we would recommend using at least , as this allows to have a short of

for regions affected by susceptibility artefacts, a middle

ms as a standard echo time for fMRI at 3T, and a longer

close to the optimal

(

). Of

course, whenever possible, acquiring more echoes would be more beneficial as long as the abovementioned constraints on the shortest and longest

are fulfilled.

CONCLUSION In this study, we developed a MEPI pulse sequence and compared it with standard EPI in a fear conditioning task. Using a normalized weighted sum approach, MEPI showed superior performance over standard EPI in detecting activation in the vmPFC at the group level, better performance in the bilateral insula and no significant difference (at the employed thresholds) in the dACC. This suggests that MEPI is better in susceptibility affected brain regions and performs as well as or better than the gold standard in control regions. MEPI performs always at least as good as the best

for a given brain region but

provides all necessary echo times to get an optimal BOLD sensitivity for the whole brain. Then, to demonstrate that the additional cluster extent (ventral extension identified in the vmPFC with the paired t31

test

for the contrast

, which corresponds to a region with no signal in

) detected in the

vmPFC with MEPI is a functionally relevant region, resting state functional connectivity analysis from the same subjects were used to demonstrate that this region was connected both to the amygdala (and more strongly with

than standard EPI) and part of the DMN. For both, the connectivity with the vmPFC was

expected from independent imaging modalities (e.g. PET) and anatomical knowledge. Our findings demonstrate that MEPI can be used effectively in the field of affective neuroscience with minimum cost to the spatio-temporal resolution. As the vmPFC is a key region for many affective regulation processes including fear learning and extinction, we expect these results to encourage a more systematic use of MEPI in affective neuroscience and psychiatry. This appears especially relevant for tasks or studies that rely on proper signal coverage of the vmPFC. ACKNOWLEDGMENTS The authors would like to thank Eric Printz for his help on the image reconstruction code and Ines Eidner for her help during data acquisition. VS acknowledges financial support from the Bavarian Academy of Sciences and Humanities. The authors also thank Catriona Wimberley for her useful comments and English corrections. We would like to thank the anonymous reviewers for the overall positive feedback and their helpful remarks and suggestions which help to enhance the manuscript.

REFERENCE Andrews-Hanna, J.R., Saxe, R. & Yarkoni, T., 2014. Contributions of episodic retrieval and mentalizing to autobiographical thought: evidence from functional neuroimaging, resting-state connectivity, and fMRI meta-analyses. Neuroimage, 91, pp.324-35. Available at: http://dx.doi.org/10.1016/j.neuroimage.2014.01.032. Behzadi, Y., Restom, K., Liau, J. & Liu, T.T., 2007. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. Neuroimage, 37(1), pp.90-101. Available at: http://dx.doi.org/10.1016/j.neuroimage.2007.04.042. Boyacioglu, R. et al., 2015. Improved sensitivity and specificity for resting state and task fMRI with multiband multi-echo EPI compared to multi-echo EPI at 7T. Neuroimage, 119, pp.352-61. Available at: http://dx.doi.org/10.1016/j.neuroimage.2015.06.089. Buckner, R.L., Andrews-Hanna, J.R. & Schacter, D.L., 2008. The brain's default network: anatomy, function, and relevance to disease. Ann N Y Acad Sci, 1124, pp.1-38. Available at: http://dx.doi.org/10.1196/annals.1440.011. 32

Button, K.S. et al., 2013. Power failure: why small sample size undermines the reliability of neuroscience. Nature reviews. Neuroscience, 14(5), pp.365-76. Chen, N. & Wyrwicz, A.M., 1999. Removal of intravoxel dephasing artifact in gradient-echo images using a field-map based RF refocusing technique. Magn Reson Med, 42(4), pp.807-12. Cho, Z.H. & Ro, Y.M., 1992. Reduction of susceptibility artifact in gradient-echo imaging. Magn Reson Med, 23(1), pp.193-200. Constable, R.T. & Spencer, D.D., 1999. Composite image formation in z-shimmed functional MR imaging. Magn Reson Med, 42(1), pp.110-17. Cusack, R. et al., 2005. An evaluation of the use of passive shimming to improve frontal sensitivity in fMRI. Neuroimage, 24(1), pp.82-91. Available at: http://dx.doi.org/10.1016/j.neuroimage.2004.08.029. Deichmann, R. et al., 2002. Compensation of susceptibility-induced BOLD sensitivity losses in echoplanar fMRI imaging. Neuroimage, 15(1), pp.120-35. Available at: http://dx.doi.org/10.1006/nimg.2001.0985. Delgado, M.R., Nearing, K.I., Ledoux, J.E. & Phelps, E.A., 2008. Neural circuitry underlying the regulation of conditioned fear and its relation to extinction. Neuron, 59(5), pp.829-38. Available at: http://dx.doi.org/10.1016/j.neuron.2008.06.029. Drevets, W.C. et al., 1997. Subgenual prefrontal cortex abnormalities in mood disorders. Nature, 386(6627), pp.824-27. Available at: http://dx.doi.org/10.1038/386824a0. Duits, P. et al., 2015. Updated meta-analysis of classical fear conditioning in the anxiety disorders. Depress Anxiety, 32(4), pp.239-53. Available at: http://dx.doi.org/10.1002/da.22353. Etkin, A. & Wager, T.D., 2007. Functional neuroimaging of anxiety: a meta-analysis of emotional processing in PTSD, social anxiety disorder, and specific phobia. Am J Psychiatry, 164(10), pp.1476-88. Available at: http://dx.doi.org/10.1176/appi.ajp.2007.07030504. Farzaneh, F., Riederer, S.J. & Pelc, N.J., 1990. Analysis of T2 limitations and off-resonance effects on spatial resolution and artifacts in echo-planar imaging. Magn Reson Med, 14(1), pp.123-39. Faul, F., Erdfelder, E., Lang, A.-G. & Buchner, A., 2007. G*Power 3: a flexible statistical power analysis program for the social, behavioral, and biomedical sciences. Behavior research methods, 39(2), pp.17591. Fullana, M.A. et al., 2015. Neural signatures of human fear conditioning: an updated and extended metaanalysis of fMRI studies. Mol Psychiatry, Available at: http://dx.doi.org/10.1038/mp.2015.88. Fullana, M.A. et al., 2015. Neural signatures of human fear conditioning: an updated and extended metaanalysis of fMRI studies. Mol Psychiatry, Available at: http://dx.doi.org/10.1038/mp.2015.88. Garrison, J., Erdeniz, B. & Done, J., 2013. Prediction error in reinforcement learning: a meta-analysis of neuroimaging studies. Neurosci Biobehav Rev, 37(7), pp.1297-310. Available at: http://dx.doi.org/10.1016/j.neubiorev.2013.03.023. Gati, J.S., Menon, R.S., Ugurbil, K. & Rutt, B.K., 1997. Experimental determination of the BOLD field strength dependence in vessels and tissue. Magn Reson Med, 38(2), pp.296-302.

33

Glover, G.H., 1999. 3D z-shim method for reduction of susceptibility effects in BOLD fMRI. Magn Reson Med, 42(2), pp.290-99. Glover, G.H. & Law, C.S., 2001. Spiral-in/out BOLD fMRI for increased SNR and reduced susceptibility artifacts. Magn Reson Med, 46(3), pp.515-22. Glover, G.H., Li, T.Q. & Ress, D., 2000. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magnetic resonance in medicine, 44(1), pp.162-67. Gonzalez-Castillo, J. et al., 2013. Effects of image contrast on functional MRI image registration. Neuroimage, 67, pp.163-74. Available at: http://dx.doi.org/10.1016/j.neuroimage.2012.10.076. Greicius, M.D. et al., 2007. Resting-state functional connectivity in major depression: abnormally increased contributions from subgenual cingulate cortex and thalamus. Biol Psychiatry, 62(5), pp.429-37. Available at: http://dx.doi.org/10.1016/j.biopsych.2006.09.020. Griswold, M.A. et al., 2002. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn Reson Med, 47(6), pp.1202-10. Available at: http://dx.doi.org/10.1002/mrm.10171. Groenewold, N.A. et al., 2013. Emotional valence modulates brain functional abnormalities in depression: evidence from a meta-analysis of fMRI studies. Neurosci Biobehav Rev, 37(2), pp.152-63. Available at: http://dx.doi.org/10.1016/j.neubiorev.2012.11.015. Hayes, D.J., Duncan, N.W., Xu, J. & Northoff, G., 2014. A comparison of neural responses to appetitive and aversive stimuli in humans and other mammals. Neurosci Biobehav Rev, 45, pp.350-68. Available at: http://dx.doi.org/10.1016/j.neubiorev.2014.06.018. Hayes, J.P., Hayes, S.M. & Mikedis, A.M., 2012. Quantitative meta-analysis of neural activity in posttraumatic stress disorder. Biol Mood Anxiety Disord, 2, p.9. Available at: http://dx.doi.org/10.1186/2045-5380-2-9. Hinks, R.S. & Xu, D., 2013. System and method of image artifact reduction using self-navigated real-time phase correction in echo planar imaging. US Patent 8,483,457. Hsu, J.-J. & Glover, G.H., 2005. Mitigation of susceptibility-induced signal loss in neuroimaging using localized shim coils. Magn Reson Med, 53(2), pp.243-48. Available at: http://dx.doi.org/10.1002/mrm.20365. Juchem, C. et al., 2014. DYNAmic Multi-coIl TEchnique (DYNAMITE) shimming of the rat brain at 11.7T. NMR Biomed, 27(8), pp.897-906. Available at: http://dx.doi.org/10.1002/nbm.3133. Kettinger, A., Windischberger, C., Hill, C. & Nagy, Z., 2016. Advanced combinations of dual-echo fMRI data provide no advantages over the simple average at group-level analyses. In Proceedings 24th Scientific Meeting, International Society for Magnetic Resonance in Medicine, Singapore., 2016. Kim, T., Zhao, T. & Bae, K.T., 2016. Enhancement of functional MRI signal at high-susceptibility regions of brain using simultaneous multiecho multithin-slice summation imaging technique. J Magn Reson Imaging, 44(2), pp.478-85. Available at: http://dx.doi.org/10.1002/jmri.25170. Kirilina, E. et al., 2016. The quest for the best: The impact of different EPI sequences on the sensitivity of random effect fMRI group analyses. Neuroimage, 126, pp.49-59. Available at: http://dx.doi.org/10.1016/j.neuroimage.2015.10.071. 34

Koch, K.M., Brown, P.B., Rothman, D.L. & de Graaf, R.A., 2006. Sample-specific diamagnetic and paramagnetic passive shimming. J Magn Reson, 182(1), pp.66-74. Available at: http://dx.doi.org/10.1016/j.jmr.2006.06.013. Koch, K.M., Rothman, D.L. & de Graaf, R.A., 2009. Optimization of static magnetic field homogeneity in the human and animal brain in vivo. Prog Nucl Magn Reson Spectrosc, 54(2), pp.69-96. Available at: http://dx.doi.org/10.1016/j.pnmrs.2008.04.001. Kundu, P. et al., 2015. Robust resting state fMRI processing for studies on typical brain development based on multi-echo EPI acquisition. Brain Imaging Behav, 9(1), pp.56-73. Available at: http://dx.doi.org/10.1007/s11682-014-9346-4. Kundu, P. et al., 2013. Integrated strategy for improving functional connectivity mapping using multiecho fMRI. Proc Natl Acad Sci U S A, 110(40), pp.16187-92. Available at: http://dx.doi.org/10.1073/pnas.1301725110. Kundu, P. et al., 2012. Differentiating BOLD and non-BOLD signals in fMRI time series using multi-echo EPI. Neuroimage, 60(3), pp.1759-70. Available at: http://dx.doi.org/10.1016/j.neuroimage.2011.12.028. Lombardo, M.V. et al., 2016. Improving effect size estimation and statistical power with multi-echo fMRI and its impact on understanding the neural systems supporting mentalizing. Neuroimage, Available at: http://dx.doi.org/10.1016/j.neuroimage.2016.07.022. Lütcke, H., Merboldt, K.-D. & Frahm, J., 2006. The cost of parallel imaging in functional MRI of the human brain. Magn Reson Imaging, 24(1), pp.1-5. Available at: http://dx.doi.org/10.1016/j.mri.2005.10.028. Maier, J.K., Vavrek, R.M. & Glover, G.H., 1992. Correction of NMR data acquired by an echo-planar technique. US Patent, (5,151,656). Mansfield, P., 1977. Multi-planar image formation using NMR spin echoes. J. Phys. C: Solid State Phys., 10(3), p.L55. Available at: http://dx.doi.org/10.1088/0022-3719/10/3/004. Milad, M.R. et al., 2009. Neurobiological basis of failure to recall extinction memory in posttraumatic stress disorder. Biol Psychiatry, 66(12), pp.1075-82. Available at: http://dx.doi.org/10.1016/j.biopsych.2009.06.026. Milad, M.R. et al., 2005. Thickness of ventromedial prefrontal cortex in humans is correlated with extinction memory. Proc Natl Acad Sci U S A, 102(30), pp.10706-11. Available at: http://dx.doi.org/10.1073/pnas.0502441102. Milad, M.R. & Quirk, G.J., 2012. Fear extinction as a model for translational neuroscience: ten years of progress. Annu Rev Psychol, 63, pp.129-51. Available at: http://dx.doi.org/10.1146/annurev.psych.121208.131631. Milad, M.R. & Rauch, S.L., 2007. The role of the orbitofrontal cortex in anxiety disorders. Ann N Y Acad Sci, 1121, pp.546-61. Available at: http://dx.doi.org/10.1196/annals.1401.006. Mitterschiffthaler, M.T. et al., 2006. Applications of functional magnetic resonance imaging in psychiatry. J Magn Reson Imaging, 23(6), pp.851-61. Available at: http://dx.doi.org/10.1002/jmri.20590. Moeller, S. et al., 2010. Multiband multislice GE-EPI at 7 tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI. Magn Reson Med, 63, pp.1144-53. Available at: http://dx.doi.org/10.1002/mrm.22361. 35

Norris, D.G., 2006. Principles of magnetic resonance assessment of brain function. J Magn Reson Imaging, 23(6), pp.794-807. Available at: http://dx.doi.org/10.1002/jmri.20587. Ogawa, S., Lee, T.M., Kay, A.R. & Tank, D.W., 1990. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc Natl Acad Sci U S A, 87(24), pp.9868-72. Olafsson, V. et al., 2015. Enhanced identification of BOLD-like components with multi-echo simultaneous multi-slice (MESMS) fMRI and multi-echo ICA. Neuroimage, 112, pp.43-51. Available at: http://dx.doi.org/10.1016/j.neuroimage.2015.02.052. Peters, A.M. et al., 2007. T2* measurements in human brain at 1.5, 3 and 7 T. Magn Reson Imaging, 25(6), pp.748-53. Available at: http://dx.doi.org/10.1016/j.mri.2007.02.014. Poldrack, R.A., 2012. The future of fMRI in cognitive neuroscience. Neuroimage, 62(2), pp.1216-20. Available at: http://dx.doi.org/10.1016/j.neuroimage.2011.08.007. Poser, B.A., Versluis, M.J., Hoogduin, J.M. & Norris, D.G., 2006. BOLD contrast sensitivity enhancement and artifact reduction with multiecho EPI: parallel-acquired inhomogeneity-desensitized fMRI. Magn Reson Med, 55(6), pp.1227-35. Available at: http://dx.doi.org/10.1002/mrm.20900. Posse, S., 2012. Multi-echo acquisition. Neuroimage, 62(2), pp.665-71. Available at: http://dx.doi.org/10.1016/j.neuroimage.2011.10.057. Posse, S. et al., 1999. Enhancement of BOLD-contrast sensitivity by single-shot multi-echo functional MR imaging. Magn Reson Med, 42(1), pp.87-97. Preibisch, C. et al., 2003. Functional MRI using sensitivity-encoded echo planar imaging (SENSE-EPI). Neuroimage, 19(2 Pt 1), pp.412-21. Preibisch, C. et al., 2008. Comparison of parallel acquisition techniques generalized autocalibrating partially parallel acquisitions (GRAPPA) and modified sensitivity encoding (mSENSE) in functional MRI (fMRI) at 3T. J Magn Reson Imaging, 27(3), pp.590-98. Available at: http://dx.doi.org/10.1002/jmri.21191. Price, J.L. & Drevets, W.C., 2010. Neurocircuitry of mood disorders. Neuropsychopharmacology, 35(1), pp.192-216. Available at: http://dx.doi.org/10.1038/npp.2009.104. Pruessmann, K.P., Weiger, M., Scheidegger, M.B. & Boesiger, P., 1999. SENSE: sensitivity encoding for fast MRI. Magn Reson Med, 42(5), pp.952-62. Sabatinelli, D. et al., 2011. Emotional perception: meta-analyses of face and natural scene processing. Neuroimage, 54(3), pp.2524-33. Available at: http://dx.doi.org/10.1016/j.neuroimage.2010.10.011. Schmidt, C.F. et al., 2005. Sensitivity-encoded (SENSE) echo planar fMRI at 3T in the medial temporal lobe. Neuroimage, 25(2), pp.625-41. Available at: http://dx.doi.org/10.1016/j.neuroimage.2004.12.002. Schmiedeskamp, H. et al., 2010. Improvements in parallel imaging accelerated functional MRI using multiecho echo-planar imaging. Magn Reson Med, 63(4), pp.959-69. Available at: http://dx.doi.org/10.1002/mrm.22222. Schwarzbauer, C. et al., 2010. Dual echo EPI--the method of choice for fMRI in the presence of magnetic field inhomogeneities? Neuroimage, 49(1), pp.316-26. Available at: http://dx.doi.org/10.1016/j.neuroimage.2009.08.032.

36

Setsompop, K. et al., 2012. Blipped-controlled aliasing in parallel imaging for simultaneous multislice echo planar imaging with reduced g-factor penalty. Magn Reson Med, 67(5), pp.1210-24. Available at: http://dx.doi.org/10.1002/mrm.23097. Stenger, V.A., Boada, F.E. & Noll, D.C., 2000. Three-dimensional tailored RF pulses for the reduction of susceptibility artifacts in T(*)(2)-weighted functional MRI. Magn Reson Med, 44(4), pp.525-31. Todd, N. et al., 2016. Evaluation of 2D multiband EPI imaging for high-resolution, whole-brain, task-based fMRI studies at 3T: Sensitivity and slice leakage artifacts. NeuroImage, 124(Pt A), pp.32-42. Tzourio-Mazoyer, N. et al., 2002. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage, 15(1), pp.273-89. Available at: http://dx.doi.org/10.1006/nimg.2001.0978. Vlieger, E.-J., Majoie, C.B., Leenstra, S. & Den Heeten, G.J., 2004. Functional magnetic resonance imaging for neurosurgical planning in neurooncology. Eur Radiol, 14(7), pp.1143-53. Available at: http://dx.doi.org/10.1007/s00330-004-2328-y. Wastling, S.J. & Barker, G.J., 2015. Designing hyperbolic secant excitation pulses to reduce signal dropout in gradient-echo echo-planar imaging. Magn Reson Med, 74(3), pp.661-72. Available at: http://dx.doi.org/10.1002/mrm.25444. Weiskopf, N., Hutton, C., Josephs, O. & Deichmann, R., 2006. Optimal EPI parameters for reduction of susceptibility-induced BOLD sensitivity losses: a whole-brain analysis at 3 T and 1.5 T. Neuroimage, 33(2), pp.493-504. Available at: http://dx.doi.org/10.1016/j.neuroimage.2006.07.029. Weiskopf, N. et al., 2007. Optimized EPI for fMRI studies of the orbitofrontal cortex: compensation of susceptibility-induced gradients in the readout direction. MAGMA, 20(1), pp.39-49. Available at: http://dx.doi.org/10.1007/s10334-006-0067-6. Weiskopf, N., Klose, U., Birbaumer, N. & Mathiak, K., 2005. Single-shot compensation of image distortions and BOLD contrast optimization using multi-echo EPI for real-time fMRI. Neuroimage, 24(4), pp.1068-79. Available at: http://dx.doi.org/10.1016/j.neuroimage.2004.10.012. Xu, J. et al., 2013. Evaluation of slice accelerations using multiband echo planar imaging at 3 T. Neuroimage, 83, pp.991-1001. Available at: http://dx.doi.org/10.1016/j.neuroimage.2013.07.055. Yang, C. et al., 2010. Four-dimensional spectral-spatial RF pulses for simultaneous correction of B1+ inhomogeneity and susceptibility artifacts in T2*-weighted MRI. Magn Reson Med, 64(1), pp.1-8. Available at: http://dx.doi.org/10.1002/mrm.22471. Yang, C., Deng, W. & Stenger, V.A., 2011. Simple analytical dual-band spectral-spatial RF pulses for B(1) + and susceptibility artifact reduction in gradient echo MRI. Magn Reson Med, 65(2), pp.370-76. Available at: http://dx.doi.org/10.1002/mrm.22725. Yang, C., Poser, B.A., Deng, W. & Stenger, V.A., 2012. Spectral decomposition of susceptibility artifacts for spectral-spatial radiofrequency pulse design. Magn Reson Med, 68(6), pp.1905-10. Available at: http://dx.doi.org/10.1002/mrm.24208. Yip, C.-Y., Fessler, J.A. & Noll, D.C., 2006. Advanced three-dimensional tailored RF pulse for signal recovery in T2*-weighted functional magnetic resonance imaging. Magn Reson Med, 56(5), pp.1050-59. Available at: http://dx.doi.org/10.1002/mrm.21048. 37

Yip, C.-Y. et al., 2009. Spectral-spatial pulse design for through-plane phase precompensatory slice selection in T2*-weighted functional MRI. Magn Reson Med, 61(5), pp.1137-47. Available at: http://dx.doi.org/10.1002/mrm.21938. Zald, D.H. et al., 2014. Meta-analytic connectivity modeling reveals differential functional connectivity of the medial and lateral orbitofrontal cortex. Cereb Cortex, 24(1), pp.232-48. Available at: http://dx.doi.org/10.1093/cercor/bhs308. Zhang, W.-N. et al., 2013. The neural correlates of reward-related processing in major depressive disorder: a meta-analysis of functional magnetic resonance imaging studies. J Affect Disord, 151(2), pp.531-39. Available at: http://dx.doi.org/10.1016/j.jad.2013.06.039. Zhan, X. & Yu, R., 2015. A Window into the Brain: Advances in Psychiatric fMRI. Biomed Res Int, 2015, p.542467. Available at: http://dx.doi.org/10.1155/2015/542467.

38