1. Neuroscience
Download icon

Threat of shock increases excitability and connectivity of the intraparietal sulcus

  1. Nicholas L Balderston  Is a corresponding author
  2. Elizabeth Hale
  3. Abigail Hsiung
  4. Salvatore Torrisi
  5. Tom Holroyd
  6. Frederick W Carver
  7. Richard Coppola
  8. Monique Ernst
  9. Christian Grillon
  1. National Institute of Mental Health, National Institutes of Health, United States
Research Article
  • Cited 11
  • Views 1,377
  • Annotations
Cite this article as: eLife 2017;6:e23608 doi: 10.7554/eLife.23608

Abstract

Anxiety disorders affect approximately 1 in 5 (18%) Americans within a given 1 year period, placing a substantial burden on the national health care system. Therefore, there is a critical need to understand the neural mechanisms mediating anxiety symptoms. We used unbiased, multimodal, data-driven, whole-brain measures of neural activity (magnetoencephalography) and connectivity (fMRI) to identify the regions of the brain that contribute most prominently to sustained anxiety. We report that a single brain region, the intraparietal sulcus (IPS), shows both elevated neural activity and global brain connectivity during threat. The IPS plays a key role in attention orienting and may contribute to the hypervigilance that is a common symptom of pathological anxiety. Hyperactivation of this region during elevated state anxiety may account for the paradoxical facilitation of performance on tasks that require an external focus of attention, and impairment of performance on tasks that require an internal focus of attention.

https://doi.org/10.7554/eLife.23608.001

eLife digest

Anxiety disorders affect around one in five Americans, and in many cases people experience anxiety so intensely that they have difficulties performing day-to-day activities. To help these people, it is important to understand how anxiety works. Current research suggests that anxiety disorders are caused when the connections in the brain that control our response to threat are either excessively or inappropriately activated. However, it was not clear what causes the anxiety to last for long periods. To better understand this phenomenon, Balderston et al. studied the brains of over 30 volunteers using two types of measurements called magnetoencephalography and fMRI.

In the each experiment, participants experienced periods of threat, where they could receive unpredictable electric shocks. In the first experiment, Balderston et al. measured the brain activity by recording the magnetic fields generated in the brain. In the second experiment, they used fMRI to record changes in the blood flow throughout the brain to measure how the different regions in the brain communicate.

The recordings identified a single part of the brain that increased its activity and changed its communication pattern with the other regions in the brain, when people are anxious. This region in a part of the brain called parietal lobe, is also important for processing attention, which suggests that anxiety might make people also more aware of their surroundings. However, this extra awareness might also make it more difficult for people to concentrate.

Future studies may be able to stimulate this area of the brain through the scalp to potentially reduce anxiety, as the affected area is close to the skull.

https://doi.org/10.7554/eLife.23608.002

Introduction

Current models of anxiety disorders suggest that pathological anxiety results from excessive or inappropriate activation of the same neural circuits that are responsible for adaptive anxiety in the face of threat (Insel et al., 2010; Insel, 2014). Although there is a long history of translational work studying neural systems mediating the acute fear response (Pavlov, 1927; Fanselow and Poulos, 2005; Fullana et al., 2016), much less is known about the neural systems mediating prolonged periods of elevated state anxiety. Closing this knowledge gap is critical because the occurrence of prolonged periods of elevated state anxiety is one of the primary symptoms of all anxiety disorders (American Psychiatric Association, 2013). Therefore, understanding neural mechanisms underlying prolonged periods of elevated anxiety has the potential to identify targets for the treatment of anxiety disorders, which are among the most prevalent psychiatric disorders (Kessler et al., 2005).

The gold-standard translational paradigm for studying elevated state anxiety in the laboratory is the threat of shock paradigm, during which subjects are exposed to periods when they are either safe, or at risk for receiving unpredictable aversive electrical stimulations (Schmitz and Grillon, 2012; Grillon, 2008; Grillon and Baas, 2003). It allows for the experimental manipulation of state anxiety within subjects (Grillon and Baas, 2003, 1998; Grillon et al., 1991, 2007, 2008, 2009), which can be quantified using psychological and physiological measures (Grillon, 2008; Grillon and Baas, 2003) and can be implemented in healthy subjects (Balderston et al., 2017a, 2017b; Cornwell et al., 2007, Cornwell et al., 2008, 2012; Lissek et al., 2007), patients (Grillon et al., 2009; Balderston et al., 2017a; Vytal et al., 2016), and non-human animals (Davis et al., 2010). A key feature of prolonged periods of threat of shock is that they induce a stable increase in anxiety that can be probed at random intervals using the acoustic startle reflex (Grillon, 2008; Grillon and Baas, 2003, 1998), suggesting that this anxious state is mediated by a fundamental sustained change in the pattern of ongoing brain activity. This is in contrast to more phasic event-related fear responses in typical cued fear conditioning (Schmitz and Grillon, 2012; Grillon, 2008; Davis et al., 2010).

Current neuroscientific models of anxiety are based in part on translational work using Pavlovian fear conditioning (Pavlov, 1927; Fanselow and Poulos, 2005; Fullana et al., 2016). Decades of work in non-human animals has shown that acute fear responding is dependent upon the amygdala (Kwapis et al., 2009; Bailey et al., 1999; Parsons et al., 2006), and functional magnetic resonance imaging (fMRI) during fear conditioning in humans has been used to identify a canonical fear network that includes the amygdala, the dorsomedial prefrontal cortex, the thalamus, and the anterior insula (Fullana et al., 2016; Schultz et al., 2012; Cheng et al., 2003, 2006). However, much less is known about the network mediating extended periods of elevated state anxiety.

In addition, cognitive scientific research in humans shows that attentional processing is profoundly influenced by both state (Vytal et al., 2013, 2012; Patel et al., 2016; Shackman et al., 2006) and trait anxiety (Derakshan et al., 2009; Eysenck et al., 2007), suggesting that multiple neural systems are affected by anxiety. Although there have been some studies investigating the neural systems that mediate anxiety, these studies often depend on an a priori focus that is centered on the regions of the canonical fear network, and typically rely on a priori methods to increase statistical sensitivity in these regions such as lowered statistical thresholds (Robinson et al., 2013a; Mobbs et al., 2010; Hooker et al., 2006; Tabbert et al., 2010), region of interest analyses (Balderston et al., 2015, 2014, 2013), and seed-based functional connectivity (Schultz et al., 2012; Vytal et al., 2014; Gold et al., 2015). Importantly, the increased sensitivity gained by using these statistical methods comes at the cost of assessing anxiety-related changes in regions not identified a priori, thus resulting in a possible under-reporting of anxiety-related changes in other areas of the brain, such as regions important for attentional processing. Therefore, the purpose of this study was to use exploratory analytical methods to identify the most prominent activity/connectivity changes induced by the threat of shock paradigm.

Toward this aim, we collected data from two complimentary imaging modalities, fMRI and magnetoencephalography (MEG) during a threat of shock paradigm. In both MEG and fMRI experiments, subjects underwent alternating blocks of safety and threat, and rated their anxiety continuously using a centrally located visual analog scale (see Figure 1). During the MEG experiment, subjects also received randomly timed white noise presentations, which served to probe the subject’s current anxiety level (via the acoustic startle response) and their ongoing brain activity (via the preceding pattern of neural oscillations).

Schematic of experimental paradigm.

(A) Subjects underwent alternating blocks of threat and safety. (B) Visual display present on the screen during the experiment. During the experiment subjects saw two circles. The color of the outer circle indicated the block type. The color of the inner circle was controlled by the subject, and reflected the subject’s then-current anxiety level.

https://doi.org/10.7554/eLife.23608.003

In both experiments, we used unbiased, data-driven, whole-brain approaches to identify changes in activity (MEG) and connectivity (fMRI) as a function of threat. To assess functional connectivity changes in the fMRI signal, we used the global brain connectivity (GBC) metric, which does not rely on a priori seed-selection for the connectivity analysis. To assess ongoing patterns of activity in the MEG study, we evaluated changes in oscillatory power during the 2 s prior to the startle probes as a function of threat. According to the translational approach, one might predict that the most prominent changes in spontaneous neural activity and connectivity would emerge in regions of the canonical fear network (Fanselow and Poulos, 2005; Kim et al., 2011). However, given that impaired attentional control is a key feature of clinical anxiety (Derakshan et al., 2009; Eysenck et al., 2007), and that threat of shock has been repeatedly shown to impact performance on tasks that require attention control (Vytal et al., 2013, 2012; Patel et al., 2016; Shackman et al., 2006), one might also expect that the most prominent changes would emerge within regions of the frontoparietal attention network (Ptak, 2012; Posner, 2012; Petersen and Posner, 2012).

Results

Threat increases subjective and physiological measures of anxiety

We began by assessing the ability of our threat manipulation to induce a sustained state of anxiety. Results from the psychological questionnaires, and the subjective rating scales can be found in Table 1. Consistent with the online ratings, subjects during both experiments reported more anxiety (MEG: t(26) = 8.65; p<0.001; fMRI: t(24) = 13.98; p<0.001) and fear (MEG: t(26) = 8.03; p<0.001; fMRI: t(24) = 9.15; p<0.001) during the threat blocks than during the safe blocks. In addition, two sample t-tests did not reveal any significant differences between experiments for either the psychological questionnaires, or the affective rating scales (all ps > 0.05). For the MEG study, we analyzed both the acoustic startle reflex and the online self-reported anxiety ratings. Because startle probes could not be presented during the MRI study, we relied only on the ratings.

Table 1

Individual differences for MEG (N = 28) and MRI (N = 25) experiments.

https://doi.org/10.7554/eLife.23608.004
MeasureMEGMRI
STAI
 State26.04 (1.37)23 (0.9)
 Trait27.12 (0.93)28.18 (1.27)
ASI11.59 (1.21)8.64 (1.18)
BAI1.37 (0.42)0.58 (0.26)
BDI0.89 (0.32)0.42 (0.19)
Shock
 Intensity (mA)5.66 (0.66)6.91 (1.01)
 Rating8.51 (0.2)9.09 (0.19)
Anxiety
 Pre2.04 (0.27)1.98 (0.25)
 Safe2.47 (0.31)1.76 (0.21)
 Threat5.41 (0.37)5.97 (0.39)
Fear
 Pre1.41 (0.15)1.5 (0.23)
 Safe1.84 (0.27)1.27 (0.12)
 Threat4.44 (0.39)4.7 (0.42)
  1. Note: Numbers reflect the mean and standard deviation of the results [M (SD)].

For each startle probe in the MEG study, we extracted the subject’s startle magnitude, and anxiety rating recorded just prior to the startle probe. Both the startle magnitude and anxiety ratings were normalized and converted to T-scores (Blumenthal et al., 2005) within subjects. These values were then averaged across trials and submitted to a paired-sample t-test (Safe vs. Threat). Both ratings (See Figure 2A and Figure 2—source data 1; t(27) = 10.03; p<0.001) and startle (See Figure 2B and Figure 2—source data 1; t(27) = 4.65; p<0.001) indicated greater anxiety during the threat blocks compared to the safe blocks. Next, we created Threat > Safe difference scores for both startle (anxiety potentiated startle; APS) and the online ratings, and correlated the values across subjects. Ratings within the MEG study were significantly correlated with startle across subjects (See Figure 2D and Figure 2—source data 1; r(26) = 0.61; p=0.001).

Behavioral results from both experiments.

(A) Anxiety ratings during the MEG study. (B) Startle magnitude during the MEG study. (C) Anxiety ratings during the fMRI study. Bars represent the mean ± within-subject SEM (Cousineau, 2005). (D) Correlations between anxiety potentiated startle (APS) and differential anxiety ratings. The black squares represent the correlation between APS and ratings during the MEG session. The red dots represent the correlation between APS during the MEG study and anxiety ratings during the fMRI study in the subset of subjects who participated in both studies.

https://doi.org/10.7554/eLife.23608.005
Figure 2—source data 1

Source data for all graphs in Figure 2.

https://doi.org/10.7554/eLife.23608.006

During the MRI study, we averaged the ratings across time in the threat and safe blocks, and converted these values to T-scores. As with the MEG study, these values were submitted to a paired-sample t-test (Safe vs. Threat), and indicated more anxiety during the threat blocks than the safe blocks (See Figure 2C and Figure 2—source data 1; t(24) = 23.06; p<0.001). We also created Threat > Safe difference scores for these values, and correlated these difference scores with startle (recorded during the MEG study) in the subjects who participated in both sessions. There was a nonsignificant small positive correlation between startle and rating during the MRI session (see Figure 2D and Figure 2—source data 1; r(16) = 0.24; p=0.344).

Threat increases whole brain GBC

Many threat of shock studies have used seed-based functional connectivity analyses to identify changes in emotional processing centers in the brain (Vytal et al., 2014; Gold et al., 2015; Satterthwaite et al., 2011; Prater et al., 2013; Hrybouski et al., 2016; Birn et al., 2014; Cha et al., 2014; Heitmann et al., 2016); however, seed-based functional connectivity methods suffer from bias because they require the experimenter to select a seed region ahead of time, while ignoring all other possible seed regions. To address this limitation, researchers have developed complementary data-driven functional connectivity metrics such as GBC, which do not rely on a priori seed-selection for the connectivity analysis. By assessing the connectivity between each voxel and every other voxel, this analysis allows the user to identify the most connected regions of the brain (Cole et al., 2010), as well as the seed regions where connectivity impacts behavior across subjects (Cole et al., 2012; Gotts et al., 2012). By identifying regions that show the largest changes in GBC during periods of threat vs. periods of safety, it is possible to identify hubs that contribute most prominently to the sustained anxious state during threat of shock.

