The relation between neuronal firing, local field potentials and hemodynamic activity in the human amygdala in response to aversive dynamic visual stimuli

The amygdala is a central part of networks of brain regions underlying perception and cognition, in particular related to processing of emotionally salient stimuli. Invasive electrophysiological and hemodynamic measurements are commonly used to evaluate functions of the human amygdala, but a comprehensive understanding of their relation is still lacking. Here, we aimed at investigating the link between fast and slow frequency amygdalar oscillations, neuronal firing and hemodynamic responses. To this aim, we recorded intracranial electroencephalography (iEEG), hemodynamic responses and single neuron activity from the amygdala of patients with epilepsy. Patients were presented with dynamic visual sequences of fearful faces (aversive condition), interleaved with sequences of neutral landscapes (neutral condition). Comparing responses to aversive versus neutral stimuli across participants, we observed enhanced high gamma power (HGP, >60 Hz) during the first 2 s of aversive sequence viewing, and reduced delta power (1-4 Hz) lasting up to 18 s. In 5 participants with implanted microwires, neuronal firing rates were enhanced following aversive stimuli, and exhibited positive correlation with HGP and hemodynamic responses. Our results show that high gamma power, neuronal firing and BOLD responses from the human amygdala are co-modulated. Our findings provide, for the first time, a comprehensive investigation of amygdalar responses to aversive stimuli, ranging from single-neuron spikes to local field potentials and hemodynamic responses.


Introduction
The human amygdala is a central hub of several networks underlying human perception and cognition. Various tasks, such as novelty detection (Balderston et al., 2013), perception of faces (Adolphs and Spezio, 2006) and emotions (Wang et al., 2014), or aversive learning (Janak and Tye, 2015), rely on neural responses within the amygdala. Although neural oscillations in the amygdala have been well studied in rodents, in humans the vast majority of evidence stems from hemodynamic studies (Janak and Tye, 2015).
In particular, previous studies have highlighted amygdalar contributions in processing static (Said et al., 2011) and dynamic information related to faces (Schultz et al., 2012;Zinchenko et al., 2018). Dynamic faces contain rich information over and beyond static ones, as they have higher ecological value stemming from motion information (Dobs et al., 2018). Presentation of dynamic faces has been shown to elicit strong electrophysiological responses, at the scalp electroencephalography (EEG) level (Recio et al., 2011), and to facilitate emotional recognition (Calvo et al., 2016). Within the face perception network, the human amygdala is an important node (Duchaine and Yovel, 2015), but its role in dynamic face perception has been mainly investigated by means of Blood Oxygen Level Dependent (BOLD) responses (Krumhuber et al., 2013;Pitcher et al., 2019;Schacher et al., 2006), leaving its electrophysiological content severely under-explored (Zheng et al., 2017).
Indeed, despite recent progress in non-invasive techniques, such as magnetoencephalography (Tzovara et al., 2019;Dumas et al., 2013), oscillatory activity from the human amygdala can be mainly recorded using invasive techniques, such as intracranial electroencephalography (iEEG) in patients with refractory epilepsy who are undergoing pre-surgical monitoring. iEEG recordings have an astonishing potential to assess local activity from ensembles of hundreds of thousands of neurons, and infer from those slow oscillations (e.g. DELTA band 1-4 Hz), as well as fast and synchronized activity, in the so-called high gamma power (HGP > 60 Hz) (Lachaux et al., 2012). In some rare cases, iEEG recordings can even give access to neuronal firing of single neurons (Fried et al., 1997;Parvizi and Kastner, 2018). Although these different types of data provide complementary pieces of information about the way that the human brain processes sensory stimuli, attempts to combine them and study their concomitant occurrence remain rare.
Past studies assessing multi-modal measurements of brain activity have contrasted BOLD responses and neural rhythms, mainly recorded with scalp EEG (Murta et al., 2015) reporting correlations over wide-spread cortical regions during rest (Laufs et al., 2003) and in particular in relation to resting alpha rhythms (Benar et al., 2007;de Munck et al., 2007;Difrancesco et al., 2008;Goncalves et al., 2006;Mulert et al., 2008). Higher specificity in the links between BOLD and electrophysiology has been achieved by contrasting fMRI with iEEG (Privman et al., 2007). Multiple studies have reported positive correlations between high gamma power (generally > 70 Hz), recorded with invasive EEG, and hemodynamic responses, either recorded simultaneously, on in separate sessions (Conner et al., 2011;Hermes et al., 2012;Khursheed et al., 2011;Kunii et al., 2013;Murta et al., 2017). By contrast, BOLD measurements and the power of low frequency oscillations exhibit a negative correlation (theta (Murta et al., 2017), beta (Conner et al., 2011)). Moreover, high frequency activity has been shown to correlate not only with BOLD responses, but also neuronal spiking (Mukamel et al., 2005;Ray et al., 2008) and mainly with neurons from superficial cortical layers (Leszczy nski et al., 2019). Overall, links between hemodynamic and electrophysiological activity have been reported for wide-spread cortical areas, and a variety of experimental tasks. Nevertheless, it still remains unknown whether similar patterns of co-modulations between BOLD and electrophysiological responses can be also observed in subcortical brain regions, such as the amygdala. Given its distinct properties and unique role in emotions and survival circuits, here we aimed at providing a fine-grained characterization of amygdalar multi-modal measurements.
We hypothesized that high frequency activity, neuronal firing, and BOLD activity show similar patterns of responses to emotional stimuli, and strong correlations, while slow oscillations show the opposite behavior, similar to what has been previously reported for the cortex. Additionally, we hypothesized that these co-modulations would be mostly apparent during conditions that have been shown to highly (A) Video still of dynamic AVERSIVE (red) and NEUTRAL (blue) sequence. Over the group of participants, delta power (DELTA, 1-4 Hz) decreased (B) and high gamma power (HGP, >60 Hz) increased during the first 2 s (C, as in (Zheng et al., 2017)). Z-score AE SEM shown as shading around the mean trace and locked to stimulus onset. Green bars represent significant differences between AVERSIVE and NEUTRAL (cluster analysis p < 0.05). (D) The modulation index (MI) for iEEG was defined as the difference in time-averaged iEEG power between the AVERSIVE and NEUTRAL conditions where iEEG power is the amplitude envelope in the DELTA or HGP during the interval where the two curves are significantly separated. Across all amygdalae (black markers within the seizure onset zone SOZ, gray circles outside SOZ), DELTA and HGP response were negatively correlated (r ¼ 0.67, p ¼ 0.01). (E) HGP but not DELTA distinguishes between amygdalae within SOZ and outside SOZ (p ¼ 0.019, not corrected). (F) The HGP response exceeds the DELTA response in specificity when indicating the amygdala's dysfunction within SOZ (p~0.2 HGP 5/5 DELTA 1/5, chi2-test).
engage the amygdala, such as viewing of emotional dynamic videos.
We employed a task that includes presentation of salient visual stimuli in a naturalistic way, in contrast to most human single-neuron studies that have only used static stimuli (Canli et al., 2002;Rutishauser et al., 2011;Wang et al., 2014). We aimed at providing direct experimental evidence on the interactions among iEEG, HGP, neuronal firing and BOLD activity within the human amygdala. Our paradigm has been previously shown to elicit strong BOLD (Schacher et al., 2006;Steiger et al., 2017) and HGP (Zheng et al., 2017) responses in the amygdala. However, there has never been a direct comparison between these two measures, and their very important counterparts (i.e. slow oscillations and neuronal firing).