We collected whole brain multi-echo echo-planar imaging data, and used the echo time-dependent independent components analysis to remove sources of noise unrelated to the blood oxygenation level dependent (BOLD) response from the timeseries (Kundu et al., 2012; Evans et al., 2015). Subjects were exposed to alternating 2-min blocks of safety and threat, without startle probes. Given that this design lacked external timing information (i.e. external stimulus presentations), we examined changes in functional connectivity as a function of block type. We opted for a data-driven GBC approach where the connectivity of every voxel was assessed. We first computed GBC maps independently for the safe and threat conditions by correlating each voxel’s timecourse with every other voxel’s timecourse, applying the Fisher’s Z transformation, and averaging across these correlation maps (See Figure 3A,B and Figure 3—source data 1). As a first pass, we averaged across all voxels to obtain the whole brain GBC for safe and threat. We then conducted a (Safe vs. Threat) paired-sample t-test on these values, and found significantly more GBC for threat blocks than for safe blocks (see Figure 3C and Figure 3—source data 1; t(24) = 2.13; p=0.044).

Overview of global brain connectivity (GBC) measure.

(A) Map showing average GBC across all safe and threat TRs. (B) Cartoon schematic of a correlation matrix. The 43204 voxel x 43204 voxel cross correlation matrix was calculated separately for each subject and each condition. Correlations were averaged across rows for the entire grey matter mask, to create a single map reflecting the average correlation between each voxel and all other voxels in the mask. (C) Graph representing the mean GBC following the Fisher’s Z transformation for safe and threat averaged across the entire grey matter mask. Bars represent the mean ± within-subject SEM (Cousineau, 2005).

https://doi.org/10.7554/eLife.23608.007

Threat increases voxelwise GBC in the intraparietal sulcus (IPS)

To follow-up the whole brain analysis, we conducted a voxelwise analysis of GBC. Using the same GBC maps created above, we conducted a voxelwise (Safe vs. Threat) paired-sample t-test. We used Monte Carlo simulations to estimate a null distribution for statistical testing, and used a cluster-based method based on this null distribution to correct for multiple comparisons. We found a significant increase in GBC in the threat blocks compared to the safe blocks in three clusters (see Table 2). The largest cluster was in the right angular gyrus. The two remaining clusters were found bilaterally in the IPS (see Figure 4A and Figure 4—source data 1). In all three clusters, we found significantly higher GBC for threat blocks than for safe blocks (see Figure 4B and Figure 4—source data 1). To determine whether these differences were affected by the delivery of the shock, or differences in motion across blocks, we repeated the threat vs. safe analysis at the cluster level after censoring the 10 TRs following shock delivery, and an equivalent number of safe TRs closely matched for motion. In addition, we covaried out any remaining differences in motion using an analysis of covariance. Using this approach, we still found a significant effect of threat on GBC in all three regions (Right angular gyrus: f(1, 23)=4.71, p=0.04; Right IPS: f(1, 23)=7.22, p=0.01; Left IPS: f(1, 23)=5.39, p=0.03), suggesting that our initial findings were not due to differences in motion, censoring, or residual neural activity evoked by the shock.

Table 2

Results from voxelwise GBC analysis.

https://doi.org/10.7554/eLife.23608.009
LabelVolumet-valuePeak activation (LPI)
xyz
Right Angular Gyrus1583.4548−5127
Right Intraparietal Sulcus833.4221−6066
Left Intraparietal Sulcus813.6−18−6366
Results from voxelwise global brain connectivity (GBC) analysis.

(A) Statistical map showing results from a threat vs. safe paired-sample t-test. (B) Graph representing average GBC values after applying the Fisher’s Z transformation for clusters shown in panel A. Bars represent the mean ± within-subject SEM (Cousineau, 2005).

https://doi.org/10.7554/eLife.23608.010

Threat increases connectivity within an IPS-centered attentional network

Given that the IPS emerged as a hub differentiating global connectivity in threat compared to safe, we used this region as a seed-region to identify changes in functional connectivity during the threat vs. safe blocks. We understand this follow-up analysis could be interpreted as circular. According to this perspective, the seed-based analysis is not independent from the GBC analysis, which was used to identify the seed. However, the purpose of the follow-up analysis (i.e. probing functional connectivity using clusters identified from the global connectivity analysis) is to limit the interpretations of the global connectivity results to those supported by seed-based functional connectivity results, which has been done in previous group-level GBC studies (e.g. [Gotts et al., 2012]).

We extracted the timecourse of activity averaged across the voxels in the bilateral IPS functional ROIs, and correlated this timecourse with the timecourse of activity across all voxels in the brain, independently for safe and threat, and applied the Fisher’s Z transformation to the resulting correlation coefficients. We conducted a voxelwise (Safe vs. Threat) paired-sample t-test on the resulting maps. We used Monte Carlo simulations and cluster thresholding to correct for multiple comparisons. Consistent with the GBC results, we found an increase in connectivity in threat blocks compared to safe blocks, in several regions of the frontoparietal attention network (See Table 3, Figure 5, and Figure 5—source data 1).

Table 3

Results from voxelwise IPS connectivity analysis.

https://doi.org/10.7554/eLife.23608.012
LabelVolumet-valuePeak activation (LPI)
xyz
Left Thalamus3423.92-9612
Right Inferior Parietal Lobule2083.6757−5739
Left Superior Medial Gyrus1843.6533642
Left Precuneus1793.593−6948
Right Middle Frontal Gyrus1373.64331560
Left Angular Gyrus1133.51−57−5430
Left Middle Frontal Gyrus963.69−241560
Left Middle Frontal Gyrus903.48−4551-3
Results from bilateral IPS seed-based connectivity analysis.

(A) Statistical map showing results from a threat vs. safe paired-sample t-test. (B) Graph representing average IPS connectivity values for clusters shown in panel A. Bars represent the mean ± within-subject SEM (Cousineau, 2005).

https://doi.org/10.7554/eLife.23608.013

Alpha oscillations dominate neuromagnetic recordings at rest

Next, we characterized the pattern of activity in the MEG study. It is well established that spontaneous neural activity at rest is dominated by oscillations in the alpha (8–12 Hz) range (de Munck et al., 2008; Sadaghiani et al., 2010; Mo et al., 2013; Mayhew et al., 2013), which are most prominent when the subject is in an alert state of restful relaxation (Doufesh et al., 2014; Kim et al., 2014; Khalsa et al., 2015; Lagopoulos et al., 2009), and that alpha asymmetries can reflect differences in arousal across subjects (Nitschke et al., 1999). Theoretical models of alpha function suggest that alpha oscillations are generated by coherent activity in local inhibitory interneurons (Klimesch et al., 2007) and that decreases in alpha power reflects increases in cortical excitability (Klimesch et al., 2007; Lange et al., 2013). Consistent with these theories, studies collecting simultaneous measures of electroencephalography (EEG) and fMRI have shown that alpha power is negatively correlated with functional connectivity (Laufs et al., 2003; Scheeringa et al., 2012; Chen et al., 2008; Chang et al., 2013).

We extracted and cleaned the 2 s of data prior to each startle probe, avoiding contamination by blink artifacts. Because the timing of the startle probes was random, the pre-stimulus recording reflected random sampling across the sustained threat period. Therefore, we collapsed across the 2 s interval and examined oscillatory activity. First, we transformed the values into the frequency domain using a Fast Fourier Transform (FFT) with upper and lower limits of 20 Hz and 1 Hz, respectively. Then, we averaged these values across sensors, trials, and subjects, and examined the spectrogram for peaks. We detected a peak at ~10 Hz, suggesting that alpha oscillations were a prominent feature of these recordings (Figure 6A and Figure 6—source data 1). We identified the largest local maxima in each subject’s spectrogram. In all but four subjects we detected a peak in the alpha frequency band (8 Hz – 12 Hz). The power within a 2 Hz band around these individual alpha frequencies (IAF)s was used in all subsequent analyses (Figure 6B and Figure 6—source data 1). For subjects without a detectable IAF, power in a narrow band around IAF averaged across subjects was used. Subsequent analyses were performed in both sensor space (Figure 6C; black dots) and source space (Figure 6C; green dots).

Overview of MEG analyses.

(A) Spectrogram representing power averaged across all subjects and all sensors with peak in the alpha frequency band. (B) Graph showing the frequency of peak alpha (individual alpha frequency) averaged across subjects. Bars represent the mean ± SEM. (C) Example of single subject alignment with sensors (black dots) source grid (green dots) and headmodel (surface) plotted together.

https://doi.org/10.7554/eLife.23608.015

Threat reduces parietal alpha oscillations

Given that alpha oscillations were the dominant feature in the MEG recordings across all blocks, we examined whether these oscillations differed as a function of threat. For the sensor space analysis, we averaged IAF across trials within conditions, and then performed a paired-sample t-test (Safe vs. Threat) on the resulting sensor space averages. We used Monte Carlo simulations and cluster thresholding to correct for multiple comparisons. Importantly, we found a large cluster of sensors over the left parietal lobe and a smaller cluster of frontal sensors showing significantly less IAF power in the threat blocks than during the safe blocks (Figure 7A–B and Figure 7—source data 1).

Alpha results from threat vs. safe t-test.

(A) Statistical map in sensor space showing a significant reduction in alpha power. Black symbols represent clusters of sensors showing significant threat vs. safe differences. (B) Graph showing average alpha power for safe and threat conditions in the largest cluster of sensors in panel A. (C) Statistical map in source space showing a significant reduction in alpha power. (D) Graph showing average alpha power for safe and threat conditions in the cluster in panel C. Bars represent the mean ± within-subject SEM (Cousineau, 2005).

https://doi.org/10.7554/eLife.23608.017

To localize the source of these power differences, we projected these signals into source space using a dynamic imaging of coherent sources (DICS) beamformer. First, we created a common filter using all trials, then we projected the safe and threat trials through the filter independently, to obtain power estimates at the source level for each condition. Finally, we conducted a (Safe vs. Threat) paired-sample t-test on the source space IAF power estimates. As with the sensor space data, we used Monte Carlo simulations and cluster thresholding to correct for multiple comparisons. Consistent with the sensor space results, we found two clusters of voxels showing significantly less IAF power in the threat and safe conditions (Figure 7C–D and Figure 7—source data 1), the larger cluster was located in the left IPS, while the smaller cluster was located in the mid-cingulate gyrus. In addition, both the sensor and the source space results held if we analyze power across the entire alpha frequency band (8 Hz to 12 Hz; See below). Finally, when comparing the MEG results with the fMRI results, we found that the left IPS cluster in the fMRI data substantially overlaps with the alpha cluster (46/81 voxels; See Figure 8).

Conjunction map from voxelwise fMRI GBC analysis and MEG alpha power differences.

Colors represent significant safe vs. threat differences from the fMRI analysis (yellow), MEG analysis (blue), and both analyses (green).

https://doi.org/10.7554/eLife.23608.019

As with the MRI data, it is important to show that our denoising steps did not affect the comparisons between the safe and threat conditions. Specifically, it is important to show that the number of trials rejected either (1) did not differ across the safe and threat blocks, or (2) did not affect the safe vs. threat comparisons. In the safe condition, there were on average 58.21 ± 5.5 trials, while in the threat condition there were on average 55.32 ± 7.6 trials, which was significantly different across subjects (t(27) = 2.78, p=0.01). Therefore, we decided to determine whether this difference in trial number impacted our results at the censor and source level. Accordingly, we included the difference in trial number across safe and threat blocks as a covariate in an ANCOVA examining the effect of threat on alpha responses. At both the sensor (f(1,26) = 7.48, p=0.01) and at the source (f(1,26) = 17.797, p<0.001), we find that even with the difference in number of trials covaried out, the effects of threat on alpha power is still significant, suggesting that this difference did not significantly impact our results.

Discussion

The purpose of this study was to use the complementary methods of fMRI and MEG to identify the regions and network hubs that contribute most prominently to sustained anxiety during threat of shock. In both cases, we used unbiased, data-driven, whole-brain measures of activity/connectivity, and identified the IPS as the region most affected by the threat of shock manipulation. Using fMRI GBC, we found that the IPS increased its connectivity during periods of threat with a set of regions important for attention control, suggesting that this region is a hub in a network important for attentional processing during threat (Du et al., 2012; Kincade et al., 2005; Rushworth et al., 1997; Bucci, 2009; Corbetta and Shulman, 2002). Using MEG, we found that the magnitude of spontaneous alpha oscillations in the same region decreased during periods of threat, suggesting that threat increases IPS cortical excitability (Klimesch et al., 2007). Together, these results suggest that threat of shock facilitates attentional processing mediated by increased excitability and connectivity of the IPS (Du et al., 2012; Kincade et al., 2005; Rushworth et al., 1997; Bucci, 2009). This enhanced attentional processing may reflect a state of hypervigilance, consistent with the findings that elevated state anxiety biases the attentional system (Waters et al., 2014; Amir et al., 2009; Lapointe et al., 2013). In addition, this hypervigilance may provide a mechanistic explanation for why threat facilitates performance on tasks that require an external focus of attention but impairs performance on tasks that require an internal focus of attention (Vytal et al., 2016, 2013; Patel et al., 2016; Torrisi et al., 2016; Grillon et al., 2016; Robinson et al., 2011; Balderston et al., 2016).

Functional connectivity

Previous studies exploring the relationship between threat and functional connectivity focused on how threat affected connectivity within networks centered on seeds placed in emotional processing regions (Vytal et al., 2014; McMenamin et al., 2014; McMenamin and Pessoa, 2015). However, in this study because we were specifically interested in identifying the regions of the brain that contributed most prominently to anxiety, we chose to forego an a priori seed-based approach, and employ data-driven connectivity measures that do not rely on choosing a seed (Cole et al., 2010; Calhoun et al., 2009). To start, we used GBC as a method to determine whether connectivity en masse changed with threat, and to identify where in the brain the largest changes occurred (Cole et al., 2010, 2012; Chu et al., 2011). We found the largest changes to occur bilaterally in the IPS, where connectivity with the rest of the brain increased as a function of threat. Then, using these regions as a seed, we found that they supported an increase in connectivity within a set of regions important for executive control (Posner, 2012; Corbetta and Shulman, 2002). These results suggest that the IPS is a critical connectivity hub in the network of regions that contribute most prominently to the sustained elevation of state anxiety induced by the threat of shock paradigm (Cole et al., 2010, Cole et al., 2014). Although our GBC approach did not reveal any hubs within the amygdala, or in other regions typically associated with emotional processing (Simmons et al., 2006; Mechias et al., 2010; Etkin and Wager, 2007), we do see an increase in connectivity between the IPS and the thalamus and dorsomedial prefrontal cortex, two regions known to be part of the canonical fear network.

In addition to the IPS, we found changes in GBC as a function of threat in the right angular gyrus. Although it is currently unclear how the angular gyrus might contribute to threat processing, there is work suggesting that this region may be a key site for multisensory integration (Seghier, 2013). One hypothesis is that enhanced GBC in this region may reflect a heightened awareness of the current context. However, this is a post hoc explanation that should be explored in future studies.

An obvious parallel to the fMRI connectivity analysis would be to conduct a similar whole brain connectivity analysis of the MEG data. Not only would corroborating evidence from an independent imaging modality strengthen the fMRI connectivity results, but results from MEG specifically would allow for a frequency-specific analysis of the effects of threat on functional connectivity (Hillebrand et al., 2012; Brookes et al., 2011). However, the current study was not optimized to reliably detect differences in MEG connectivity. In the current study, we included the white noise probes as a way to obtain a quantitative measure of anxiety throughout the safe and threat blocks (Grillon et al., 1999). These white noise probes trigger the acoustic startle reflex, which varies as a function of an individual’s current level of anxiety (Grillon, 2008). Unfortunately, these reflexive blinks also inject an artifact into the MEG signal, and because the magnitude of these blinks differs across conditions, the blink artifact also differs across conditions. Although adaptive beamforming can theoretically remove the artifact induced by the blink response (Van Veen et al., 1997), the only way to ensure that the blink artifact does not differentially affect the MEG signal is to remove the contaminated time periods from the analysis, or remove the startle probes from the design at the outset. In the current study, we chose to address this limitation by extracting two-second time windows prior to each startle presentation to minimize the effect of the blink artifact on our power estimates. However, it has been shown that reliability of the MEG connectivity estimates increases as the duration of the recording increases, and durations of ~10 min or greater may be needed to maximize reliability (Liuzzi et al., 2016). Therefore, using such short intervals did not allow for the ability to obtain reliable estimates of MEG connectivity. Future studies should be conducted to address this limitation. In addition, it will be important to use appropriate correction methods to account for signal leakage between source space regions (Colclough et al., 2015), and verify that the resulting connectivity estimates match previously published work (i.e. strong alpha connectivity in occipital regions and strong beta connectivity linking the parietal cortex with other frontal, temporal, and occipital regions; [Hunt et al., 2016]).

fMRI and alpha

In addition to the increases in functional connectivity measured with fMRI, we also found decreases in alpha as a function of threat of shock in the same parietal region using MEG. Although the focus of the paper was on alpha, our initial approach (described in the introduction) was to examine all frequency bands independently (other bands not shown). The focus on alpha emerged out of the observations that, (1) alpha was the strongest signal in the recordings, (2) alpha showed the largest power changes as a function of threat, (3) alpha was the only frequency band that showed consistent results at both the sensor and the source level, and (4) the source space results aligned nicely with the corresponding fMRI GBC data.

These results are also consistent with several previous studies using simultaneous recordings of fMRI and EEG (Sadaghiani et al., 2010; Mo et al., 2013; Mayhew et al., 2013; Scheeringa et al., 2012; Chang et al., 2013; Wu et al., 2010; Walz et al., 2015; Scheeringa et al., 2011). For instance, functional connectivity at the whole brain level (Chang et al., 2013), and within the visual system (Scheeringa et al., 2012) is negatively correlated with posterior alpha power. In addition, alpha power is negatively correlated with activity with dorsal attention network (Sadaghiani et al., 2010), and positively correlated with activity in the default mode network (Mo et al., 2013). Finally, intrinsic connectivity within networks is negatively correlated with alpha during eyes open vs. eyes closed resting state studies (Wu et al., 2010), and both the phase and power of pre-stimulus alpha affects event-related BOLD responses in sensory regions (Mayhew et al., 2013; Walz et al., 2015; Scheeringa et al., 2011). Together these results suggest that the increase in IPS connectivity and decrease in IPS alpha power may reflect a common process.

Alpha oscillations

We also show that threat reduces the power of resting state alpha oscillations. These oscillations are thought to reflect cortical inhibition, driven by inhibitory interneurons, and triggered by top-down modulatory control (Klimesch et al., 2007). According to this view, our current findings may be driven by a release from top-down inhibition, resulting in a net increase in cortical excitability. Importantly, reductions in alpha power (i.e. increases in excitability) are associated with enhanced sensory and motor processing (Cornwell et al., 2007; Baas et al., 2006). For instance, reduced pre-stimulus alpha is associated with enhanced visual and illusory visual perception (Lange et al., 2013), and increased pre-stimulus alpha is associated with reduced transcranial magnetic stimulation-induced motor-evoked potentials (Klimesch et al., 2007). Accordingly, alpha reduction (i.e. cortical excitation) may provide a mechanistic explanation for one of the most commonly reported symptoms of elevated state anxiety, namely anxiety-potentiated startle (Schmitz and Grillon, 2012; Grillon and Ameli, 1998; Grillon et al., 1991). That is, threat of shock potentiates the startle reflex by increasing cortical excitability.

Cognitively, alpha oscillations are thought to play a key role in selective attentional processes (Klimesch, 2012), such that increases in alpha typically reflect inhibition or filtering of information that is to-be-ignored (Bonnefond and Jensen, 2013; Kelly et al., 2006; Händel et al., 2011). For instance, alpha oscillations have been shown to filter out noise during distracting listening conditions (Strauß et al., 2014). Similarly, increases in pre-stimulus alpha to predictable stimuli are lateralized to to-be-ignored locations (Horschig et al., 2014; Rihs et al., 2009). In addition, alpha power increases during the maintenance of items in working memory (WM) (Klimesch et al., 2007; Meyer et al., 2013), and this increase in alpha serves to protect these items from distractors (Bonnefond and Jensen, 2012; Manza et al., 2014). Taken together, these results suggest that the decreases in alpha power observed in our study may reflect greater perceptual sensitivity, consistent with the idea that threat may increase vigilance (Eilam et al., 2011). In addition, alpha-gamma coupling is thought to be important for allocating attention to unattended salient stimuli (Jensen et al., 2012).

The reductions in alpha power observed here were lateralized to the left hemisphere, leading to a right dominant parietal asymmetry. According to prominent models, left dominant frontal alpha is associated with positive affect and/or approach behavior, while right dominant alpha is associated with negative affect and/or avoidance behavior (Davidson, 2004; Harmon-Jones et al., 2010). Like the valence model of alpha asymmetry, one prominent model on parietal alpha asymmetry is rooted in the arousal-valence model of emotional processing (Heller, 1993). According to this model, right dominant parietal alpha is associated with increased arousal. Consistent with this theory, we find that threat of shock, which increases arousal, also reduces left parietal alpha, resulting in a right dominant profile. However, more research should be conducted to specifically test this hypothesis.

Cognition/anxiety interactions

Although our results only provide indirect evidence for the hypothesis that threat facilitates orienting, a bias toward hyper-orienting during periods of elevated state anxiety could explain many of the conflicting findings related to the cognition/anxiety interaction. The relationship between cognition and anxiety has been extensively studied (For reviews, see [Eysenck et al., 2007; Robinson et al., 2013b]). Importantly, threat of shock improves performance on sustained attention tasks (Torrisi et al., 2016; Grillon et al., 2016; Robinson et al., 2011), while impairing performance on working memory tasks (Vytal et al., 2016, 2013; Patel et al., 2016; Balderston et al., 2016). These tasks can be distinguished from one another based on the locus of attention required for performance. In sustained attention tasks, subjects are constantly monitoring the external environment for odd-ball stimuli, while in working memory tasks, subjects are required to maintain information internally. Our hypothesis is that threat-induced hyper-orienting improves performance on sustained attention tasks by reducing lapses in attention (Torrisi et al., 2016). In contrast, we hypothesize that threat-induced hyper-orienting impairs performance on working memory tasks by lowering the threshold for distractors to gain access to WM resources (Balderston et al., 2016). Consistent with these hypotheses, high pre-stimulus alpha is associated with better memory performance while low pre-stimulus alpha is associated with better perception performance (Klimesch et al., 2007).

Strengths and limitations

This study had a number of strengths. First, we used the gold standard threat of shock paradigm, which has been extensively applied to manipulate state anxiety within subjects (Grillon, 2008; Grillon and Baas, 2003, 1998; Grillon et al., 1991, 2007, 2008, 2009; Balderston et al., 2017a, 2017b; Cornwell et al., 2007, Cornwell et al., 2008, 2012; Lissek et al., 2007). Unlike previous studies, we also collected anxiety ratings throughout the recordings, and found that these ratings were similar across studies and correlated with anxiety-potentiated startle within session. Similarly, retrospective anxiety ratings were similar across studies, suggesting a comparable anxiety induction in both studies. We also collected data from multiple imaging modalities that both support the finding of enhanced IPS processing during periods of threat. In addition, we used the state-of-the-art multi-echo fMRI acquisition, and accompanying echo time-dependent independent components analysis to eliminate non-BOLD sources of noise from our fMRI data.

Although there were a number of strengths to our study, there were also several limitations. First, as mentioned above, our MEG study was not optimally designed to reliably detect differences in connectivity, which would have been an obvious parallel to the fMRI connectivity analysis. Second, because we had a relatively small sample size, this study was not optimally designed to study how individual variability in anxiety affected IPS activity/connectivity. Third, our initial goal was to collect both the fMRI and the MEG data in all subjects, but a number of subjects could not participate in both, making it difficult to compare the results at the single subject level. Fourth, given our interest in the relationship between alpha and connectivity, it might have been better to collect simultaneous fMRI and EEG. However, because of our interest in the startle data, it was important to present white noise probes in at least one of the experiments, and our MRI scanner does not currently have that capability. To overcome that limitation, we collected startle data in the MEG study, which we used to validate the continuous anxiety ratings, which were collected during both studies. Finally, another limitation with the MEG data is that we do not have a measure of within-block movement, making it difficult to determine whether motion differed between the safe and threat blocks; however, it should be noted that any trials contaminated by movement (i.e. muscle) artifact were removed.

Conclusions

Current translational neuroscientific models of anxiety focus on regions of the canonical fear network as anatomical hubs for anxiety. However, support for these models in humans is often based on studies that focus on these regions a priori (i.e. seed-based functional connectivity), at the expense of rigorous, unbiased, whole-brain approaches. In this work we conduct two separate experiments, using complimentary imaging modalities (fMRI and MEG) to assess the effect of elevated state anxiety on functional connectivity and cortical excitability across the entire brain. In these studies, we identify effects in a common region, the IPS, which is a key node in the frontoparietal attention network. These results suggest that threat enhances processing in this region, possibly facilitating attentional processing, leading to increased vigilance. However, these results have broader implications for future research and treatment. First, our results suggest that elevated anxiety in humans is primarily a cognitive state that is not fully captured by the current translational models. Second, because the most prominent region contributing to elevated state anxiety is cortical, it may be possible to target this region with neuromodulatory methods and reduce anxiety.

Materials and methods

Participants

Forty-two (28 female; age: M = 28.45, SD = 6.23) healthy volunteers from the Washington DC metropolitan area were recruited to participate in the current study. Of these, 32 (18 female; age: M = 27.82, SD = 5.11) participated in the magnetoencephalography (MEG) study, and 30 (19 female; age: M = 28.52, SD = 6.75) participated in the functional magnetic resonance imaging (fMRI) study. Among the subjects included in the final analysis, 18 (13 female; age: M = 27.67, SD = 4.76) participated in both. For the MRI study, two subjects were removed due to technical malfunction with the scanner/data, three subjects were removed because their anxiety ratings were two standard deviations below the mean. For the MEG study, three subjects were removed from the final analysis for sleeping, moving, or not paying attention to the task. One subject withdrew in the middle of the recording. The final sample included 38 total subjects (MEG N = 28; MRI N = 25): 18 subjects participated in both experiments, 10 participated only in the MEG experiment, and seven participated only in the MRI experiment. Although no explicit power analysis was conducted prior to the experiment, we chose a sample size of approximately 25 participants, to ensure enough power to detect a behavioral threat of shock effect based on previous studies (Schmitz and Grillon, 2012; Balderston et al., 2017b, 2016).

Following an initial telephone screen, participants visited the National Institutes of Health Clinical Center for a comprehensive screening by a clinician. Inclusion criteria for healthy volunteers were: (1) no current Axis I psychiatric disorder as assessed by SCID-I/NP (First et al., 2012), (2) no first-degree relative with a known psychotic disorder, (3) no interfering acute or chronic medical condition, (4) no brain abnormality on MRI as assessed by a licensed radiologist, (5) negative urine drug screen, and (6) right-handedness. All participants gave written informed consent approved by the National Institute of Mental Health (NIMH) Combined Neuroscience Institutional Review Board and received compensation for participating.

Stimulus presentation

Request a detailed protocol

Presentation software package (version 14.6, Neurobehavioral Systems, Berkeley, CA) was used to present the stimuli via projection systems in both the MEG (front) and MRI (rear projection) studies.

Shock

Request a detailed protocol

A 100 ms, 200 Hz train of stimulation was administered via a Digitimer constant current stimulator (DS7A; Digitimer, Letchworth Garden City, UK). The transistor–transistor logic (TTL) pulse train used to trigger the train of stimulation was generated via a Grass stimulator (SD9, Warwick, RI). In the MEG experiment, this shock was delivered to the subjects’ right wrist via two 8 mm Ag/AgCl surface cup electrodes (EL258-RT, Goleta, CA), filled with electrolyte gel (GEL100, Goleta, CA). In the MRI experiment, the same shock was delivered to the subjects’ right wrist via two 11 mm Ag/AgCl MRI-safe disposable sticker electrodes (EL508, Goleta, CA). The intensity of the shock could range from 0 mA to 100 mA and was calibrated prior to the experiment to a level that the subject rated as ‘uncomfortable but not painful’.