iEEG power in delta and high gamma across participants
We recorded intracranial EEG in 9 epilepsy patients while they watched video sequences of dynamic fearful faces (AVERSIVE condition) interleaved with sequences of landscapes (NEUTRAL condition) (Fig. 1  A). Across the group of participants, we recorded iEEG from 7 healthy amygdalae (non Seizure Onset Zone -nSOZ) and 7 amygdalae within the SOZ (Table 1). In nSOZ amygdalae, DELTA band power was consistently reduced throughout the duration of the AVERSIVE sequence (Fig. 1B). The HGP was enhanced in response to aversive stimuli during the first 2 s of the AVERSIVE sequence presentation (Fig. 1C, unpaired t-tests, corrected for multiple comparisons using random permutations, p < 0.05), in accordance to an earlier report using the same experimental manipulation but on different patients (Zheng et al., 2017). Thus, the aversive visual stimulation modulated slow and high frequency amygdalar activity consistently across the group of participants. Since no effect was observed in the other spectral bands (see Methods 7.8), in the following we focus on DELTA power and HGP.
To quantify the direction of response in different frequency bands, we defined the modulation index (MI) as the difference in iEEG power between the AVERSIVE and NEUTRAL conditions in the temporal window of task-dependent modulation, which occurred at stimulus onset in HGP and was sustained throughout the trial in DELTA ( Table 2). The DELTA-MI and HGP-MI were negatively correlated across all amygdalae (Spearman's r ¼ À0.67, p ¼ 0.01, Fig. 1D). Specifically for amygdalae not included in the SOZ, we observed that the stronger the HGP enhancement, the stronger the reduction in DELTA power (Table 2). Therefore, we next asked whether this inverse relationship between DELTA power and HGP ( Fig. 1 E) might be affected by the epileptogenicity of the amygdalar tissue. DELTA power was reduced in 12 of 14 amygdalae (negative MI in Table 2), irrespective of whether the amygdalae were within the SOZ or not. This robust response suggests that residual neural processing might take place within the SOZ, but, at the same time, changes in DELTA power were irrelevant for the identification of epileptogenic tissue. To the contrary, the modulation of HGP was more specific for the epileptogenicity of the amygdalar tissue. All amygdalae outside SOZ showed a positive HGP-MI in response to AVERSIVE videos, while this response was absent in amygdalae within the SOZ (p ¼ 0.009). With sensitivity (7/7 ¼ 100% CI [66-100]%), high specificity (6/7 ¼ 86% CI [42 99]%) and accuracy 93% (13/14, CI [66-99]%), HGP enhancement indicated amygdalar functional integrity.
2.2. The iEEG and the BOLD are co-modulated As hemodynamic measures are widespread in humans, it is desirable to base them on electrophysiological activity. Our video material was initially designed to assess the functionality of amygdalar structures based on BOLD responses in pre-operative fMRI (Schacher et al., 2006;Steiger et al., 2017). Among the 9 patients of the current study, 7 had also viewed the video clips during fMRI several months previous to the electrode implantation (a subset of participants included in (Schacher et al., 2006)). In 9/14 amygdalae of these participants, we compared response patterns in BOLD and iEEG activity. Fig. 2 shows BOLD and iEEG responses in the left amygdala of Participant 4, in two separate sessions of presentation of the same sequence of videos.
We then examined whether hemodynamic and electrophysiological responses showed any commonalities. Interestingly, we found a dissociation in the correlation of BOLD responses to slow and high frequency iEEG oscillations: BOLD response and DELTA power were negatively correlated during presentation of video sequence (correlation of BOLD-DELTA time-courses r ¼ À0.48 p < 0.01), while BOLD and HGP were positively correlated (correlation BOLD-HGP Spearman's r ¼ 0.41 p < 0.01) for the participant displayed in Fig. 2.
A similar pattern was also observed across the group of participants. DELTA power and BOLD responses showed a negative correlation in 4/5 nSOZ amygdalae, while the correlation of HGP and BOLD was positive in the 3/5 the nSOZ amygdalae and 1/4 SOZ amygdalae (Table 2, Fig. 2C). Interestingly, the two nSOZ amygdalae that did not show a significant correlation between HGB and BOLD, had a remarkably low MI for HGP (MI < 0.05, Table 2). In summary, the positive correlation between HGP and BOLD points to a common origin in the underlying physiology of these two signals, and an interesting dissociation with respect to slow oscillations.

Neuronal firing and oscillatory activity are co-modulated
To increase the spatial specificity of our analysis, we recorded activity from microwires protruding from the tip of the amygdala macroelectrodes ( Fig. 3A). We identified 50 putative clusters of neurons (Niediek et al., 2016) from 9 amygdalae and included those with firing rate >1 Hz (N ¼ 41 neurons) in further analysis (see Supplemental Fig. S2 for the waveforms and inter-spike intervals for neurons identified for each patient).
First, we investigated whether viewing of video sequences might modulate neuronal firing. Fig. 3A shows the raster plot of a representative neuron and its firing rate averaged over trials. All 41 clusters of neurons increased their firing rate during the AVERSIVE condition (Firing Rate Variation, FRV ¼ 0.26 AE 0.15). Contrary to our expectation of an immediate spiking onset, the neurons increased their firing gradually and sustained enhanced firing throughout the AVERSIVE video sequences.
Second, we asked how the neuronal firing contributed to HGP. We computed the correlation between neuronal firing and HGP, recorded from macro-contacts (still referred in the following as HGP), and locally at the level of the microwires (μHGP). We observed a high covariation of neuronal firing with HGP from both micro-wires (μHGP, Fig. 3B upper inset) and macro contacts (HGP, Fig. 3B lower inset). Across all neurons, the covariation of neuronal firing with HGP was not significantly different between the two experimental conditions (Spearman's r ¼ 0.17 Table 1 Patients and recordings. For each patient, we report age, gender, pathology, hemisphere of the amygdala implanted with iEEG (L ¼ left, R ¼ right, SOZ ¼ seizure onset zone, marked with *), the presence of significant BOLD response following fMRI recordingaccording to Schacher et al. (2006) -and multi-unit activity following micro-wire implantation. . This is evidence that μHGP has sufficient spatial precision to describe the modulation of firing in local neuronal assemblies. Third, we investigated the synchronization among neurons. Fig. 3C shows the firing rate of 6 neuronal clusters recorded from different microwires in the same amygdala of one exemplar participant. For each neuron across the group of participants, we then computed the spike synchronization index (SSI, averaged pairwise correlation of the firing rate of this neuron with that of all neighboring neurons), which is a commonly used metric to quantify the co-variation of firing in a neuronal assembly. We found a strong positive correlation across all participants between the SSI and FRV (Spearman's r ¼ 0.78, p < 1e-7). Together with the increased neuronal firing in the AVERSIVE condition (positive FRV), this would suggest that the correlation is driven by neuronal assembly synchronized firing during the AVERSIVE, rather than the NEUTRAL condition.

Table 2
The AVERSIVE condition modulates iEEG power and BOLD. The modulation index (MI) for iEEG is defined as the difference in time-averaged iEEG power between the AVERSIVE and NEUTRAL conditions where iEEG power is the amplitude envelope in the delta or high gamma band during the interval where the two curves are significantly separated. For all 6/7 SOZ amygdalae, the MI of HGP is negative, which indicates that no amygdalar response was elicited by AVERSIVE video clips. The BOLD was negatively correlated to in DELTA (4/5 with p < 0.05) and positively correlated to HGP (5/5 with p < 0.05). p-values are corrected by FDR.  To understand to which extent the variation in neuronal firing reflects the neural dynamics of amygdalar oscillations at a macroscale, we compared the firing rates of isolated clusters to μHGP, HGP and BOLD during the video clip presentation. Fig. 3D shows an example from an exemplar participant, where the firing rate, μHGP, HGP, and BOLD covary in the timescale of the task (smoothing 15 s window). Across all neurons and participants, the ongoing firing rate positively correlated with μHGP (Spearman's r ¼ 0.27 AE 0.15), HGP (Spearman's r ¼ 0.19 AE 0.17) and BOLD (Spearman's r ¼ 0.23 AE 0.20) (Fig. 3 D inset). Across all recording modalities, the correlation of the firing rate and other measures was lowest with HGP (Wilcoxon's paired test, p < 1e-2). Importantly, this covariation reflects common correlates of amygdalar activity across recording modalities, ranging from the neuronal firing, to high frequency micro-and macro-scopic oscillations, up to hemodynamic responses.

Neuronal firing, HGP and BOLD are co-modulated
The variability in neural firing driven by AVERSIVE stimulation (high FRV) correlated with the presence of tightly synchronized neuronal assemblies (high SSI) as depicted in Fig. 4A. We asked whether this microscale activity is also represented in the macroscopic measurements of HGP and BOLD. To this aim, we computed a multivariate linear regression model, to examine how neuronal features as FRV and SSI are contributing to describe the association between firing rate and HGP, or firing rate and BOLD. We found that the co-modulation between neuronal firing and HGP is described by both FVR and SSI (Fig. 4B, r ¼ 0.59, p < 1e-4), while the co-modulation between neuronal firing and BOLD is described by FRV (Fig. 4C, Spearman's r ¼ 0.69, p < 1e-4), but not SSI. Therefore, local synchronization (SSI) is a prerequisite for a correlation between neuronal firing and iEEG, but not for correlation between neuronal firing and BOLD. By contrast, the task-dependent modulation of neuronal firing contributes to the link between neuronal firing and HGP as well as neuronal firing and BOLD.
Finally, we observed a correspondence across all three recording modalities neuronal firing, macroscopic HGP and BOLD. Fig. 4D displays the covariation of neuronal firing with BOLD (x-axis) and HGP (y-axis) (Spearman's r ¼ 0.72, p < 1e-3). Thus, neuronal firing that covaries with HGP, most likely also covaries with the hemodynamic response. We thereby provide a unified picture of how dynamic visual stimulation elicits consistent neurophysiological responses in our three recording modalities: macroscopic high gamma power, BOLD, and neuronal firing.

Discussion
We examined commonalities across multimodal representations of electrophysiological and hemodynamic responses in the human amygdala during dynamic visual stimulation.

iEEG waves in delta and high gamma
While previous studies investigating the role of amygdala in processing fearful faces have focused on static images (Murray et al., 2014), here we used a paradigm with higher ecologic value that presents faces dynamically, in video clips. We observed a decrease in DELTA power and an increase in HGP during viewing of AVERSIVE compared to NEUTRAL videos. While the HGP increase is consistent with an earlier report (Zheng et al., 2017), we additionally show that this HGP increase is accompanied by a persistent decrease in DELTA. The enhanced HGP is most evident at the onset of the trial within the first 2 s of the onset of the video clip. An early amygdala response is barely surprising, as previous studies have shown that amygdala responses to images of fearful faces can occur as fast as <100 ms after the image onset (Krolak-Salmon et al., 2004). Interestingly, early components have been linked to explicit attention (Pourtois et al., 2010), while later components (600-800 ms) have been associated with both implicit and explicit attention to emotional faces (Murray et al., 2014). In our case, the suppression of delta amplitude lasted for most of the AVERSIVE sequence, suggesting that faster dynamics, peaking at the onset of the sequence, account for a switch in activity level, while lower frequencies, which were sustained over the whole sequence might be responsible for its continued visual perception. In the context of our study, DELTA power may rather be associated to inhibitory regulation of the sensory input to optimize local computation (Harmony, 2013).

iEEG and BOLD
As the most consistent result across all participants, the HGP increase and DELTA power reduction in response to AVERSIVE videos appear as the electrophysiological counterpart of the observed increase in BOLD response. In search for an electrophysiological counterpart for BOLD responses, animal studies have reported a theta (1-8 Hz) decrease concurrent with high gamma power increase in regions of BOLD activation, which mirrors our findings (Burke et al., 2014;Buzsaki et al., 2012;Ekstrom, 2010;Logothetis and Wandell, 2004;Niessing et al., 2005). In humans, the co-localization of electrophysiological and BOLD response shows an inversion of correlations between high gamma and low frequency fluctuations (Manning et al., 2009;Mukamel et al., 2005;Nir et al., 2007;Scheeringa et al., 2011) in visual and auditory cortex. Interestingly, compared to theta oscillations in rodents, human theta oscillations are expected in lower frequency bands, generally in 1-4 Hz (Jacobs, 2014). Previous studies have shown temporal correlations between a delta power reduction and BOLD responses during sleep in anterior cingulate and orbitofrontal cortex (Hofle et al., 1997), where delta is generated by the mesotelencephalic system (Murphy et al., 2009). During wakefulness, delta oscillations have been found to negatively correlate with BOLD responses in fronto-parietal cortical structures (Marawar et al., 2017) where they have been linked to motivation (Knyazev, 2007). Consistent with previous studies in humans, we show here a negative correlation between BOLD signal and DELTA. We may speculate that a decrease in DELTA indicates the gating that facilitates the enhanced HGP connectivity between the amygdala and the hippocampus (Zheng et al., 2017).

Neuronal firing, macroscopic oscillations and BOLD
Neuronal firing, HGP and BOLD responses were first compared in the macaque visual cortex (Logothetis et al., 2001). In humans, multimodal analysis was performed in several cortical areas of the brain, from the visual to the associative cortex (Ojemann et al., 2013) or in the hippocampus (Ekstrom, 2010). To the best of our knowledge, we provide here the first study evaluating common correlates in neuronal firing, iEEG and BOLD in the human amygdala. We observed that the firing rate of isolated neurons in the amygdala increases during the viewing of fearful faces. It was previously shown that single neurons in the amygdala competitively respond to face features (Rutishauser et al., 2011) and encode a perceptual judgement of emotional faces (Wang et al., 2014). While emotional judgment might be subject to inter-individual differences, amygdala activation to fearful expression is highly reproducible across individuals (Canli et al., 2002). Consistent with previous studies in the auditory cortex, the correlation between HGP and neuronal firing is proportional to the degree of interneuronal correlation across neurons , while HGP also explains the variability of BOLD responses (Conner et al., 2011;Kujala et al., 2014). Considering that pharmacologically induced depression of efferent neuronal firing in the macaque visual cortex does not affect BOLD and the local field potential (Rauch et al., 2008), the correlation between neuronal firing, HGP and BOLD might reflect the processing of afferent information. While our experimental design does not target the distinction between input and output of amygdalar processing, we could identify a tight connection across modalities (Burns et al., 2010) over the coarse timescale of several seconds.

Limitations
We need to mention some limitations of the present study. First, we cannot assess whether the observed effects of the amygdala are limited to fearful faces, since our control condition is videos of landscapes and not neutral faces. This is because we targeted the activation of amygdala using a complex ecological visual input that was optimized for functional clinical evaluation (Schacher et al., 2006;Steiger et al., 2017). Second, we could not assess whether neuronal firing carries any clinical value for detecting epileptogenic tissue, as we had no microwire recordings from SOZ amygdalae. Third, we compared BOLD activations with iEEG that were recorded in the same participant but not simultaneously. While simultaneous recordings might have revealed additional fine-grained correlations across BOLD and iEEG signals, our approach is still valid, because the test-retest reliability of EEG spectra while performing a cognitive task has been shown to be high even after more than a year (N€ apflin et al., 2008).

Clinical relevance
While epilepsy is increasingly regarded as a disease affecting a network-wide organization in the brain (Lagarde et al., 2018;Steiger et al., 2017), we show here that we can address the functionality of a specific area, reconciling different spatiotemporal scales and modalities of investigation. The functional evaluation prior to epilepsy surgery is mostly based on non-invasive assessment, i.e. the hemodynamic response in tailored fMRI tasks. Nevertheless, each recording modality has its limitations and might be misguiding during the clinical evaluation (Muthukumaraswamy and Singh, 2009). Here, we have shown that HGP in the iEEG may serve as a cognitive marker to indicate the functional integrity of the human amygdala, which in turn could be used to support the clinical treatment.

Participants and electrodes
We included 9 patients with drug-resistant mesial temporal lobe epilepsy. All patients were implanted with intracranial electrodes for presurgical monitoring in order to localize their epileptogenic zone ( Table 1). The placement of electrodes was exclusively guided by clinical criteria. All patients were recruited at Schweizerische Epilepsie-Klinik, Zurich. Intracranial depth electrodes were implanted stereotactically at locations planned according to the results of the non-invasive presurgical work-up. Electrode placement was guided exclusively by clinical needs. Macro electrodes had 1.3 mm diameter, 8 contacts of 1.6 mm length, spacing between contacts centers 5 mm, and contained a bushel of microwires that protruded from their tip to record spiking of single neurons (AD-Tech, www.adtechmedical.com).

Ethics statement
The study was approved upfront by the institutional ethics review board (Kantonale Ethikkommission PB-2016-02055) and participants gave written informed consent before the measurements.

Depth electrode localization
Electrode position was verified in participants' native space by the neurosurgeon (LS) after merging preoperative MR with postimplantation CT images of each individual patient (iPlan Stereotaxy 3.0, Brainlab, Germany).

Stimulus material
We used a paradigm comprising of a series of dynamic videos, which has been already validated in previous clinical investigations (Schacher et al., 2006). The videos were all silent and consisted of dynamic fearful faces and dynamic neutral landscapes, presented in an alternating order, in a block design. The paradigm included eight blocks of 75 short video clips (2-3 s) of fearful faces and nine blocks of 72 short video clips (2-3 s) of neutral landscapes. Each block lasted 24 s in total, and contained short video clips without any intermission between consecutive videos. Video clips of fearful faces were extracted from thriller and horror movies and contained faces of actors showing fear, without being violent or aggressive. Video clips of neutral landscapes were chosen as a control condition, and were matched to the duration of the fearful faces videos (2-3 s). They included domestic landscapes which are posited to have a low emotional content and visual properties comparable to the emotional videos (Schacher et al., 2006). All videos were only included once. A panel of psychologists had evaluated the stimuli to ensure that they are suitable for the patients and that they do not include any episodes of violence or aggression (Schacher et al., 2006). In particular, we started with a set of 120 videos of fearful faces and reduced that to 72, by excluding videos where: (a) the actor's face was not continuously visible (b) fear was not clearly recognized on the actor's face (c) no other emotion was displayed (e.g. anger/surprise) and (d) the display of fear was intense.
During electrophysiological recordings the videos were presented to the patients via a laptop screen, while during the fMRI scan they were presented through a tilted overhead mirror. In both cases, patients were instructed to pay attention to the videos and focus on the eyes of the actors during the clips containing faces. For the electrophysiological recordings, blocks were separated by a repeated baseline of 2 s taken from a NEUTRAL condition.

fMRI data acquisition
Structural and functional scans were acquired on a 3.0-T Philips scanner (Philips Medical Systems, Best, The Netherlands) with a 32-channel head coil. Foam padding was fixed around the participant's head to reduce head motion. Anatomical scans were acquired with a T1-weighted 3D magnetization-prepared rapid gradient echo (MPRAGE) sequence with the following parameters: 176 sagittal slices, thickness ¼ 1 mm, skip ¼ 0 mm, repetition time (TR) ¼ 8.1 ms, echo time (TE) ¼ 3.7 ms, flip angle ¼ 8 , field of view (FOV) ¼ 240 Â 240 mm 2 , isotropic voxel size ¼ 1 mm 3 . Functional data was obtained with a gradient echo planar imaging (EPI) sequence. To maximize signal to noise ratio, 18 coronal slices were orientated orthogonally to the axis of the hippocampi. Slices covered temporal and frontal lobes with TR ¼ 1500 ms, TE ¼ 35 ms, flip angle ¼ 75 , FOV ¼ 220 Â 220 mm 2 , voxel size ¼ 2.75 Â 2.75 Â 4 mm.

fMRI analysis
BOLD responses were extracted using BrainVoyager QX version 4.8 (BrainInnovation, Maastricht, The Netherlands) for each participant separately as part of their clinical pre-surgical evaluation. Patients were presented with the same set of stimuli in the MRI scanner and were instructed to pay attention to the videos and focus on the eyes of the actors during aversive clips. They were asked to rate their emotional involvement on a continuous scale from 0 to 10 after the scan.
Preprocessing of the functional data involved slice time correction to account for the interleaved slice acquisition, rigid-body motion correction to the first functional image, linear trend removal by temporal fast Fourier-transform based high-pass filtering and spatial smoothing with a 4 mm full width at half maximum Gaussian kernel. BOLD responses were extracted from volumes of interest (VOIs), defined at the single-patient level, and ranging between 30 and 100 voxels (3.9 Â 3 Â 5 mm). The anatomic borders were the uncal recess -caudally-, the optical chiasm -rostally-and the white matter -superiorly-. Within these volumes of interest we extracted BOLD responses from all activated voxels (Schacher et al., 2006). Subsequently, functional scans were co-registered to the structural image. To avoid further spatial distortion, the activity analyses were carried out in native space instead of normalization to a standard space.
A general linear model (GLM) was built for each participant with regressors for the two conditions and the time courses of the regressors were convolved with the hemodynamic response function. For each subject, an F-test was carried out to determine regions more active during dynamic faces than during perception of dynamic landscapes. Results were considered significant at a threshold p < 0.05. This liberal threshold is chosen in clinical settings to avoid type two error that in the case of presurgical planning would underestimate BOLD activity.
Clinical evaluation of the amygdala's functional integrity largely depends on the patient's general level of activity and asymmetry between the two amygdalae. It was shown repeatedly that the applied paradigm leads to bilateral activity of the amygdala (Schacher et al., 2006;Toller et al., 2015), thus asymmetry between left and right amygdala indicates a loss of functional integrity within the afflicted amygdala. In order to best account for the individual functional integrity of the amygdala, an ordinal judgement of amygdala activity was conducted. Two members of the study team (HJ and BS) rated the amount of amygdalar activity as either 0 ¼ no amygdala activation, þ ¼ weak amygdala activity and þþ ¼ strong activity in amygdala. Ratings were carried out for left and right amygdala independently and raters were blind to the side of lesion as well as to the results in the recording of the depth electrodes.
The fMRI data included in the present work are part of a previously published dataset, including a larger cohort of epileptic patients and healthy controls, and performing a comprehensive contrast between BOLD responses to neutral vs. aversive faces (Schacher et al., 2006). In our present analyses we included a subset of these patients, for whom iEEG recording of the same task was available. Our goal was to examine the co-occurrence of BOLD and electrophysiological metrics, which is why for this study we kept the same BOLD analyses and responses as the previous study (Schacher et al., 2006). Here, we focused on the time-course of the BOLD response.

Intracranial data acquisition and preprocessing
The intracranial data was recorded starting from the day after electrode implantation. The recording was performed in the intensive monitoring unit under continuous video surveillance. Data from macroelectrodes was acquired at 4000 Hz with an ATLAS recording system (0.5-1000 Hz pass-band, Neuralynx, ww.neuralynx.com), imported to Fieldtrip  and down-sampled to 2000 Hz for further analysis. The iEEG was recorded against a common reference and then transformed to a bipolar montage for further analysis. To compare with a previous publication (Zheng et al., 2017), we also referenced to a contact in white matter but found no qualitative difference in results. Data from microelectrodes was acquired with 32000 Hz. We visually excluded time intervals showing epileptiform activity.

Spectral power analysis
Power envelopes were computed by filtering the iEEG signal with an IIR order 2 Butterworth filter for the following frequency bands: delta (1-4 Hz), theta (4-8 Hz), alpha (8-15 Hz), beta (15-30 Hz), gamma (30-60 Hz), high gamma (HGP, >60 Hz). We then applied the Hilbert transform to extract the analytic amplitude of the power envelope for each frequency band. The time-course of power was then smoothed with 1 s window in the time domain and baseline corrected using 1 s window before the beginning of each sequence to eliminate slow drifts. To account for inter-patient variability, we z-scored all data within each trial.
The baseline was a 2 s landscape clip that was presented before each block. To control whether the choice of the baseline could affect the outcome of the analysis, different choices of baseline were tested, i.e. taking as baseline the mean of all epochs, or the mean of the previous epoch. Results remained quantitatively similar across these different baseline choices.

Spike sorting
To characterize the axonal spikes of putative neuronal units, we used the program package Combinato (Niediek et al., 2016) (https://gith ub.com/jniediek/combinato). The following steps were performed in a fully automated way: 1) the signal was band-pass filtered (passband 300 Hz-1000 Hz, elliptic filter); 2) an amplitude threshold was set 5/0.6745 of the median of the rectified filtered signal; 3) events crossing the threshold were marked as putative neuronal spikes; 4) for each marked event, a segment of 2 ms duration was extracted from the signal after band-pass filtering (passband 300 Hz-3000 Hz, elliptic filter); 5) wavelet coefficients were obtained from the extracted segments and passed to superparamagnetic clustering; 6) well-defined data clusters from different "temperature" parameters of the superparamagnetic clustering were identified; 7) large clusters were subjected to iterative clustering, noise clusters were removed, and highly similar clusters were merged. After spike sorting, putative units were visually inspected and cluster with artifactual shapes rejected. Single units were visually identified by distinct peaks whose standard error was separable from the standard error of the tails on the y-axis (amplitude). Examples are shown in Fig. S2.

Correlation analysis
We quantified similarities between BOLD responses and delta/high gamma iEEG power at the single trial level, by correlating time-courses of activity, using Spearman's r. In order to capture the temporal evolution of both hemodynamic and electrophysiological responses at a similar time scale, the time-course of the iEEG power envelope was temporally smoothed and down-sampled to the BOLD-time points, i.e. at 1.5 s, and then smoothed over 5 time-points, i.e. 7.5 s, since the BOLD signal was recorded at 1.5 s time frames. This smoothing might penalize our metrics, as it sacrifices some of the high temporal resolution of iEEG and neuronal firing metrics but it has the advantage of taking into account the time-lag of BOLD responses with respect to iEEG, which peaks typically around 4-5 s after a stimulus is presented (Silva et al., 2005).

Statistical analysis
We statistically contrasted single-trial oscillatory power in response to AVERSIVE vs. NEUTRAL conditions in different frequency bands. To account for differences in variance across the group of patients while taking advantage of single-trial measurements, we used a linear mixed effects models which takes into account inter-individual variability.
Using this model, we tested time-point by time-point differences in DELTA and GAMMA power between NEUTRAL and AVERSIVE trials with the formula: Power~condition þ (1jPatients) The formula includes a fixed factor of condition (NEUTRAL/AVERSIVE) and a random factor of patients. For each patient, we included all aversive and neutral single trials, and we took into account inter-patient differences through the random factor of 'Patients'.
We corrected for multiple comparisons across time using a clusterbased permutation approach, where we permuted the labels of neutral and aversive trials across all participants, with 500 permutations .
To compare across different modalities, we extracted the difference in delta and high gamma power in response to aversive vs. neutral stimuli over the temporal clusters where the two conditions differ significantly, and termed this difference modulation index MI ¼ (AVERSIVE-NEUTRAL). Correlation coefficients were computed according to Spearman's definition. Statistical significance was established for p < 0.05. Reported p-values have been corrected for multiple comparisons using the false discovery rate (FDR).
We also quantified the specificity, sensitivity and accuracy of the high gamma MI with respect to seizure onset zone (SOZ). The effect was either HGP enhancement in the first 2 s or sustained delta decrease (see Results). We defined as true positive (TP) the occurrence of this effect in healthy amygdalae. We defined as false positive (FP) the occurrence of the effect in the amygdala within SOZ. We defined as false negative (FN) the absence of the effect in nSOZ amygdalae. We defined as true negative (TN) the absence of the effect in the amygdala within SOZ.
Finally, we computed correlations among neuronal firing and across modalities (firing rates, delta power, HGP and BOLD) using Spearman's r. In order to extract features of neural firing across participants we defined the Firing Rate Variability (FRV) and the Spatial Synchronization Index (SSI). The FRV is the difference in average firing rate between AVERSIVE and NEUTRAL conditions, and captures the effect of the modulation induced by the task on the firing rate of the single neuron. Essentially, it is a similar metric as the MI, but it is computed over the whole timeinterval of neuronal firing, while the MI is computed over significant clusters of slow and fast oscillatory power. The SSI is defined as the average pairwise correlation of each neuron with neighboring neurons recorded in the same amygdala but on different microwires. This correlation is therefore a direct estimation of the local synchronization among neurons.
Multiple linear regression was implemented to quantify the association between neuronal firing and macroscale measurement (HGP, BOLD response) as dependent variable, with FRV and SSI as independent variables.

Conclusions
In summary, we have shown, for the first time, common correlates across different physiological signals of the human amygdala, spanning from single-neuron, to local field potentials and hemodynamic responses, showing that delta and high gamma power are locally correlated with single neuron activity, and more globally linked to hemodynamic responses.