Acoustic startle stimulus

Request a detailed protocol

During the MEG study, subjects were exposed to several presentation of an acoustic white noise, to trigger an acoustic startle reflex used to assess anxiety. In order to avoid artifact due to the magnets and moving metal in traditional headphone drivers, we engineered a custom pneumatic system for generating white noise with air pressure. First we used a 3D printer to create plastic over-the-ear air vortex generators (vortices). Then, we attached these to tubes connected to a solenoid, which was connected to an air tank. When triggered by a TTL pulse, the solenoid allowed air from the tank to pass through tubes to the air vortex generators. In doing so, the air generated a white noise with an intensity proportional to the air pressure released. Therefore, we calibrated the volume of the white noise to 103 dB by adjusting the pressure on the air tank regulator, and testing the intensity with a sound pressure level meter.

Response devices

Request a detailed protocol

During the experiment, subjects had continuous access to an online visual analogue scale that they could use to continuously update their anxiety rating during the task. Subjects controlled this scale using the response device provided for each experiment. In the MEG experiment, subjects used a custom-built fiber optic joystick. In the MRI experiment, subjects used a 4-button fiber optic response device (Current Designs, Philadelphia, PA).

Eyeblink reflex

Request a detailed protocol

The acoustic startle reflex was measured during the MEG study via electromyographic (EMG) activity of the eyeblink reflex recorded via 2 8 mm Ag/AgCl surface cup electrodes (EL258-RT, Goleta, CA), filled with electrolyte gel (GEL100, Goleta, CA) placed below the right eye over the orbicularis oculi muscle (Blumenthal et al., 2005). EMG was recorded at 600 Hz via the MEG system amplifier and analyzed using a custom MATLAB script (see Source code 1). The EMG signal was extracted from the recordings, bandpass filtered (30–500 Hz), rectified, and smoothed with a 20 ms time constant. The peak startle/eyeblink magnitude was determined during the 20–100 ms after the onset of the noise presentation. The peaks were then transformed to z-scores and converted to t-scores within-subjects to reduce large inter-individual differences in the overall magnitude of the startle reflex (Blumenthal et al., 2005), and to facilitate comparison with the online anxiety ratings.

Online ratings

Request a detailed protocol

Subjects reported their anxiety level continuously throughout the experiment using the response device. They also received continuous feedback via a colored circle in the center of the screen, updated at 60 Hz (A colored circle, as opposed to a moving cursor [Schultz et al., 2012; Balderston and Helmstetter, 2010], was chosen to minimize eye movements in the MEG study). As they updated their rating it updated the color of the circle, which ranged from white to red with 256 possible intermediate hues. Values representing these hues were used to numerically represent the subjects’ current anxiety level. During the MEG study, these ratings were sampled just prior to each startle probe for comparison with the startle magnitudes. Because startle probes were not presented during the MRI study, online ratings were sampled once per repetition (TR). The values were then transformed to z-scores and converted to t-scores within-subjects to facilitate comparison with the startle responses.

Procedure

MEG

Request a detailed protocol

On the day of the appointment, the subject arrived at the NIH clinical center, and completed the informed consent form, and pre-experiment questionnaire packet in the waiting room. Next, subjects were escorted to the MEG suite and given the instructions for the task. Afterward, the subject was prepped to enter the magnetically shielded room (MSR). Electrodes to deliver the shock and measure heartbeats/eyeblinks were attached and secured to the subject. In addition, head position indicator (HPI) coils were attached to the subject. One each was attached 1.5 cm anterior to the left and right tragi, and one was attached 1.5 cm superior to the nasion. Finally, the vortices were affixed to the subject’s hairnet, and secured over the ears. Once setup was completed, the subject was escorted into the MSR, given the shock workup procedure (Balderston et al., 2017a, 2017b), and raised into place. Recording began with a 2-minresting run to ensure that there were no visual artifacts (breathing, jaw movement, etc.) in the data. After the resting run, the subject was asked the pre-experiment affective rating scales over the intercom. Once complete, the subject began the experiment.

During the experiment, the subject viewed two concentric circles (See Figure 1). The color of the outer circle indicated the type of block (orange = threat; blue = safe). The color of the inner circle indicated the subject’s current level of anxiety. The subject was instructed to keep the color scale updated continuously throughout the experiment with their current level of anxiety. They were instructed that red meant that they were extremely anxious, while white meant that they were not anxious at all. They were also instructed that they could adjust the color to any shade between red and white, depending on their anxiety level.

The experiment consisted of four runs, each of which were ~6 min long, determined dynamically by the length of the trials. Each run began with four presentations of the startle probes to facilitate habituation, and ensure no head movement. Next, we use the HPI coils to localize the subject’s head relative to the sensors. After head localization, the experimental portion of the run began. The experimental portion of the run consisted of two blocks each of safety and threat in alternating order. Each block contained eight presentations of the startle probe. Probes were separated by a variable interprobe interval (min = 6 s; max = 14 s), which was randomly determined for each probe. Shocks were presented during the threat blocks in randomly selected interprobe intervals, 2–4 s after the preceding probe. The number of shocks was randomly determined as well. Each run could contain between 0 and 2 shocks, and each interprobe interval had a 1 in 12 chance of containing a shock (unless the two shock ceiling had been reached).

MRI

Request a detailed protocol

As with the MEG experiment, the subject arrived at the NIH clinical center, and completed the informed consent form, and pre-experiment questionnaire packet in the waiting room. Next subjects were escorted to the MRI control room and given the instructions for the task. Afterward, the subject was prepped for scanning. Electrodes to deliver the shock and record skin conductance were attached and secured to the subject. In addition, vitamin E capsules were attached to the subject in the same locations as the HPI coils during the MEG study, to facilitate coregistration of the structural MRI and MEG. Next the subject was escorted into the scan room, and given the shock workup procedure (Balderston et al., 2017a, 2017b). Subjects were then given ear plugs, situated on the table, and connected to the pulse oximeter and breathing belt. Scanning began with a localizer, a T1-weighted structural scan, and two 30 s echo planar imaging (EPI) runs. These short EPI runs had opposing phase-encoding directions, and were used for EPI distortion correction (See below). After these, the subject was asked the pre-experiment affective rating scales over the intercom. Once complete, the subject began the experiment.

As in the MEG study, the subject viewed two concentric circles. The color of the outer circle indicated the type of block (orange = threat; blue = safe). The color of the inner circle indicated the subject’s current level of anxiety. The instructions for the task were the same as the MEG study as well, with the exception of the difference in response device.

The experiment consisted of four runs, each of which were ~8 min (490 s) long. Each run consisted of two blocks each of safety and threat in alternating order. The number of shocks was randomly determined as well. Each run could contain between 0 and 3 shocks at random intervals. There were no startle probes presented.

Affective rating scales

Request a detailed protocol

Prior to each experiment, subjects completed several standard psychological questionnaires, including the Spielberger State-Trait Anxiety Index (STAI) (Spielberger, 1987), the Anxiety Sensitivity Index (ASI) (Peterson and Heilbronner, 1987), the Beck Anxiety Inventory (BAI) (Beck et al., 1988), and the Beck Depression Inventory (BDI) (Beck and Steer, 1987). At the start of both experiments and after each magnetoencephalography (MEG)/functional magnetic resonance imaging (fMRI) run, subjects were given a set of affective rating scales: (1) How anxious are you (1 = not anxious, 9 = extremely anxious)? (2) How afraid are you (1 = not afraid, 9 = extremely afraid)? (3) How would you rate the intensity of the electrical stimulation (1 = not painful at all, 9 = uncomfortable but not painful)?

MRI acquisition

Request a detailed protocol

We collected four runs, each containing 245 multi-echo EPI images, using a 3T Siemens MAGNETOM Skyra (Erlangen,Germany) fMRI system, and a 32-channel head coil. For each image, we collected 32 interleaved 3 mm slices (matrix = 64 mm×64 mm; FOV = 192 × 192) parallel to the AC-PC line (TR = 2 s; TEs = 12 ms, 24.48 ms, 36.96 ms; flip angle = 70°). These 32 slices covered the entire cerebrum, but did not cover the most posterior parts of the brainstem and cerebellum. Slices were collected with an anterior-to-posterior phase encoding direction. We also collected two, 10 image multi-echo EPI series with the same parameters and the same field of view; however, one of these series was collected with a posterior-to-anterior phase encoding direction, and these series were used to correct for EPI distortion in the phase encoding direction (Morgan et al., 2004). We also acquired a multi-echo T1-weighted MPRAGE (TR = 2530 ms; TEs = 1.69 ms, 3.55 ms, 5.41 ms, 7.27 ms; flip angle = 7°). We acquired 176, 1 mm axial slices (matrix = 256 mm × 256 mm; field of view (FOV) = 256 mm × 256 mm), which were later co-registered to the EPI images.

MRI preprocessing

Request a detailed protocol

Functional images were preprocessed and analyzed using the AFNI software package (see Source code 1 for processing and analysis scripts) (Cox, 1996). EPI images for each run and each echo were first reconstructed, despiked (i.e. single voxel outliers were truncated), slice-time corrected, and then deobliqued. Then each volume in the series was registered to the first volume, and skull-stripped.

Preprocessed images were then entered into a multi-echo-independent components analysis using the meica.py script distributed with the AFNI software package (Kundu et al., 2012). This analysis uses the T2* decay of BOLD signals, measured across the echoes to denoise the timeseries. The analysis first decomposes the timeseries into independent components using FastICA. Then it determines whether signal intensity across echoes decays in a manner consistent with what is expected from BOLD data. Components that fit the model are kept, components that do not (i.e. components where the signal intensity does not decay across echoes) are discarded. A new denoised timeseries is then synthesized from the components not discarded. This technique has been previously shown to robustly remove sources of noise corresponding to motion, physiology, and scanner artifact (Kundu et al., 2012).

The denoised timeseries for all runs are then registered to the first run, scaled, and further denoised using a general linear model with regressors of no interest. Regressors of no interest included the six motion parameters from the volume registration step, up-through third-order polynomials to model baseline drift, and hemodynamic response functions (HRF)s corresponding to button presses and shock deliveries. In addition, images where the derivative of the motion regressors from volume registration step had a Euclidean norm above 0.5 mm were censored (‘scrubbed’) from further analyses. All remaining images from the safe and threat blocks were concatonated into separate timeseries; however, as part of the denoising procedure, we removed neural responses related to both button presses and transitions from one block type to another.

To correct for geometric distortion of the EPI images, the forward and reverse phase-encoding blips are first averaged across time, skull-stripped, and then rigid-body aligned with the reference image for the other EPI timeseries. Then, the two blips are non-linearly aligned to each other using the ‘plusminus’ flag in the AFNI program 3dQwarp, so that the resulting image is ‘in the middle’ of the two. This reference image can then be registered to the T1 image, and the voxelwise displacement map for the forward blip is then saved, so that it can be applied to the EPI timeseries. Although a standard linear registration approach would have possibly yielded similar results at the group level, the nonlinear approach we used leads to better T1/EPI within-subject registration, and has been shown to perform even better than when using fieldmaps (Hong et al., 2015).

To align the EPI data to Montreal Neurological Institute (MNI) space and to mask out non-grey matter voxels, the T1 data was processed as follows. First, the volumes for the T1 echoes were averaged, and the resulting volume was run through the standard Freesurfer processing pipeline (Desikan et al., 2006; Fischl et al., 2004). Next, the skull-stripped anatomy is non-linearly registered to MNI space using the MNI_avg152T1 template distributed with AFNI. In addition, masks that included all cortical and sub-cortical grey matter atlas regions for each subject. These were warped to MNI space, downsampled to the EPI resolution, and dilated by 1 voxel. A group grey matter mask was created by averaging these binary masks, and thresholding out voxels with less than 2/3 overlap (Torrisi et al., 2015). Next the original space, skull-stripped anatomy was aligned to the distortion corrected EPI reference image in two steps. The first was a simple affine transformation. The second was a non-linear transformation using a criterion based on the local Pearson correlation (Saad et al., 2009), and the inverse displacement map was saved. Finally, the following warps were applied to the EPI timeseries data: inverse warp from the distortion correction step, affine transformation matrix from the T1 alignment step, the inverse warp from the T1 alignment step, and the non-linear warp from the T1 to MNI transformation. Finally, the EPI data were masked with the group grey matter mask, and blurred within this mask using a 6-mm FWHM Gaussian kernel.

MRI analysis

Request a detailed protocol

To analyze the MRI recordings, we used a global brain correlation (GBC) approach (Cole et al., 2010). In brief, we correlated each voxel in our grey matter mask with each other voxel in the mask, applied the Fisher’s Z transform, and summed across correlations. The result was a map where the value in each voxel reflected the strength of the correlation between that voxel and the rest of the brain. We conducted this GBC analysis for each subject and each condition, and used these values for further analysis. First, to identify changes in GBC across the entire brain, we averaged across voxels for safe and threat. We then conducted a paired-sample t-test on these averages. Next, to identify changes in GBC at the regional level, we analyzed the voxel-wise GBC for safe and threat. For this, we conducted a voxel-wise paired-sample t-test on the GBC values. Finally, to identify changes in connectivity with highly connected regions, we identified regions of interest (ROI)s from the voxel-wise GBC threat > safe t-test, and correlated the timecourse of activity in these ROIs with every voxel in the brain for safe and threat, and applied the Fisher’s Z. We then conducted paired-sample t-tests on the resulting connectivity maps.

We used Monte Carlo simulations and a cluster-based method to correct for multiple comparisons across voxels. First, we estimated the smoothness in our residual timeseries using a Gaussian plus mono-exponential shaped function implemented by the ‘-acf’ option in the AFNI program 3dFWHMx, which addresses recent concerns over inflated Type one error in studies using the cluster correction method (Cox et al., 2016). We calculated smoothness for each subject, and averaged this across subjects. Next, we simulated 10,000 random statistical parametric maps in 3dClustSim with a smoothness matching that of the original timeseries. For each simulation, we thresholded at a voxel-wise alpha of 0.005, and extracted the largest cluster. We then compared our test statistics to the distribution of clusters across all simulations to identify a minimum cluster size threshold of 80, corresponding to a two-tailed alpha of 0.05.

MEG acquisition

Request a detailed protocol

We recorded neuromagnetic activity at 600 Hz from 271 radial first-order gradiometers using a 275 channel CTF-OMEGA whole-head magnetometer (VSM MedTech, Ltd., Canada). Recording took place in a magnetically shielded room (Vacuumschmelze, Germany), and Synthetic third-order gradient balancing was used for active noise cancellation (Vrba and Robinson, 2001).

MEG preprocessing

Request a detailed protocol

MEG recordings were preprocessed and analyzed using the Fieldtrip toolbox in MATLAB (See Source code 1 for processing and analysis scripts) (Oostenveld et al., 2011). First, movement within runs was checked by comparing the position of the HPI coils from the beginning of the run to the end of the run, and any run with a root mean square movement value above 1 cm was excluded from the analysis. Next, startle probe onsets were identified, and the 2 s window prior to each trial was extracted, demeaned, and detrended. Next, muscle artifacts were identified using the ft_artifact_muscle function in Fieldtrip. Trials where muscle movements were identified were removed. Next the recordings were low-pass filtered with a 90 Hz cutoff, and notch filtered at 60 Hz to remove line noise. The recordings were then downsampled to 300 Hz, and submitted to an independent components analysis. Components were visually inspected, and those with a topography, and timecourse consistent with either blinks or heartbeats were identified for removal (typically 1–2 per artifact). Rejected components were projected out of the dataset, and the singular value decomposition of the covariance matrix was inspected to determine the regularization factor (lambda). A lambda of 5% was found to be sufficient to regularize the covariance matrix for source analysis following the removal of the rejected components. The result was a dataset cleaned for blink and heartbeat artifacts.

MEG frequency analysis

Request a detailed protocol

Because we were specifically interested in alpha oscillations, we identified each subject’s individual alpha frequency (IAF). We began by transforming the timeseries data into the frequency domain using a multi-taper fast Fourier transform (mtmfft) based on a set of discrete slepian sequences. In this initial transformation, we used a frequency window of 1 Hz to 20 Hz, and a single taper per frequency. We then averaged these spectrograms across sensors and across trials for each subject, and identified the largest local maxima in the average spectrogram for each subjects. For the majority of subjects (24/28), the largest peak occurred in the alpha frequency band (8–12 Hz). For all other subjects, the average IAF was used for further analyses. Next we conducted a second mtmfft, using a 2 Hz window centered around each subject’s IAF, with two tapers per frequency and averaged the resulting power estimates across frequency for each sensor and each trial. This mean IAF power estimate was used in both the sensor space analyses and the source space analysis.

MEG forward model

Request a detailed protocol

Single subject T1 images were used to generate the forward model for the MEG source analysis. First the T1 images were aligned to a single subject template in MNI space. Next the brain surface was extracted, and a single shell head model was generated from this surface. Then a source model was created using a single subject template with current dipoles placed along a regular 8 mm grid inside the brain surface. The single subject MNI space images were aligned to CTF space (i.e. coregistered to the sensors) manually, guided by the vitamin E capsules placed over the fiducial points. The resulting transformation matrix was applied to both the head model and the source model, and alignment between the sensors, head model, and source model was visually inspected. Finally, leadfields were then created using the location of the sensors, head model, and source model dipole locations.

MEG inverse model

Request a detailed protocol

Because we were interested in frequency information (as opposed to time), we used the dynamic imaging of coherent sources (DICS) (Gross et al., 2001) technique to localize the sources of our recordings. We began by computing the cross spectral density matrix (CSD) from the frequencies of interest from all trials in the analysis. We then estimated the beamformer filter using the CSD, leadfields (with fixed orientations), headmodel, and gradiometer locations. Once estimated, this common filter was applied to the safe and threat conditions independently.

MEG sensor-level analyses

Request a detailed protocol

We compared IAF power in safe vs. threat. For this, we averaged the IAF power across trials independently for safe and threat for each subject and each sensor, and conducted a paired sample t-test on these averages. For both analyses, we used Monte Carlo simulations and a cluster-based method to correct for multiple comparisons across sensors. We calculated 1000 random permutations, where condition labels were shuffled across subjects. For each permutation, we thresholded the shuffled results at a sensor-level alpha of 0.005, summed the t-value across sensors in each cluster, and extracted the largest summed t-value. We then compared our test-statistic to the distribution of summed t-values and discarded any clusters where the summed t-value was smaller than the summed t-value corresponding to a two-tailed alpha of 0.05.

MEG source-level analyses

Request a detailed protocol

As with the sensor-level analyses, to compare power in the safe and threat conditions, we projected the average IAF power into source space independently for each condition using the common filter calculated from all trials. We then conducted a paired-sample t-test on these power estimates at each voxel within the source model. As with the sensor-space data we used Monte Carlo simulations and a cluster-based method to correct for multiple comparisons across voxels. As before we calculated 1000 random permutations, where condition labels were shuffled across subjects. For each permutation, we thresholded the shuffled results at a source-level alpha of 0.005, summed the t-value across sensors in each cluster, and extracted the largest summed t-value. We then compared our test-statistic to the distribution of summed t-values and discarded any clusters where the summed t-value was smaller than the summed t-value corresponding to a two-tailed alpha of 0.05.

References

  1. 1
    Diagnostic and Statistical Manual of Mental Disorders: Dsm-5
    1. American Psychiatric Association
    (2013)
    Amer Psychiatric Pub Incorporated.
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
    BDI, Beck Depression Inventory: Manual
    1. AT Beck
    2. RA Steer
    (1987)
    San Antonio: Psychological Corporation.
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
    Structured Clinical Interview for DSM-IV® Axis I Disorders (SCID-I), Clinician Version, Administration Booklet
    1. MB First
    2. RL Spitzer
    3. M Gibbon
    4. JBW Williams
    (2012)
    American Psychiatric Publishing.
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83
  84. 84
  85. 85
  86. 86
  87. 87
  88. 88
  89. 89
  90. 90
  91. 91
  92. 92
  93. 93
  94. 94
  95. 95
  96. 96
  97. 97
  98. 98
  99. 99
  100. 100
  101. 101
  102. 102
  103. 103
  104. 104
  105. 105
  106. 106
  107. 107
    Conditioned Reflex: An Investigation of the Physiological Activity of the Cerebral Cortex
    1. I Pavlov
    (1927)
    Oxford Univ Press.
  108. 108
  109. 109
  110. 110
  111. 111
  112. 112
  113. 113
  114. 114
  115. 115
  116. 116
  117. 117
  118. 118
  119. 119
  120. 120
  121. 121
  122. 122
  123. 123
  124. 124
  125. 125
    The angular gyrus: multiple functions and multiple subdivisions
    1. ML Seghier
    (2013)
    The Neuroscientist : A Review Journal Bringing Neurobiology, Neurology and Psychiatry 19:43–61.
    https://doi.org/10.1177/1073858412440596
  126. 126
  127. 127
  128. 128
    State-trait anxiety inventory
    1. CD Spielberger
    (1987)
    Anxiety 19:2009.
  129. 129
  130. 130
  131. 131
  132. 132
  133. 133
  134. 134
  135. 135
  136. 136
  137. 137
  138. 138
  139. 139
  140. 140
  141. 141

Decision letter

  1. Matthew J Brookes
    Reviewing Editor; University of Nottingham, United Kingdom

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Threat of shock increases excitability and connectivity of the intraparietal sulcus" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by David Van Essen as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Krish Singh (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Your paper was assessed positively by all three reviewers, with all three saying that this is a valuable addition to the literature and of broad scientific interest. In particular, the multi-modal aspect was well received. I am therefore happy to provisionally recommend publication, however a number of significant revisions must be successfully carried out first. On receiving your revised manuscript with explanatory comments, the Reviewing Editor will decide whether the manuscript needs to be seen again by the reviewers.

Essential revisions:

For fMRI, all of the technical concerns raised by reviewer 2 must be addressed comprehensively. These are copied in full below.

For MEG, both reviewers 1 and 3 suggest that a functional connectivity analysis would be a better candidate for comparison with your fMRI results than the current approach. I therefore strongly recommend that you undertake such analyses; here, a key confound is signal leakage between source space regions and an appropriate correction method must be employed. Though multiple schemes are available I would suggest using the approach published by Colclough, Brookes, Smith and Woolrich, 2015. (referenced by reviewer 1) but perhaps applied to a denser parcellation (e.g. the AAL parcellation). Again for completeness the reviewers comments are copied in full below.

Finally, a quantitative analysis of the similarity of the regions identified by MEG and fMRI must be given.

Reviewer #1:

The paper by Balderston et al. describes a study in which both fMRI and MEG are used to investigate the neural correlates of anxiety. To this end, a threat of shock paradigm is employed to increase anxiety in subjects; this increase is then assessed behaviourally in multiple ways. Both fMRI and MEG data are acquired during the paradigm and analysed in an unbiased way; briefly, in fMRI an 'all-to-all' connectivity approach is used to measure connectivity among voxels. In MEG, data are analysed in multiple frequency bands and changes in oscillatory power between task blocks (safe versus threat of shock) characterised. The interesting result is that both the fMRI connectivity analysis and the α band power assessment both implicate the intraparietal sulcus as a key hub in processing anxiety.

Overall I believe this to be a very interesting and exciting paper. The primary finding is not only impactful for basic science but may also, potentially, have significant clinical relevance. The paper is well written and the analyses, for the most part, solid. However I do have a number of suggestions that I think would improve the paper.

Major comments

1) It is assumed that the fMRI and MEG responses occur in the same location, however this was never analysed quantitatively. The distance between the fMRI cluster and the MEG peak location should be quantified. What are the chances that this difference could occur by chance?

2) It is not made clear in the Introduction why complementary analysis is not used – e.g. the MEG is used to assess oscillatory changes whereas fMRI is used to assess connectivity. Why? I realise of course the links that have been made between oscillations and connectivity but this is not well discussed in the paper (I personally didn't find the relation to the simultaneous EEG-fMRI literature very convincing). Why not just do a standard GLM approach to the BOLD analysis to get regions of increased BOLD response during threat and compare that with changes in α oscillations? In short, the analysis pipeline used should be better motivated.

3) Related to the above question, I found the MEG analysis quite limited. Wouldn't it have been better to take advantage of the huge steps forward in MEG connectivity methods that have been made recently and employ one such approach to measure whole brain α band connectivity – and then compare this to the fMRI connectivity. Its true that you couldn't do this at the voxel level spatial scale but with a brain parcellation it should be easily possible. See papers by e.g. Colclough, Brookes, Smith and Woolrich, 2015 or O'Neill et al., http://iopscience.iop.org/article/10.1088/0031-9155/60/21/R271

4) I felt there should be more justification for looking at the α band, which seems to have been chosen specifically, with the other frequency bands treated as an afterthought. Given the apparent close links between β oscillations and functional connectivity why was the β band not chosen? Again just saying that α was the strongest signal at many sensors didn't really convince me.

5) Can the authors give some explanation as to why the α response was lateralised?

Reviewer #2:

General

This study is a nice example of an investigation into neuronal mechanisms of anxiety that uses two different complementary modalities. However, some of the fMRI analysis approaches need further investigation to ensure that the results were not driven by noise, and the manuscript is longer than necessary in some places.

Major comments

1) The within-network results do not pass multiple comparison correction and do not add much to the story. I would therefore suggest that the authors remove these results. In general, the manuscript should be shortened. Specifically the Discussion would benefit from being more concisely written, and the Materials and method could also be reduced (specifically the subject numbers and the 'on the day of the appointment…' text).

2) Was there a difference in motion parameters between the safe and threat blocks in fMRI? If so, this could lead to a shift between short-distance and long-distance connectivity, which could drive the results. The scrubbing that was performed most likely avoids this possibility, but it would be of interest to present the difference in motion parameters, and to potentially run further analyses to explicitly exclude this possibility (for example, by selecting periods or subjects in which the motion was matched between the two conditions).

3) How many time points were removed during scrubbing, and was there a significant difference between the safe and shock blocks for this? If so, the difference in power could bias the results, and it would be worth excluding this possibility by matching the number of timepoints included in shock and safe blocs (within participants).

4) The change in global connectivity during the shock period could occur as a result of the actual shock stimuli that were not present in the safe periods. I know the HRF responses to shocks and button presses were modeled and removed, but this may not capture the changes in functional connectivity resulting from shock stimuli. Would it be possible to repeat the analysis using only shock blocks that contained zero shocks to ensure that this is not driving the results?

5) Cluster-based corrections with an initial cluster forming threshold above 0.001 have been shown to suffer from inflated false positives (Eklund et al., 2016). The authors should at least point this out in the limitations, and should ideally repeat analyses using the latest guidelines.

6) The methods used in this study are described as unbiased and multimodal, which is a strong statement to make given that some of the whole-brain statistics might suffer from a false positive bias, and the different modalities are analysed separately, rather than in a joint multimodal approach.

Reviewer #3:

I really enjoyed reading this paper, which is a valuable contribution to the field. It uses a multi-modal approach of combining fMRI and MEG to reveal changes in both connectivity across the brain (using fMRI) and oscillatory power changes (using MEG) while participants were either in a state of 'Safe' or 'Threat'

To me the results are interesting both in terms of probing the neurophysiological underpinnings of threat/anxiety but also the relationship between fMRI and MEG measures of brain function.

The methods are thorough and well-executed and reveal a good understanding of the state-of-the-art in fMRI and MEG analyses. However I think more could be done with the MEG data (see my comment 1 below).

The paper is well-written and readable and the authors nicely discuss the weaknesses of their approach.

Here are my main comments:

1) It seems as if there is a missed opportunity here in terms of directly comparing connectivity measures in fMRI to those extracted with MEG (albeit at the group-level). Using a sub-sampled atlas approach (say the AAL atlas) the authors could easily compare 'GBC"-type connectivity matrices extracted with fMRI with those extracted with MEG, either using amplitude-amplitude coupling within frequency bands or phase-phase coupling. There are several recent papers using this approach, so it seems odd that the authors did not do this.

2) It is not absolutely clear in the Abstract/Introduction that when the authors talk about neural activity and connectivity in the IPS, they are really talking about haemodynamics measures (i.e. FMRI). This should be made clearer.

3) When presenting the fMRI-GBC (for example in Figure 3) I think it would be useful to show the spatial distribution across the brain of the connectivity measure i.e. after summing across rows but before global summing across the brain. This would be a companion figure to Figure 4, which shows the voxelwise difference in connectivity strength between Threat v. Safe. Based on my comment 1 above it would be great to see similar visualisations based on frequency-specific MEG connectivity maps.

4) In the MEG analysis, it looks like some trials are rejected if contaminated with artefacts. I just wanted to check that the statistical tests performed later (between Safe and Threat) properly account for differing number of trials.

5) In the MEG analysis, did the authors check that there was no difference in head movement between the Safe and Threat blocks?

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Threat of shock increases excitability and connectivity of the intraparietal sulcus" for further consideration at eLife. Your revised article has been favorably evaluated by David Van Essen (Senior editor), a Reviewing editor, and two reviewers.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

We thank the authors for their careful consideration of the comments of all three reviewers, and for the extra analyses that have been undertaken. The paper is much improved as a result of these changes and is now close to publication. However, a few points remain to be addressed.

1) I consider the conjunction map, provided in the response to reviewers, to be critical to the paper and I would very much like to see this in the main paper, perhaps added to Figure 7.

2) Its rather a shame that the MEG connectivity analysis didn't work, but nevertheless I thank the authors for attempting it. I wonder however if some extra clarification could be given: What frequency band was connectivity computed in (sincere apologies if I missed it)? I assume α? I'm also unsure as to what is meant by "we then computed the Hilbert transform using a 500 ms sliding window"? The Hilbert transform shouldn't need a sliding window – when undertaking this analysis we usually compute the HT over all time? Please clarify why this was done. Could the authors please show the MEG adjacency matrices for the safe and thread blocks independently, as well as the difference? This would allow the reader to confirm that the two separate matrices look sensible (one expects large occipital connectivity in the α band, see e.g. Hunt et al., 2016). Finally, given the strong finding of MEG connectivity in the β band, I would like to see the connectivity analysis attempted in this band. I of course read the authors argument that "one might expect a high degree of stationarity in coherence of the β oscillations across time" but I disagree, specifically because multiple papers (e.g. O'Neill et al., 2015; O'Neill et al., 2016; Baker et al., 2014) show that, in fact, β band functional connectivity is highly dynamic.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Threat of shock increases excitability and connectivity of the intraparietal sulcus" for further consideration at eLife. Your revised article has been favorably evaluated by David Van Essen (Senior editor) and the Reviewing editor.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

I thank the authors for going to such lengths on the MEG connectivity analysis. However, something has clearly gone wrong here in the analysis.

By my assessment of the adjacency matrices presented, the case for the α and β bands look virtually identical. However this should not be the case (again see Hunt et al, 2016). It is quite hard for me to judge the spatial signature of connectivity as, unfortunately, in the pdf version the figure quality is too low resolution to read the region names on the axes. However, these matrices do not look correct to me and I would urge the authors to find out why.

In their response, the authors have suggested that they use blocks of only 2 seconds of envelope data. However they also apply downsampling prior to connectivity estimation. This to me seems silly – why apply the downsampling if you have such short blocks. It is well established that connectivity works without the downsampling so perhaps removing that step might improve the adjacency matrices. I also don't understand why such short data setments were used when in fact the blocks are quite long. Why not just try using the whole block?

I remain enthusiastic about this article, but I maintain that the MEG connectivity analysis has not yet been properly carried out. I would want to see this estimated reliably (or a concrete argument on why it cannot be undertaken) put forward before recommending publication.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Threat of shock increases excitability and connectivity of the intraparietal sulcus" for further consideration at eLife. Your revised article has been favorably evaluated by David Van Essen (Senior editor) and the Reviewing editor.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

I would like to thank the authors for their attempt at addressing my concerns. I'm pleased that the authors were able to spot errors in their pipeline (use of incorrect spatial filters) although it does appear that the major difference came from the down-sampling, as suggested in my previous review. Downsampling over such a small time window is obviously incorrect so these analyses should now be discarded.

Unfortunately, the adjacency matrices still don't look right – without seeing them plotted on a brain it's hard to judge exactly where the primary pathways of connectivity are (and I still can’t read the labels clearly in the pdf). As I have said previously, the highest connectivity in the α band should be in the occipital lobe whilst the highest connectivity in the β band should encompass bilateral parietal and occipital connections alongside tempero-parietal and fronto-parietal networks (again as in Hunt et al, 2016). In fact this does not seem to be the case. In what the authors have provided the α network looks like pure noise. There is clearly some structure to the β band connectivity matrix although this doesn't really look like one would expect. So I strongly suspect something is still wrong with the analysis.

Without looking at the data directly it’s hard for me to judge but I suspect that the problem is the short time windows (unless there are other more basic errors (similar to the incorrect spatial filter) which the authors have not discovered). It is known that a reasonable amount of data is required to make MEG connectivity analyses work reliably (Luizzi et al., Optimising experimental design for MEG resting state functional connectivity measurement, NeuroImage 2016) and perhaps this is the reason the adjacency matrices look so poor. In the light of these failed attempts, I suggest that a paragraph is added to the Discussion stating that MEG connectivity analysis would have been a useful means to probe these data; however the study design was poorly set up to make such analyses work.

Please note that this must now be your final attempt to satisfy the Board.

https://doi.org/10.7554/eLife.23608.031

Author response

Essential revisions:

For fMRI, all of the technical concerns raised by reviewer 2 must be addressed comprehensively. These are copied in full below.

Reviewer 2 raised 6 technical concerns in their major issues section. Technical concern 1 was addressed as part of a larger effort to shorten the manuscript, and therefore will be described below in a separate section about manuscript length. Technical concerns 2-4 were all related to ruling out alternative explanations for our global connectivity results. These will be addressed together in the paragraphs below. Finally, technical concerns 5 and 6 were separate issues, and are addressed inline with the reviewer’s original comments.

In technical concerns 2-4, reviewer 2 questioned whether the following could explain our global connectivity results: 2) the difference in motion across block types (i.e. safe vs. threat), 3) the difference in the number of TRs censored across block types, and 4) the presentations of the shocks in the threat blocks. Before we describe the additional steps conducted to rule out these possibilities, we would like to note that we took very thorough steps in the original analysis to account for possible sources of noise in the fMRI data. First, we used the novel multi-echo ICA cleaning approach which identifies and eliminates sources of noise that do not exhibit the expected TE-dependence of the BOLD signal. Second, we used a strict threshold for head motion scrubbing (i.e. censoring TRs with a framewise displacement > 0.5 mm). Third, we used a GLM to remove noise related to the 6 motion parameters from the EPI timeseries. Fourth, we modeled the HRF corresponding to each shock delivery, and included this as a regressor of no interest in above mentioned GLM.

To further address these questions, we began by examining whether motion and TRs censored differed across block type in the original analysis. To quantify motion, we calculated the Euclidean norm of framewise displacement across the six motion parameters after scrubbing, and averaged this across the TRs included in the analysis for each block type. We found that this indeed differed across the safe and threat blocks (t(24) = 15.32; p < 0.001). Similarly, when we count the total TRs included in the analysis (Total TRs – Censored TRs) we found that this also differed across threat and safe blocks (t(24) = -6.89; p < 0.001). Therefore to account for this difference, and to account for the potential effects of the shock delivery, we recalculated global brain connectivity (GBC) using a more selective set of TRs.

In this analysis, we censored the 10 TRs during the threat periods following each shock delivery, and then censored safe TRs so that the total number of TRs was equivalent across the safe and threat blocks, and the Euclidean norm of framewise displacement was as equivalent as possible across the safe and threat blocks (See Author response image 1 for total number of TRs and framewise displacement before and after this correction). Although we were able to fully match the total number of TRs across blocks, the framewise displacement was not fully matched (t(24) = -3.75; p = 0.001). However, it should be noted that the average framewise displacement after this correction was very low (< 1mm). To address this remaining difference in framewise displacement, we calculated the difference in framewise displacement across the safe and threat blocks, and used this value as a covariate in an ANCOVA looking at the effect of threat on GBC across this more selective set of TRs.

Since the purpose of this analysis was to determine whether our original results remained after this more stringent scrubbing, we sampled GBC within the same ROIs identified by the original whole brain GBC analysis. We found that even after censoring out 10 TRs following the shock, matching the total number of TRs censored, closely matching framewise displacement, and factoring out the remaining differences in framewise displacement, we still see a significant difference in GBC as a function of threat in all 3 ROIs identified in the original analysis. These results suggest that our initial findings were not due to differences in motion, censoring, or residual neural activity evoked by the shock. We understand that these are important controls, and have therefore revised the original version of the manuscript to include the details of this analysis.

For MEG, both reviewers 1 and 3 suggest that a functional connectivity analysis would be a better candidate for comparison with your fMRI results than the current approach. I therefore strongly recommend that you undertake such analyses; here, a key confound is signal leakage between source space regions and an appropriate correction method must be employed. Though multiple schemes are available I would suggest using the approach published by Colclough Brookes, Smith and Woolrich, 2015. (referenced by reviewer 1) but perhaps applied to a denser parcellation (e.g. the AAL parcellation). Again for completeness the reviewers comments are copied in full below.

We agree with the reviewers that this work would benefit from a connectivity analysis of the MEG data, and recognize that signal leakage is a serious issue affecting such analysis. Therefore, we have chosen to follow the reviewers’ advice and conduct an analysis of the differences in MEG connectivity as a function of threat, while accounting for signal leakage using the symmetric orthogonalisation method detailed by Colclough, Brookes, Smith and Woolrich, 2015.

We began by computing the virtual MEG channels using the single trial timecourses and leadfields computed in the original manuscript. Then, using the toolbox provided by Dr. Colclough (https://github.com/OHBA-analysis/MEG-ROI-nets), we downsampled the source space MEG data to the AAL atlas, and computed the symmetric ortrhogonalization of the atlas space signals to remove any 0-phase correlations. We then computed the Hilbert transform using a 500 ms sliding window, low-pass filtered the data to 0.5 Hz, and downsampled the data to 1 Hz. Next we concatonated these envelope timeseries across all trials for safe and threat separately for each subject, and computed the cross correlation matrix for each subject and each block separately. To determine whether there were any differences in connectivity between safe and threat, we computed paired-sample t-tests at each cell of the cross correlation matrix. We then conducted 10,000 permutation tests and extracted the maximum t-value for each permutation in order to build a null distribution for these t-tests. Finally, we identified a corrected t-value at the 95 percentile of this distribution (t = 5.5), and thresholded the original t-map using this corrected t-value. As can be seen from the graphs in Author response image 2, there were no significant changes in pairwise functional connectivity as a function of threat based on this analysis and correction method.

Given that the overall suggestion was to reduce the length of the manuscript, we have decided not to include this analysis in the revised submission. We have however revised the limitations section of the manuscript to include a discussion of this point. We are happy to make additional changes if requested.

Finally, a quantitative analysis of the similarity of the regions identified by MEG and fMRI must be given.

We agree with this reviewer that it is important to demonstrate the correspondence between the regions identified by the MEG and fMRI analyses. Accordingly, we have prepared a conjunction map of these two results, Author response image 3. There are two major points illustrated by this map. First, although the fMRI results are bilateral, the MEG results are unilateral. This point will be addressed below. Second, the fMRI cluster in the left IPS substantially overlaps with the source space MEG cluster (46 of 81 voxels). Due to the recommendation to reduce the length of the manuscript, we did not include this figure in the revised version of the manuscript. However, we did revise the Results section to reflect this finding. If requested, we would be happy to add a figure with this result to the manuscript.

Reviewer #1:

[…]

Major comments

1) It is assumed that the fMRI and MEG responses occur in the same location, however this was never analysed quantitatively. The distance between the fMRI cluster and the MEG peak location should be quantified. What are the chances that this difference could occur by chance?

This point has been addressed in the general comments to the editor, see above.

2) It is not made clear in the Introduction why complementary analysis is not used – e.g. the MEG is used to assess oscillatory changes whereas fMRI is used to assess connectivity. Why? I realise of course the links that have been made between oscillations and connectivity but this is not well discussed in the paper (I personally didn't find the relation to the simultaneous EEG-fMRI literature very convincing). Why not just do a standard GLM approach to the BOLD analysis to get regions of increased BOLD response during threat and compare that with changes in α oscillations? In short, the analysis pipeline used should be better motivated.

We agree with this reviewer that a standard GLM approach to analyze the BOLD data might be a better match for the analysis of the MEG oscillation data. However, our blocks were quite long for fMRI data ~2 min, and therefore subject to baseline drifts. Although there are approaches to minimize the contributions of these baseline drifts to activation estimates (including the multi-echo ICA denoising technique used here), it is unclear whether these techniques are robust enough to warrant conducting the additional GLM approach suggested. Given the recommendation to reduce the overall length of the manuscript, we have decided not to include this analysis. However, we are happy to conduct the GLM analysis, and include it in the paper if deemed necessary.

3) Related to the above question, I found the MEG analysis quite limited. Wouldn't it have been better to take advantage of the huge steps forward in MEG connectivity methods that have been made recently and employ one such approach to measure whole brain α band connectivity – and then compare this to the fMRI connectivity. Its true that you couldn't do this at the voxel level spatial scale but with a brain parcellation it should be easily possible. See papers by e.g. Colclough, Brookes, Smith and Woolrich, 2015 or O'Neill et al., http://iopscience.iop.org/article/10.1088/0031-9155/60/21/R271

We agree that MEG connectivity is an interesting topic, and have undertaken the analysis described in Colclough, Brookes, Smith and Woolrich, 2015. This analysis is described in the general comments to the editor, see above.

4) I felt there should be more justification for looking at the α band, which seems to have been chosen specifically, with the other frequency bands treated as an afterthought. Given the apparent close links between β oscillations and functional connectivity why was the β band not chosen? Again just saying that α was the strongest signal at many sensors didn't really convince me.

We agree that there is a strong connection between functional connectivity and β oscillations. Indeed, a careful examination of spontaneous β oscillations can uncover a set of resting state networks that is remarkably similar to those identified from spontaneous fluctuations in BOLD data. However, given that resting state networks identified from BOLD data are highly correlated with independent measures of anatomical connectivity, one might expect a high degree of stationarity in coherence of the β oscillations across time. This might be acceptable if the purpose of this study were to identify correlations between resting state networks and anxiety across subjects. However, the purpose of this study was to determine how a within-subject experimental manipulation of anxiety affects spontaneous measures of neural activity/connectivity.

Although the focus of the paper was on α, our initial approach (described in the Introduction) was to examine all frequency bands independently. The focus on α emerged out of the observations that, 1) α was the strongest signal in the recordings, 2) α showed the largest power changes as a function of threat, 3) α was the only frequency band that showed consistent results at both the sensor and the source level, 4) the source space results aligned nicely with the corresponding fMRI GBC data.

It should also be noted that there are several theoretical reasons to focus on α as well. For instance, it is well established that spontaneous neural activity at rest is dominated by oscillations in the α (8-12 Hz) range, which are most prominent when the subject is in an alert state of restful relaxation, and that α asymmetries can reflect differences in arousal across subjects. Theoretical models of α function suggest that α oscillations are generated by coherent activity in local inhibitory interneurons, and that decreases in α power reflects increases in cortical excitability. Consistent with these theories, studies collecting simultaneous measures of electroencephalography (EEG) and fMRI have shown that α power is negatively correlated with functional connectivity.

Given these points, we believe that it is appropriate to focus primarily on the α results. However, we agree that this motivation should be clarified. This is especially true given that we have chosen to remove the results from the other frequency bands, in order to respond to the request to shorten the manuscript. Therefore we have revised the Discussion section of the manuscript to clarify our motivation for focusing primarily on the α results.

5) Can the authors give some explanation as to why the α response was lateralised?

There is actually a robust literature on α asymmetry and emotional processing. According to prominent models, left dominant frontal α is associated with positive affect and/or approach behavior, while right dominant α is associated with negative affect and/or avoidance behavior (Davidson 2004; Harmon-Jones, Gable and Peterson, 2010). Like the valence model of α asymmetry, one prominent model on parietal α asymmetry is rooted in the arousal-valence model of emotional processing (Heller, 1993). According to this model, right dominant parietal α is associated with increased arousal. Consistent with this theory, we find that threat of shock, which increases arousal, also reduces left parietal α, resulting in a right dominant profile. However, more research should be conducted to specifically test this hypothesis. We understand that this is an important point, and have revised the manuscript to include a brief discussion of this aspect of the findings.

Reviewer #2:

General

This study is a nice example of an investigation into neuronal mechanisms of anxiety that uses two different complementary modalities. However, some of the fMRI analysis approaches need further investigation to ensure that the results were not driven by noise, and the manuscript is longer than necessary in some places.

Major comments

1) The within-network results do not pass multiple comparison correction and do not add much to the story. I would therefore suggest that the authors remove these results. In general, the manuscript should be shortened. Specifically the Discussion would benefit from being more concisely written, and the Materials and method could also be reduced (specifically the subject numbers and the 'on the day of the appointment…' text).

This point has been addressed in the general comments to the editor, see above.

2) Was there a difference in motion parameters between the safe and threat blocks in fMRI? If so, this could lead to a shift between short-distance and long-distance connectivity, which could drive the results. The scrubbing that was performed most likely avoids this possibility, but it would be of interest to present the difference in motion parameters, and to potentially run further analyses to explicitly exclude this possibility (for example, by selecting periods or subjects in which the motion was matched between the two conditions).

As with point 1, this point has been addressed in the general comments to the editor, see above.

3) How many time points were removed during scrubbing, and was there a significant difference between the safe and shock blocks for this? If so, the difference in power could bias the results, and it would be worth excluding this possibility by matching the number of timepoints included in shock and safe blocs (within participants).

As with point 1 and 2, this point has been addressed in the general comments to the editor, see above.

4) The change in global connectivity during the shock period could occur as a result of the actual shock stimuli that were not present in the safe periods. I know the HRF responses to shocks and button presses were modeled and removed, but this may not capture the changes in functional connectivity resulting from shock stimuli. Would it be possible to repeat the analysis using only shock blocks that contained zero shocks to ensure that this is not driving the results?

As with points 1, 2, and 3, this point has been addressed in the general comments to the editor, see above.

5) Cluster-based corrections with an initial cluster forming threshold above 0.001 have been shown to suffer from inflated false positives (Eklund et al., 2016). The authors should at least point this out in the limitations, and should ideally repeat analyses using the latest guidelines.

Indeed, we agree that several common approaches to multiple comparisons suffer from inflated Type 1 error, as indicated by Eklund et al., 2016. Cluster thresholding, as implemented by the AFNI program 3dClustSim, had, until recently, used Gaussian estimates of the noise smoothness in the monte carlo simulations, which do not appropriately model the noise smoothness of BOLD data. The AFNI team has since provided new option that models the noise smoothness using a Gaussian plus mono-exponential spatial AutoCorrelation Function (ACF), which addresses the criticisms raised by the Eklund paper.

See:

http://biorxiv.org/content/early/2016/07/26/065862

In the current work, we used the newly created autocorrelation function to estimate the smoothness of the noise in our data, which was used to identify minimum cluster-size threshold reported in the manuscript. We realize that this was not sufficiently explained in the original version of the manuscript, and have therefore clarified this methodological detail in the revised version.

In addition, although our voxelwise p-value is higher than recommended by the reviewer, this is offset by a larger minimum cluster-size threshold. This, combined with the fact that our results are bilateral and consistent with the MEG data, leads us to conclude that it is unlikely that our results are due to Type 1 error.

6) The methods used in this study are described as unbiased and multimodal, which is a strong statement to make given that some of the whole-brain statistics might suffer from a false positive bias, and the different modalities are analysed separately, rather than in a joint multimodal approach.

The main finding of the current work is that aspects of IPS functioning are affected by the threat of shock manipulation. We show this in two ways: first we show that α oscillations are decreased by threat of shock (MEG), second we show that global IPS functional connectivity is increased by threat of shock (fMRI). These two independent findings are significant at the whole brain level, corrected for multiple comparisons. For that reason, we feel that it is appropriate to refer to these findings as both unbiased and multimodal.

Reviewer #3:

[…]

Here are my main comments:

1) It seems as if there is a missed opportunity here in terms of directly comparing connectivity measures in fMRI to those extracted with MEG (albeit at the group-level). Using a sub-sampled atlas approach (say the AAL atlas) the authors could easily compare 'GBC"-type connectivity matrices extracted with fMRI with those extracted with MEG, either using amplitude-amplitude coupling within frequency bands or phase-phase coupling. There are several recent papers using this approach, so it seems odd that the authors did not do this.

This point has been addressed in the general comments to the editor, see above.

2) It is not absolutely clear in the Abstract/Introduction that when the authors talk about neural activity and connectivity in the IPS, they are really talking about haemodynamics measures (i.e. FMRI). This should be made clearer.

In the Introduction and Abstract we refer to activity and connectivity. When we use the term activity, we are referring to reductions in α power (MEG) that reflect increases in excitability. When we use the term connectivity, we are referring to the functional connectivity measured with fMRI. These points have been clarified in the revised version of the manuscript (See Abstract and Introduction).

3) When presenting the fMRI-GBC (for example in Figure 3) I think it would be useful to show the spatial distribution across the brain of the connectivity measure i.e. after summing across rows but before global summing across the brain. This would be a companion figure to Figure 4, which shows the voxelwise difference in connectivity strength between Threat v. Safe. Based on my comment 1 above it would be great to see similar visualisations based on frequency-specific MEG connectivity maps.

We agree that the revised version of the manuscript would benefit from this figure. Therefore, we have updated Figure 3A to reflect this change. As for the MEG connectivity maps, we conducted the functional connectivity analysis recommended by the editor, which yielded no significant results. Given the overall recommendation to shorten the manuscript, we have decide not to include these results in the final work. However, we are happy to add this work if deemed necessary.

4) In the MEG analysis, it looks like some trials are rejected if contaminated with artefacts. I just wanted to check that the statistical tests performed later (between Safe and Threat) properly account for differing number of trials.

Indeed, when we examine the number of trials included in the analysis across conditions, we do see a difference in the number of trials for the safe and threat conditions. In the safe condition there were on average 58.21 ± 5.5 trials, while in the threat condition there were on average 55.32 ± 7.6 trials, which was significantly different across subjects (t(27) = 2.78, p = 0.01). Therefore, we decided to determine whether this difference in trial number impacted our results at the censor and source level. Accordingly, we included the difference in trial number across safe and threat blocks as a covariate in an ANCOVA examining the effect of threat on α responses. At both the sensor (f(1,26) = 7.48, p = 0.01) and at the source (f(1,26) = 17.797, p < 0.001), we find that even with the difference in number of trials covaried out, the effects of threat on α power is still significant, suggesting that this difference did not significantly impact our results. This analysis has been included in the revised version of the manuscript.

5) In the MEG analysis, did the authors check that there was no difference in head movement between the Safe and Threat blocks?

Unfortunately we are unable to determine this from our current data. In this study we determined the amount of motion by comparing the subject’s position before a run to their position after the run. Accordingly, we did not have an online measure of motion within a given run. Each run contained four blocks: two safe and two threat, so we do not have condition specific measures of motion for the MEG study. We agree that this is a limitation of the current work, and have added a statement to the Discussion to address this limitation.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

We thank the authors for their careful consideration of the comments of all three reviewers, and for the extra analyses that have been undertaken. The paper is much improved as a result of these changes and is now close to publication. However, a few points remain to be addressed.

1) I consider the conjunction map, provided in the response to reviewers, to be critical to the paper and I would very much like to see this in the main paper, perhaps added to Figure 7.

We agree that this conjunction map is an important part of the work, and are happy to add it to the main paper. We have added it as Figure 8 in the revised version of the manuscript.

2) Its rather a shame that the MEG connectivity analysis didn't work, but nevertheless I thank the authors for attempting it. I wonder however if some extra clarification could be given: What frequency band was connectivity computed in (sincere apologies if I missed it)? I assume α? I'm also unsure as to what is meant by "we then computed the Hilbert transform using a 500 ms sliding window"? The Hilbert transform shouldn't need a sliding window – when undertaking this analysis we usually compute the HT over all time? Please clarify why this was done. Could the authors please show the MEG adjacency matrices for the safe and thread blocks independently, as well as the difference? This would allow the reader to confirm that the two separate matrices look sensible (one expects large occipital connectivity in the α band, see e.g. Hunt et al., 2016). Finally, given the strong finding of MEG connectivity in the β band, I would like to see the connectivity analysis attempted in this band. I of course read the authors argument that "one might expect a high degree of stationarity in coherence of the β oscillations across time" but I disagree, specifically because multiple papers (e.g. O'Neill et al., 2015; O'Neill et al. 2016; Baker et al. 2014) show that, in fact, β band functional connectivity is highly dynamic.

There were several points raised in the above comment. For clarity, we have numbered these, and addressed them below.

a) This is correct. Given that the other analyses in the paper were focused on the α band, we conducted the connectivity analyses in the α band as well.

b) This statement was made in error. We were describing the functionality of the power envelope calculation implemented in the MEG-ROI-nets toolbox, and mistakenly mixed up the steps. Upon reviewing the methods of the Colclough, Brookes, Smith and Woolrich, 2015, and the matlab functions provided by Dr. Colclough on github (https://github.com/OHBA-analysis/MEG-ROI-nets), it is clear that these functions compute the power envelope by first taking the Hilbert transform, then low-pass filtering the data, and finally downsampling the data using a sliding window approach.

c) We have attached the adjacency matrices for the safe and threat conditions in Author response image 4. Although the matrices do not exactly match those from Hunt et. al 2016, there appears to be α connectivity in the occipital lobe for both safe and threat conditions. One possible reason for the difference across studies is that we calculated connectivity across a series of discrete events (the 2 seconds prior to each startle probe). We understand that these short blocks are not optimal for studying connectivity; however, the study was initially designed to investigate pre-stimulus activity and startle magnitude under threat.

Author response image 4
Α connectivity across AAL regions during 2 second baseline prior to startle probe.

A) Mean α connectivity across AAL regions during safe periods. B) Mean α connectivity across AAL regions during threat periods. C) Unthresholded T-test results comparing α connectivity during safe and threat conditions. Thresholded T-test results comparing α connectivity during safe and threat conditions. Labels on Y-axis correspond to regions of the AAL axis. Labels on the X-axis correspond to groups from AAL atlas (frontal, limbic, occipital, parietal, subcortical, temporal, cerebellum). Boxes in A and B represent α connectivity in the occipital cortices.

https://doi.org/10.7554/eLife.23608.025

d) As requested, we have conducted an initial analysis of connectivity in the β band. We used the exact approach as in the previous study. Briefly, we bandpass filtered the data, downsampled it to the AAL atlas, applied the orthogonalization procedure detailed in Colclough, Brookes, Smith and Woolrich, 2015, computed the Hilbert transform, downsampled the power envelopes, and computed pairwise t-tests at each node of the adjacency matrices. Note that, like α, there were no significant differences in β connectivity from safe to threat after correcting for multiple comparisons. Although it may be interesting to explore changes in MEG connectivity as a function of threat, perhaps future studies can be conducted to examine this question. See Author response image 5.

Because the connectivity analysis involves downsampling the power envelopes, one recommendation for future studies would be to ensure longer windows of clean data by removing the startle probes from the experimental design. We included the startle probes as a way to probe anxiety throughout the safe and threat blocks, however, because blink magnitude differs across conditions, the blink artifact also differs across conditions. Although beamforming can theoretically remove the artifact induced by the blink response, the only way to ensure that the blink artifact does not differentially affect the MEG signal is to remove the contaminated time periods from the analysis, or remove the startle probes from the design to begin with.

Author response image 5
Β connectivity across AAL regions during 2 second baseline prior to startle probe.

A) Mean β connectivity across AAL regions during safe periods. B) Mean β connectivity across AAL regions during threat periods. C) Unthresholded T-test results comparing β connectivity during safe and threat conditions. Thresholded T-test results comparing β connectivity during safe and threat conditions. Labels on Y-axis correspond to regions of the AAL axis. Labels on the X-axis correspond to groups from AAL atlas (frontal, limbic, occipital, parietal, subcortical, temporal, cerebellum).

https://doi.org/10.7554/eLife.23608.026

[Editors' note: further revisions were requested prior to acceptance, as described below.]

I thank the authors for going to such lengths on the MEG connectivity analysis. However, something has clearly gone wrong here in the analysis.

By my assessment of the adjacency matrices presented, the case for the α and β bands look virtually identical. However this should not be the case (again see Hunt et al., 2016). It is quite hard for me to judge the spatial signature of connectivity as, unfortunately, in the pdf version the figure quality is too low resolution to read the region names on the axes. However, these matrices do not look correct to me and I would urge the authors to find out why.

We would like to thank you for noticing this abnormality, and encouraging us to troubleshoot our MEG connectivity scripts. Indeed, we did notice an error in the analysis that affected the estimation of connectivity in the β band. When switching from α to β band we initially neglected to change the leadfield file during the construction of the virtual timeseries. This lead to an estimation of connectivity in the β band using the spatial filter optimized for the α band. This was an oversight, and has been corrected in the revised analyses. In addition, we understand that the original figures were difficult to read in the auto-generated PDF, and have included the original picture files for all adjacency matrices as a separate zip attachment.

After correcting this mistake, we re-ran the β connectivity analysis using the same method as in the previous revision. When looking at the adjacency matrices from the downsampled power envelopes, we noticed that the α and β adjacency matrices still looked very similar (See Author response image 6).

Author response image 6
Adjacency matrices constructed from downsampled timeseries for the α and β bands.
https://doi.org/10.7554/eLife.23608.027

Therefore, as recommended in your decision letter we re-ran the connectivity analyses for both the α and β bands using the non-downsampled, Hilbert-transformed, orthogonalized timeseries. (Author response image 7) This analysis yielded adjacency matrices for the α and β band that looked more distinct than the analysis using the downsampled timeseries. However, the resulting correlation coefficients were much smaller from the non-downsampled timeseries.

Author response image 7
Adjacency matrices constructed from non-downsampled timeseries for the α and β bands.
https://doi.org/10.7554/eLife.23608.028

One possibility for this reduction in the magnitude of the correlations is that some of the structure in the initial adjacency matrices resulted from variability across trials, and this variability represented a smaller proportion of the overall variability when the analysis was based on the non-downsampled timeseries. This might also explain why the adjacency matrices for the α and β bands looked so similar when constructed from the downsampled timeseries. To test this possibility, we z-scored the individual trial data before concatenating the trial timeseries and conducting the correlations. If across trial variability is sufficient to explain the larger correlation coefficients and similarity between the adjacency matrices, these effects should be eliminated by the z-scoring. This does not appear to be the case. Based on the resulting adjacency matrices, the downsampled data yield higher correlation coefficients, and the adjacency matrices from the downsampled data look similar between the α and β bands. See Author response image 8.

Author response image 8
Adjacency matrices constructed from downsampled timeseries that have been converted to z-scores for the α and β bands.
https://doi.org/10.7554/eLife.23608.029

Furthermore, the adjacency matrices from the non-downsampled data seem to lose their structure (Author response image 9).

Author response image 9
Adjacency matrices constructed from non-downsampled timeseries that have been converted to z-scores for the α and β bands.
https://doi.org/10.7554/eLife.23608.030

Together, these results suggest that the best approach is probably to use the non-downsampled, Hilbert-transformed, orthogonalized timeseries to create the adjacency matrices. However, it should be noted that this approach did not yield any significant differences in connectivity as a function of threat (See the high-resolution images included in the attached zip file). Furthermore, although the other approaches seemed to be sub-optimal, we tested whether there were any differences in connectivity as a function of threat with these as well, and were unable to identify any significant differences. In summary, we feel we have conducted a thorough analysis of the connectivity data, using methods recommended by both yourself and the reviewers, and were unable to identify any changes in connectivity as a function of threat.

Although we are happy to conduct additional analyses if requested, we feel that these analyses will likely be unfruitful, and only delay a decision on the manuscript. We understand that you are enthusiastic about the manuscript, and hope that you find it acceptable in its current form.

In their response, the authors have suggested that they use blocks of only 2 seconds of envelope data. However they also apply downsampling prior to connectivity estimation. This to me seems silly – why apply the downsampling if you have such short blocks. It is well established that connectivity works without the downsampling so perhaps removing that step might improve the adjacency matrices. I also don't understand why such short data setments were used when in fact the bliocks are quite long. Why not just try using the whole block?

We agree that downsampling the timeseries prior to conducting the correlation analyses was not the best approach, and have addressed this limitation (See response above).

We also agree that analyzing the entire block (rather than a series of 2 s periods) might yield more power to detect a potential difference in connectivity as a function of threat. However, we initially chose to segment the data into 2 s intervals corresponding to the pre-startle period for 2 reasons. 1) We wanted to characterize the effect of threat on pre-stimulus oscillations during the otherwise well-studied threat of shock paradigm. 2) We wanted to minimize the effect of the blink artifact on our power estimates. This is especially important in our study because startle-elicited blinks (and thus the blink artifact) are larger under threat. Although it could be argued that the adaptive beamformer spatial filter should remove noise associated with blinks from the data, we prefer a more conservative approach, and therefore chose to include only verifiably clean data in our analysis.

Again, we understand that using the data from the entire block may yield more power for the connectivity analysis, this would require re-analyzing and cleaning the data from scratch. Not only would this be time consuming, but it would also require re-running the data-cleaning ICA and the beamformer analysis, which would bring the connectivity data out of sync with the power data reported in the manuscript. If there were some indication that our current null connectivity results may be due to type 2 error, it could be argued that this re-analysis would be justified. However, given that several previous analyses yielded no effect of threat on connectivity in the MEG data, we do not feel that our results are due to type 2 error, and feel that this re-analysis of the whole-block data would not be fruitful.

Although we are happy to undertake the analysis if requested, it would require a substantial amount of time, and would significantly delay decision on the manuscript.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

I would like to thank the authors for their attempt at addressing my concerns. I'm pleased that the authors were able to spot errors in their pipeline (use of incorrect spatial filters) although it does appear that the major difference came from the down-sampling, as suggested in my previous review. Downsampling over such a small time window is obviously incorrect so these analyses should now be discarded.

[…]

In the light of these failed attempts, I suggest that a paragraph is added to the Discussion stating that MEG connectivity analysis would have been a useful means to probe these data; however the study design was poorly set up to make such analyses work.

We would like to thank the reviewer for this suggestion, and we agree that this addition is necessary. In order address this concern we have added the following paragraph to the functional connectivity section of the Discussion, and edited the limitations section to be consistent with this addition.

An obvious parallel to the fMRI connectivity analysis would be to conduct a similar whole brain connectivity analysis of the MEG data. Not only would corroborating evidence from an independent imaging modality strengthen the fMRI connectivity results, but results from MEG specifically would allow for a frequency specific analysis of the effects of threat on functional connectivity (Brookes et al., 2011; Hillebrand et al., 2012). However, the current study was not optimized to reliably detect differences in MEG connectivity. In the current study we included the white noise probes as a way to obtain a quantitative measure of anxiety throughout the safe and threat blocks (Grillon et al., 1999). These white noise probes trigger the acoustic startle reflex, which varies as a function of an individual’s current level of anxiety (Grillon, 2008). Unfortunately, these reflexive blinks also inject an artifact into the MEG signal, and because the magnitude of these blinks differs across conditions, the blink artifact also differs across conditions. Although adaptive beamforming can theoretically remove the artifact induced by the blink response (Van Veen et al., 1997), the only way to ensure that the blink artifact does not differentially affect the MEG signal is to remove the contaminated time periods from the analysis, or remove the startle probes from the design at the outset. In the current study, we chose to address this limitation by extracting two-second time windows prior to each startle presentation to minimize the effect of the blink artifact on our power estimates. However, it has been shown that reliability of the MEG connectivity estimates increases as the duration of the recording increases, and durations of ~10 min or greater may be needed to maximize reliability (Liuzzi et al., 2016). Therefore, using such short intervals did not allow for the ability to obtain reliable estimates of MEG connectivity. Future studies should be conducted to address this limitation. In addition, it will be important to use appropriate correction methods to account for signal leakage between source space regions (Colclough, Brookes, Smith and Woolrich, 2015), and verify that the resulting connectivity estimates match previously published work (i.e. strong α connectivity in occipital regions and strong β connectivity linking the parietal cortex with other frontal, temporal, and occipital regions; (Hunt et al., 2016)).

https://doi.org/10.7554/eLife.23608.032

Article and author information

Author details

  1. Nicholas L Balderston

    Section on Neurobiology of Fear and Anxiety, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    NLB, Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    nicholas.balderston@nih.gov
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8565-1544
  2. Elizabeth Hale

    Section on Neurobiology of Fear and Anxiety, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    EH, Investigation, Methodology, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
  3. Abigail Hsiung

    Section on Neurobiology of Fear and Anxiety, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    AH, Investigation, Methodology, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
  4. Salvatore Torrisi

    Section on Neurobiology of Fear and Anxiety, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    ST, Conceptualization, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
  5. Tom Holroyd

    MEG Core Facility, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    TH, Conceptualization, Software, Methodology, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
  6. Frederick W Carver

    MEG Core Facility, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    FWC, Conceptualization, Software, Methodology, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
  7. Richard Coppola

    MEG Core Facility, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    RC, Conceptualization, Resources, Software, Methodology, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
  8. Monique Ernst

    Section on Neurobiology of Fear and Anxiety, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    ME, Conceptualization, Supervision, Methodology, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
  9. Christian Grillon

    Section on Neurobiology of Fear and Anxiety, National Institute of Mental Health, National Institutes of Health, Bethesda, United States
    Contribution
    CG, Conceptualization, Supervision, Investigation, Methodology, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.

Funding

National Institute of Mental Health (ZIAMH002798)

  • Christian Grillon

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors report no conflicts of interest. We would like to thank the NIMH Section on Instrumentation for fabricating the custom acoustic startle stimulus devices. Financial support of this study was provided by the Intramural Research Program of the National Institute of Mental Health, ZIAMH002798 (ClinicalTrial.gov Identifier: NCT00047853: Protocol ID 02 M-0321).

Ethics

Human subjects: All participants gave written informed consent approved by the National Institute of Mental Health (NIMH) Combined Neuroscience Institutional Review Board and received compensation for participating.

Reviewing Editor

  1. Matthew J Brookes, University of Nottingham, United Kingdom

Publication history

  1. Received: November 28, 2016
  2. Accepted: May 29, 2017
  3. Accepted Manuscript published: May 30, 2017 (version 1)
  4. Version of Record published: June 20, 2017 (version 2)

Copyright

This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

Metrics

  • 1,377
    Page views
  • 154
    Downloads
  • 11
    Citations

Article citation count generated by polling the highest count across the following sources: Scopus, Crossref, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Neuroscience
    María Fernanda López-Gutiérrez et al.
    Research Article

    Previous studies have related pair bonding in Microtus ochrogaster, the prairie vole, with plastic changes in several brain regions. However, the interactions between these socially-relevant regions have yet to be described. In this study, we used resting state magnetic resonance imaging to explore bonding behaviors and functional connectivity of brain regions previously associated with pair bonding. Thirty-two male and female prairie voles were scanned at baseline, 24h and 2 weeks after the onset of cohabitation By using network based statistics, we identified that the functional connectivity of a cortico-striatal network predicted the onset of affiliative behavior, while another predicted the amount of social interaction during a partner preference test. Furthermore, a network with significant changes in time was revealed, also showing associations with the level of partner preference. Overall, our findings revealed the association between network-level functional connectivity changes and social bonding.

    1. Developmental Biology
    2. Neuroscience
    Alessia Caramello et al.
    Research Article Updated

    During embryonic development, radial glial cells give rise to neurons, then to astrocytes following the gliogenic switch. Timely regulation of the switch, operated by several transcription factors, is fundamental for allowing coordinated interactions between neurons and glia. We deleted the gene for one such factor, SOX9, early during mouse brain development and observed a significantly compromised dentate gyrus (DG). We dissected the origin of the defect, targeting embryonic Sox9 deletion to either the DG neuronal progenitor domain or the adjacent cortical hem (CH). We identified in the latter previously uncharacterized ALDH1L1+ astrocytic progenitors, which form a fimbrial-specific glial scaffold necessary for neuronal progenitor migration toward the developing DG. Our results highlight an early crucial role of SOX9 for DG development through regulation of astroglial potential acquisition in the CH. Moreover, we illustrate how formation of a local network, amidst astrocytic and neuronal progenitors originating from adjacent domains, underlays brain morphogenesis.