Neural computations underlying contextual processing in humans.

Context shapes our perception of facial expressions during everyday social interactions. We interpret a person's face in a hostile situation negatively and judge the same face under pleasant circumstances positively. Critical to our adaptive fitness, context provides situation-specific framing to resolve ambiguity and guide our interpersonal behavior. This context-specific modulation of facial expression is thought to engage the amygdala, hippocampus, and orbitofrontal cortex; however, the underlying neural computations remain unknown. Here we use human intracranial electroencephalograms (EEGs) directly recorded from these regions and report bidirectional theta-gamma interactions within the amygdala-hippocampal network, facilitating contextual processing. Contextual information is subsequently represented in the orbitofrontal cortex, where a theta phase shift binds context and face associations within theta cycles, endowing faces with contextual meanings at behavioral timescales. Our results identify theta phase shifts as mediating associations between context and face processing, supporting flexible social behavior.

INTRODUCTION

Memories of our daily experiences are closely linked to the context of the event. Context provides situation-specific framing to resolve ambiguity and generate predictions that guide our interpersonal behavior. A critical function of contextual framing is to shape our perception of faces as we interpret ongoing social interactions and determine appropriate actions. For example, we interpret the face of a person with a gun as “threatening” but perceive the face of the same person with a birthday cake as “happy.” These inferences allow us to resolve uncertainties by integrating background information (e.g., the gun versus the cake) with sensory perceptions of the face to construct a coherent representation. Interpreting situation-specific facial signals is critical for human communication because we use contextual information to predict our environment, take actions based on anticipated outcomes, and adaptively modify our behavior.

Despite the importance of contextual processing in human social communication, how the brain encodes context to modulate the perception of facial emotion is a long-standing question (Hammerschmidt et al., 2018; Weiss et al., 2020). Convergent evidence in animals and humans implicates the hippocampus (HPC) for processing contextual (Davachi and DuBrow, 2015; Maren et al., 2013) and mnemonic information (Lisman et al., 2017), whereas the amygdala (AMY) codes for events’ emotional significance (Morrison and Salzman, 2010; Phelps and LeDoux, 2005; Zheng et al., 2017). Human AMY neurons have a privileged role in social cognition, with cells selectively responding to faces (Mormann et al., 2015), and AMY lesions impair recognition of facial emotions (Adolphs et al., 1994). Thus, the AMY and HPC are hypothesized to bind contextual information with faces to produce social meanings (Cao et al., 2021; Vuilleumier et al., 2004), but the circuit interactions underlying this cognitive ability have yet to be elucidated in humans (Hortensius et al., 2016; Wieser and Brosch, 2012; Yang and Wang, 2017). Context also provides contingencies to predict upcoming stimuli, but how humans accomplish this cognitive function is also debated. Recent theoretical and experimental evidence suggests that the orbitofrontal cortex (OFC) plays an essential role in tracking inferred states and transmitting information about expectations, rewards, and outcomes (O’Doherty et al., 2001; Saez et al., 2018; Stalnaker et al., 2014; Wikenheiser and Schoenbaum, 2016; Zangemeister et al., 2016). These findings imply that the OFC is well suited to abstract contextual information from medial temporal lobe regions such as the AMY and HPC (Sep et al., 2020; Zhang et al., 2018). The OFC and HPC construct cognitive maps with different but overlapping encoding properties (Wikenheiser and Schoenbaum, 2016), and the OFC and AMY encode decision-making variables, such as choice value and feedback (Morrison and Salzman, 2010; O’Doherty et al., 2001; Zangemeister et al., 2016). Such a functional overlap suggests the existence of an AMY-HPC-OFC tripartite network supporting contextual integration during emotional processing. Simultaneous recordings from these pathways are needed to define how context-face association signals interact across these frontal and temporal brain regions.

To address these questions, we employed human intracranial stereo-electroencephalography (SEEG) recordings obtained simultaneously from the AMY, HPC, and OFC, coupled with an experimental task that sequenced the encoding, maintenance, and integration of context-to-face associations. We show behaviorally that subjects implicitly link the preceding image’s valence with the subsequent sensory information of the face. Indeed, the context’s valence robustly modulates subjects’ perception of facial expression because a face paired with an aversive image is interpreted as negative, and the same face following a pleasant image is considered positive.

At the circuit level, we show that theta oscillations flexibly mediate frontotemporal communications, tracking encoding and maintenance of context within the AMY-HPC and OFC-HPC network, respectively, followed by context-face associations within the OFC-AMY network. Cross-regional theta-gamma phase-amplitude coupling (PAC) provides a means of directional interactions within this tripartite network that facilitates temporal contextual modulation. We show that the high gamma activity (HGA)-theta phase shift promotes integration of the emotional temporal context into facial expression. The dynamics and behavioral relevance of this phenomena are reminiscent of theta phase precession, a temporal coding mechanism for associative learning (O’Keefe and Recce, 1993; Salman et al., 2020; Terada et al., 2017). Our data support a model where theta-gamma interactions flexibly shape functional connections within the AMY-HPC-OFC network and the sequential binding of temporal context and face within theta cycles. By conferring faces with contextual valence, our report reveals a putative neural mechanism in humans that infers facial expressions at behavioral timescales.

RESULTS

Emotional context biases subsequent facial perception

We recorded intracranial SEEGs from the AMY (20 electrodes), HPC (18 electrodes), and OFC (19 electrodes) in 8 pre-surgical individuals with epilepsy (5 males, 3 females; see Table S1 for subject information and Figure 1D for electrode locations) while they performed a decision-making task that required integration of contextual information and facial expression. Each trial consisted of five task stages (Figure 1A): a fixation baseline (fixation; 0.5 s), display of a context image (context; 1 s), a maintenance period with a blank screen (maintenance; 0.5 s), presentation of a neutral face (face; 1 s), and the valence rating of the face (decision; self-paced). At the decision stage, subjects rated the emotional valence of each face by selecting the corresponding color (warm color = positive valence, cold color = negative valence). Each neutral face (N = 35 faces; 17 males and 18 females) was paired with three context images (one positive, one neutral, and one negative). This three-to-one mapping between context and face provided the ability to test the flexible accommodation of different contextual valences within the same face. We then investigated whether the context image’s valence parametrically modulated the valence rating of the face. We found that, consistent across all subjects (Figures S1A and S1B), subjects’ valence ratings of faces were positively correlated with the valence scores of the preceding context images, which had been independently derived from the International Affective Picture System (IAPS; Coan and Allen, 2007; Figure 1B; r = 0.511, p = 2.236 × 10−11, Spearman correlation), suggesting that subjects implicitly integrated the contextual valence with neutral faces. Neutral faces could also contain valence information because subjects tended to rate specific neutral faces more positively or negatively (Lin et al., 2021; Oosterhof and Todorov, 2008; Sutherland et al., 2013; Figure S1E), which also confirmed that subjects did not ignore the faces or rate only based on context valences. To confirm that such contextual modulation is not solely driven by the valence information embedded in the neutral faces, we subtracted the valence rating of each neutral context-face pair from other positive or negative context-face pairs. The correlation remained significant after adjusting for the intrinsic face valence, and the adjusted correlation was stronger for negative contexts than positive contexts (Figures S1C and S1D; positive contexts: r = 0.31, p = 7 × 10−3; negative contexts: r = 0.46, p = 4 × 10−8).

Figure 1.

Contextual modulation task design, behavioral performance, and electrode locations

(A) Five-stage contextual modulation task: baseline (fixation, 0.5 s), context image (context, 1 s), context maintenance (maintenance, 0.5 s), face image (face, 1 s) and valence rating of the face (decision, self-paced). Each face was paired with a positive, neutral, or negative valence context. After the face display, subjects rated the valence of each face by moving the blue triangle to the appropriate position on the color bar (warmer colors, more positive valence; cooler colors, more negative valence).

(B) Subjects’ valence ratings of the faces were positively correlated with the valence of the context images. The valence of the context images and faces were normalized within each subject (STAR Methods), with larger/smaller numbers indicating more positive/negative valence.

(C) Distributions of contextual modulation strengths expressed as a ratio of face valence rating over context valence (STAR Methods). The red area denotes negatively modulated trials (e.g., a face rated as negative when preceded by a negative context image), and the blue area contains positively modulated trials (e.g., a face rated as positive when preceded by a positive context image). A gray area indicates the opposite modulation condition (e.g., a face rated as negative despite a prior positive context image). The pink line at zero represents no modulated trials (e.g., a face rated as neutral regardless of context image valence).

(D) The locations of electrodes from all 8 subjects were rendered onto a Colin27 template brain (Holmes et al., 1998) based on a high-resolution atlas template (red, AMY; blue, HPC; green, OFC). A, anterior; P, posterior; D, dorsal (D); V, ventral; L, left; R, right.

To quantify the magnitude of contextual modulation (Figure 1C), we computed a contextual modulation strength for each trial based on the subject’s valence rating of the face relative to the IAPS valence score of the corresponding contextual image (STAR Methods). We found that 87% of trials (822 trials across subjects) demonstrated a congruent contextual modulation direction. Specifically, subjects rated neutral faces paired with pleasant images more positively. In contrast, faces that followed aversive images were perceived as more negative. We focused our analyses on trials demonstrating congruent modulation, and the small numbers of incongruent and non-modulated trials were collapsed into a vector of non-interest and excluded from further analyses (Figure 1C; Nincongruent modulation = ~8 trials per subject; Nno-modulation = ~2 trials per subject). Trials with negative context valence showed a stronger modulation strength compared with trials with positive valence (F (1, 820) = 8.642, p = 0.003), but the range still allowed us to examine emotional valence across the positive and negative valence continuum. These behavioral findings demonstrate a robust carryover effect where the valence information from context images was sequentially integrated to faces and, thus, influenced subjects’ subsequent facial perception.

Theta oscillations promote frontotemporal interactions during contextual modulation

Building on the behavioral effects, we examined the functional connectivity across the AMY-HPC-OFC network during contextual modulation. Connections among brain regions are thought to be dynamically modulated by synchronization of low-frequency neural oscillations, reflecting rhythmic fluctuations of neuronal excitability (Wang, 2010). Theta oscillations, prominent in the AMY and HPC, play a critical role in promoting cross-regional communication during aversive memory retrieval and contextual fear processing (Seidenbecher et al., 2003; Zheng et al., 2017, 2019). We also found task-evoked theta power modulation across different task stages (Figure S2A).

We then measured theta phase synchrony among all possible electrode pairs between the AMY, HPC, and OFC using the weighted phase lag index (WPLI). The WPLI quantifies the cross-regional phase synchrony using the imaginary components of the cross-spectrum between signals from each electrode pair. This method is insensitive to spurious increases in synchrony induced by volume conduction and random noise (Vinck et al., 2011). For each subject, low-frequency (4–12 Hz) phase synchrony was computed, normalized to the baseline, and projected to a time-frequency plot (Figure 2A). We then measured the averaged phase synchrony across all subjects within the theta band (4–12 Hz) as a function of time (Figure 2B). Significant theta phase synchrony (p < 0.01, permutation test; STAR Methods) emerged within the AMY-HPC-OFC circuit at different task stages. During context presentation, the network was dominated by increased theta synchrony in the AMY-HPC network and then transitioned to an enhanced temporal-frontal (AMY-OFC and HPC-OFC) interaction during context maintenance and face display stages. Finally, during the decision stage, we observed a sustained rise in AMY-OFC theta synchrony. These findings highlight theta oscillations as prominent mediators of cross-regional interactions, selectively scaling up and down specific cross-regional functional connectivity during emotional sequence processing. During periods of significantly increased phase synchrony (4–6 Hz for AMY-HPC and 8–11 Hz for HPC-OFC; Figure 2B), the phase synchrony simultaneously decreased at the nearby frequencies (7–8 Hz for AMY-HPC and 4–5 Hz for HPC-OFC). To identify directional influences within this tripartite network, we computed Granger causality indices (Ding et al., 2000; Schelter et al., 2006) for each cross-regional electrode pair. Significant directional influence (p < 0.01, permutation test; STAR Methods; Figure 2C) across subjects was plotted individually at different task stages (fixation, context, maintenance, face, and decision), with arrows indicating the directional influence between connected regions and the line thickness representing the ratio of electrode pairs showing significant Granger causality. We observed that the tripartite network started with bidirectional interactions between the AMY and the HPC during context presentation and shifted to a unidirectional influence from the HPC to the OFC at the maintenance stage. When the face was presented, bidirectional communications were established between the AMY and OFC and were later dominated by a unidirectional influence from the OFC to the AMY and HPC at the decision stage. Consistent directional results were also observed when using the cross-frequency directionality analysis (Figures S2C–S2F).

Figure 2.

Theta oscillations promote fronto-temporal interactions during contextual modulation

(A) Low-frequency (4–12 Hz) phase synchrony quantified as WPLI and averaged across trials for cross-regional electrode pairs (top, AMY-HPC; center, HPC-OFC; bottom, AMY-OFC) of a representative subject. Warm/cold colors denote increased/decreased phase synchrony relative to the baseline, respectively. Significant power increases/decreases are highlighted with white contours. Vertical dashed lines mark task stages.

(B) Theta phase synchrony averaged across all subjects, plotted as a function of time for each cross-regional pair (purple, AMY-HPC; blue, HPC-OFC; green, AMY-OFC; shaded color, SEM). Color bars at the top indicate time windows with significant theta phase synchrony.

(C) Directional maps (i.e., Granger causality) between cross-regional electrode pairs for all subjects during five task stages. Colored dots represent different brain regions (red, AMY; blue, HPC; green, OFC). Regions with at least one electrode pair showing significant Granger Causality (p < 0.05, permutation test) are connected by lines, with arrows indicating the directional influence between connected regions. The line thickness increases with the ratio of electrode pairs with significant Granger causality.

Cross-regional phase amplitude coupling predicts perceptual bias of faces

Having shown dynamic modulations of effective connectivity in the AMY-HPC-OFC network during contextual processing, we wanted to determine whether enhanced synchronous neural activities would predict behavior outcomes. Guided by the functional role of theta-gamma phase amplitude coupling (PAC) in supporting context-item associations (Tort et al., 2009) and linking items in sequences (Axmacher et al., 2010; Heusser et al., 2016; Lisman and Jensen, 2013), we reasoned that this neural mechanism may serve to bind context and face associations. To test this, we first measured cross-regional PAC in the tripartite network (as shown in Figures S3A and S3B, theta-gamma PAC is most prominent across the full frequency band) and then examined the trial-by-trial relationship between PAC and subjects’ contextual modulation strength. For each electrode pair, we first calculated the cross-regional PAC for both directions (Atheta -> BHGA: theta phase from A modulating HGA [70–250 Hz] from B; Btheta -> AHGA: vice versa) and identified electrode pairs with significant PAC (permutation test; STAR Methods; AMY-HPC, n = 31 pairs; HPC-OFC, n = 28 pairs; AMY-OFC, n = 27 pairs). We then plotted the PAC coupling strength across all of these PAC significant electrode pairs as a function of time. The strength of PAC was quantified as the Euclidean norm of the complex signal (Samiee and Baillet, 2017) (i.e., the low-frequency phase from A/B and the high gamma amplitude from B/A; STAR Methods) and Z scored relative to the fixation baseline. We observed enhanced PAC strength emerging across different task stages within the AMY-HPC-OFC network (Figure 3A). We then quantified the correlation between theta-gamma PAC and contextual modulation strengths at the single-trial level and plotted correlation values as a function of time (Figure 3B). Two sequential task stages, maintenance and face presentation, were characterized by significant positive correlations between theta-gamma PAC and contextual modulation strength (Figures 3C and 3D). Contextual modulation strength was related to the cross-regional PAC strength between HPC theta and OFC HGA during the maintenance stage (Figure 3C; rHPC -> OFC = 0.494, pHPC -> OFC = 2.19 × 10−11, Spearman correlation) and to PAC strength between the OFC theta and AMY HGA during the face presentation stage (Figure 3D; rOFC -> AMY = 0.393, pOFC -> AMY = 5.11 × 10−7, Spearman correlation), which is also valid at the single-subject level (Figures S3F and S3G). Intra-regional PAC (e.g., theta phase from the HPC modulating HGA from the HPC) was not significantly correlated with contextual modulation strength (Figures 3E and 3F; Figures S3C–S3E). These results suggest that the HPC -> OFC network is modulated during maintenance of contextual information for subsequent association with facial information via the OFC -> AMY pathway.

Figure 3.

Cross-regional phase amplitude coupling (PAC) correlates with contextual modulation strength

(A) Cross-regional PAC averaged across subjects, and Z score normalized to baseline and plotted as a function of time. Theta-gamma PAC for cross-regional pairs (purple, AMY-HPC; blue, HPC-OFC; green, OFC-AMY) are plotted as solid and dashed lines (shaded area, SEM) for opposing modulation directions (e.g., theta oscillations from the AMY modulating HGA from the HPC and vice versa). Colored bars at the top indicate time periods with significant cross-regional PAC.

(B) Correlations between cross-regional theta-gamma PAC and contextual modulation strength averaged across subjects and plotted as a function of time. The plotting format is similar to (A).

(C) Averaged HPC-> OFC PAC (Z-scored over baseline) at context maintenance (y axis) and the contextual modulation strength (x axis). Each dot refers to one trial, with PAC strength averaged across all pairs within the same region from the same subject. Correlation plots for individual subjects are shown in Figure S3F.

(D) Averaged OFC-> AMY PAC (Z-scored over baseline) during face display (y axis) and contextual modulation strength (x axis) (similar format as in C). Correlation plots for individual subjects are shown in Figure S3G.

(E and F) Correlation coefficients between theta-gamma PAC and contextual modulation strength for intra-regional (diagonal values) and cross-regional (the remaining values) interactions during the maintenance (E) and face stages (F), with significant correlations marked by asterisks. Significance level for (A), (B), (E), and (F): p < 0.05, non-parametric permutation test.

Theta-gamma phase shift supports temporal binding of emotional sequences

Given that the two consecutive stages of the task, context maintenance and face presentation, selectively showed behaviorally relevant theta-gamma coupling, we next investigated how theta-gamma coupling supports the temporal coding of context and face associations. Classic work has identified theta phase precession, a phenomenon where neurons systematically shift their phase of firing relative to theta rhythm while an animal traverses the place field, as a critical mechanism for compression of temporal sequences (O’Keefe and Recce, 1993; Skaggs et al., 1996; van der Meer and Redish, 2011). At the neural population level, transient gamma bursts have been proposed to represent multiple items in a sequence that are temporally organized by ascending and descending phases of theta cycles (Lisman and Jensen, 2013). Therefore, we hypothesize that the timing of gamma bursts relative to theta phase shift across event transitions, which we here refer to as “theta-gamma phase shift,” is a potential mechanism for integrating contextual sequences. Motivated by methods used to detect gamma bursts (Lundqvist et al., 2018) in nonhuman primates (NHPs), we first extracted high-gamma events at each frequency bin (70–250 Hz) that exceeded 2 standard deviations of the mean high-frequency power within each trial, lasting for a minimum of 3 cycles (Figures S4A and S4B). Using these criteria, we detected bursting dynamics in the high-gamma band (70–250 Hz) with a narrow frequency bandwidth (22.486 ± 4.232 Hz, mean ± SEM) and short duration (0.123 ± 0.037 s, mean ± SEM) (Figures S4C and S4D), in line with characteristics found in NHPs (Lundqvist et al., 2018). We then assessed for the presence of theta-gamma phase shifts by examining whether gamma bursts at later time points during the task occur at earlier theta phases (see the single trial examples in Figures S4E–S4H and group-level assessments in Figures S5A–S5C), resulting in a negative correlation between gamma bursts and theta phases and time (O’Keefe and Recce, 1993; Salman et al., 2020; Skaggs et al., 1996; van der Meer and Redish, 2011). First, we separately calculated the occurring phase of HGA bursts in OFC relative to HPC theta oscillations (Figures 4A and 4B) and the HGA burst in the AMY relative to the OFC theta oscillations (Figures 4E and 4F; STAR Methods). Second, we quantified the circular-linear correlation between the occurring phase of HGA bursts relative to theta oscillations and time (solid gray lines in Figures 4A, 4B, 4E, and 4F) and tested the significance of the pattern (p < 0.05, permutation test; STAR Methods). We observed a robust theta-gamma phase shift during the context maintenance and face presentation stages. The observed phase shifts were not driven by transient frequency changes in HGA or theta oscillations (Figures S5D–S5I), changes in theta/gamma power (Table S2), or bimodal phase locking across the maintenance and face stages (Figures S5J–S5L; Table S3). The OFC HGA bursts were nested around late HPC theta phases (180°–360°) during context maintenance and shifted to early phases (0°–180°) while the face was presented (pmaintenance = 0.006, pface = 0.003, Rayleigh test; Figures 4C and 4D). The AMY HGA bursts also showed phase clustering that shifted from late phases during context maintenance to coalesce at early OFC theta phases (0°–180°) after face presentation (pmaintenace = 0.032, pface = 0.008, Rayleigh test; Figures 4G and 4H). The theta-gamma phase shift across the event boundary between the maintenance and the face presentation stages was consistent across subjects, as evidenced by the strongest and weakest representative single subject examples (Figures 4C, 4D, 4G, and 4H) and the group phase distribution plots (Figures 4I–4L). In total, we found significant theta phase shifts in 57.1% (24 of 42) of HPC theta and OFC gamma electrode pairs and 57.9% (22 of 38) OFC theta and AMY gamma electrode pairs. Similar results were observed when using different phase estimation approaches (Figures S4E–S4H and Table S3). The characteristic features of observed phase shifts (Figures S4I–S4K) are comparable with rodent studies of phase precession (Feng et al., 2015; Schmidt et al., 2009; Terada et al., 2017).

Figure 4.

Theta-gamma phase shift

(A and B) Theta phase plots of OFC HGA bursts across all trials in two representative electrode pairs with significant theta-gamma phase shift (A, strongest; B, weakest) in the maintenance and face stages. Each dot represents an OFC HGA burst relative to HPC theta phases. Gray lines mark the circular-linear correlation between the occurring phases of HGA bursts relative to theta oscillations and the time across the maintenance and face stages (results computed for each task stage are shown in Figures S5J–S5L). Numbers in the parentheses indicate the circular-linear correlation coefficient, and the asterisks indicate significance level (*p < 0.05, **p < 0.01, ***p < 0.001, permutation test).

(C and D) Distributions of OFC HGA bursts relative to HPC theta phases as plotted in (A) and (B) for the maintenance (dark blue) and face stages (light blue). Vertical lines mark the mean occurring phase.

(E and F) Theta phase plots of AMY HGA bursts across all trials in two representative electrode pairs with significant theta-gamma phase shift (E, strongest; F, weakest). Each dot represents an AMY HGA burst relative to the OFC theta phase (format similar to A and B).

(G and H) Distributions of AMY HGA bursts relative to OFC theta phases as plotted in (E) and (F) for the maintenance (dark green) and face stages (light green). Vertical lines marked the mean occurring phase.

(I and J) Distributions of OFC HGA burst occurrences across subjects relative to HPC theta phases for the maintenance (I) and face stages (J).

(K and L) Distributions of AMY HGA burst occurrences across subjects relative to OFC theta phases for the maintenance (K) and face stages (L).

Although theta-gamma phase shift was evident at event transitions, whether early and late theta phases encode distinct information from the event sequence remained unknown. To address this question, we performed Bayesian decoding analysis (Terada et al., 2017) to reveal different task information and its associated neural representations. First, we decoded the time information based on HGA bursts from the OFC and AMY (i.e., whether HGA bursts at time A encoded similar HGA bursts features at different time points). Figure 5A demonstrates the reconstructed time information at each 5-ms time step in a single example trial. Consistent with the theta-gamma phase shift results, context maintenance and face presentation were the only task stages that encoded future neural states. Specifically, OFC HGA bursts during context maintenance reflected neural representations in the context maintenance and face display stages, whereas AMY HGA bursts at face display mirrored neural representations in the face display and decision stages. The encoded current and future states cycled at approximately 8 Hz (~125 ms, 8 cycles/s), suggesting a theta oscillation-mediated temporal coding mechanism. We then decoded the dual-state cycling based on the timing of HGA bursts relative to theta oscillations during the context maintenance (Figure 5B) and face presentation stages (Figure 5C). In both task stages, the dual states were compressed in a theta sequence, with early theta phases (0°–180°) reflecting current neural states in which the HGA was sampled (dashed rectangles in Figures 5B and 5E) and late theta phases (180°–360°) representing the future neural states (solid rectangles in Figures 5B and 5E).

Figure 5.

Early and late theta phases encoded distinct events

(A) Probability density of time information computed using OFC (top) or AMY (bottom) HGA from a representative subject.

(B) Probability density of time information as a function of HPC theta phase, estimated at the context maintenance stage.

(C and D) Decoded probabilities of context valence (brown) and face rating (gray) based on late HPC theta phase (C, solid rectangle in B) and early HPC theta phases (D, dashed rectangle in B) for trials with strong and weak contextual modulation. Horizontal dashed lines denote the decoding chance level.

(E) Probability density of time information as a function of OFC theta phase, estimated at the face stage from the same trials as plotted in (B).

(F and G) Decoded probabilities of context valence (brown) and face rating (gray) based on late OFC theta phase (F, solid rectangle in E) and early OFC theta phases (G, dashed rectangle in E) for trials with strong and weak contextual modulation (similar format as C and D).

Error bars in (C), (D), (F), and (G) indicate SEM. *p < 0.05; **p < 0.01.

We next estimated the type of information (context valence versus face rating) represented in the current and future states. During context maintenance, we found that late HPC theta phases (180° – 360°) encoded context valence information for the future face presentation stage (Figure 5C; p = 0.007, permutation test; STAR Methods). During the subsequent face presentation stage, context valence information was encoded at early OFC theta phases (0°–180°) (Figure 5G; p = 0.0023, permutation test; STAR Methods), whereas late OFC theta phases (180° – 360°) reflected face rating (p = 0.008) and context valence information (p = 0.031) for the future decision stage. This phase-specific encoding mechanism appeared only in the strong contextual modulated trials (split by contextual modulation strength, the top half of the trials; left panels in Figures 5C, 5D, 5F, and 5G) and was absent for the weak contextually modulated condition (the bottom half of the trials; right panels in Figures 5C, 5D, 5F, and 5G). These findings suggest that contextual modulation of facial expression is mediated by frontotemporal theta oscillations, with context valence information carried across event transitions via theta phase advancement encoded from HPC late theta phase to OFC early theta phase and compressed with face rating information in theta cycles.

DISCUSSION

It is widely acknowledged that contextual information modulates how we perceive facial expressions. Although prior work has demonstrated behavioral effects (Barrett and Kensinger, 2010; Mobbs et al., 2006; Wieser and Brosch, 2012) and their associated brain regions (Lee and Siegle, 2014; Wieser and Brosch, 2012), less is known about how these regions coordinate to induce affective biases. Here, to examine the neural signatures underlying contextual modulation of facial perception, we recorded direct neural signals from the AMY-HPC-OFC circuit, a well-established contextual processing network (Maren et al., 2013), while subjects performed valence ratings of neutral faces preceded by emotional contexts. Contextual modulation in the present task refers to temporal context based on the emotional stimulus (contextual modulator) preceding face presentation and affecting the face valence rating. Our results reveal a dynamic functional connectivity map within the AMY-HPC-OPC network, with theta oscillations (4–12 Hz) selectively tuning up and down specific cross-regional interactions (i.e., phase synchrony) at different task stages (Figure 2). Such cross-regional influence via theta oscillations modulates local neuronal activity (i.e., high-frequency activity), with the strength of cross-frequency modulation at maintenance and during face presentation reflecting subjects’ perceptual bias of the face valence (Figure 3). The degree of this perceptual bias relies on successful encoding of the context and face information by the timing of HGA (Figures 4 and 5), which is precisely coordinated by cross-regional theta-gamma phase shifts. These findings provide evidence of a network-level mechanism in the human brain to encode contextual information for interpreting facial expressions, a function that is critical for our daily social interactions.

In this study, we find that theta synchrony dynamically shapes interactions in the AMY-HPC-OFC network to support contextual processing (Figure 2), modulating anatomical connections to adapt to different task demands (Akam and Kullmann, 2014). We observed that the OFC and medial temporal lobe (MTL) networks dynamically engage (i.e., enhanced theta synchrony) and disengage (i.e., lack of increased theta synchrony) throughout the task to optimize context-dependent decision-making. The periods of significantly increased phase synchrony between the HPC and AMY (4–6 Hz) or OFC (8–11 Hz) were accompanied by nearly simultaneous decreases in phase synchrony for the same region pairs (7–8 Hz for AMY-HPC and 4–5 Hz for HPC-OFC). This dynamic is consistent with the oscillatory multiplexing framework, where the same structure can use different frequencies for communicating with different other structures, and the phase synchrony decrease can contribute to segregating these distinct communication channels (Akam and Kullmann, 2014; Watrous et al., 2013; Zheng et al., 2019).

During processing of context images, AMY-HPC theta synchrony dominates, whereas OFC interactions with the AMY-HPC network is downmodulated. Selective enhancement of theta coherence could be crucial for the AMY-HPC network to compute behaviorally relevant information, such as the context’s valence (Redondo et al., 2014; Zheng et al., 2019). When keeping the information online is essential for face-context integration and decision-making, the theta-mediated network switches to an OFC- and MTL-centric pathway between the OFC-HPC and OFC-AMY network (Hoover and Vertes, 2007; Jay et al., 1989; Yizhar and Klavir, 2018). Previous studies of rodents have linked theta-associated segregation and integration of frontal-MTL connectivity with contextual and fear processing. In context-dependent working memory tasks, HPC-to-prefrontal theta-mediated connectivity is prominent during memory maintenance or when cued by a context (Place et al., 2016; Siapas et al., 2005), whereas retrieving the maintained memory or planning the choice response engages prefrontal-to-HPC influence (Hallock et al., 2016; Place et al., 2016). Bidirectional OFC-AMY interactions support valence assessment (Likhtik and Paz, 2015) during fear learning, whereas the processed valence information is thought to transmit from the OFC to the AMY for modifying behavioral responses (Yizhar and Klavir, 2018). Our results specify theta synchrony as a selective mechanism that adjusts online and offline connectivity between the OFC and MTL, depending on task demands, boosting MTL-OFC interactions during context processing and maintenance and augmenting the OFC-MTL network during context-face integration.

How do bidirectional interactions between the OFC and MTL promote sequential binding of context and face? An influential theory posits that theta oscillations facilitate temporal coding and integration of sequential information via spike-timing-dependent plasticity (Lisman and Jensen, 2013). The dynamics of theta-gamma phase shifts in the present study is reminiscent of theta phase precession, a phenomenon where hippocampal place cells systematically shift to earlier phases of firing relative to the theta rhythm as the animal traverses the place field. Theta phase precession is thought to be a critical mechanism for compression of temporal sequences (O’Keefe and Recce, 1993; Skaggs et al., 1996; van der Meer and Redish, 2011). In rodents, theta phase precession is also observed between the HPC and medial prefrontal cortex (Jones and Wilson, 2005; van der Meer and Redish, 2011) and, more recently, has been extended to non-spatial domains such as binding of odor and sound (Terada et al., 2017). We examined whether a similar neural mechanism can be distilled at a population level with local field potentials rather than single neurons. Guided by theoretical models (Lundqvist et al., 2011) that predict that gamma bursts represent cell ensembles in discrete attractor states and by empirical works (Lopes-Dos-Santos et al., 2018; Lundqvist et al., 2016) that show that gamma bursts correspond to neuronal spiking conveying informative sensory information, we tested the hypothesis that the gamma burst timing will shift to earlier phases of theta oscillation between the OFC and MTL during context-face integration. Our findings of theta-gamma phase shift between the HPC and OFC as well as the OFC and AMY are consistent with the rodent literature and extend this putative mechanism for integrating affective representations into humans. Theta oscillations in bats (Ulanovsky and Moss, 2007), monkeys (Stewart and Fox, 1991), and humans (Bohbot et al., 2017; Jacobs et al., 2007; Zheng et al., 2017) are transient and non-sinusoidal (Figures S4E–S4H) compared with the continuous theta oscillations reported in rodents. Therefore, consistent with the asymmetric depolarization model proposed by Mehta et al. (2002) and experimental observations in bats (Eliav et al., 2018) and humans (Qasim et al., 2021), phase precession might serve as a universal neural code across different species and different cognitive dimensions (spatial and non-spatial tasks) regardless of the presence of sinusoidal oscillations.

A critical phase precession prediction is that the current state is anchored at early phases of the theta oscillation, whereas the future state is encoded at late phases (Dragoi and Buzsaki, 2006; Foster and Wilson, 2007; Skaggs et al., 1996). The transition from current to future hypothetical states should occur rapidly at behavioral timescales. Using Bayesian decoding, we demonstrate that HGA bursts in the OFC and AMY encoding current and future states are swiftly alternating at ~125-ms cycles. This repetitive sampling of current and future states at ~8 Hz is reminiscent of sub-second cycling of possible future states in rodents during navigation (Kay et al., 2020) and suggests that a similar dynamic neural representation occurs in humans for coding of prospective information. During the maintenance and face presentation stages of the sequential task, temporal coding of contextual information shifts from late theta phases in the HPC (180°–360°) to early theta phases in the OFC (0°–180°). These results suggest that theta dynamics provide a temporal organization to “hand off” contextual information from the HPC to the OFC in sub-second cycles, possibly facilitating “offloading” of encoding spaces for upcoming events while maintaining the context information for further decision (Eichenbaum, 2017). Theta phase shift between the HPC and in its anatomically connected OFC in humans implies a neural representation capable of facilitating the precise timing of long-range information transfer (Harris and Gordon, 2015). Theta precession, classically implicated in organizing spatial information, is also capable of sequencing the current and future representations of non-spatial events, such as context and face.

At the population level, an essential requirement for the theta precession-like dynamics would be a sufficient contribution of precessing individual neurons to a net population activity. Such dynamics might not be observable in the spatial context because of sparse spatial encoding of place cells (Barnes et al., 1990; Jung et al., 1994), but in non-spatial tasks, the proportion of simultaneously precessing neurons is larger (Terada et al., 2017). Also, a population-level phase shift is compatible with the theta-gamma PAC when the phase-shift is accompanied by phase-dependent activity modulation. Theta-gamma phase shift in the present study shows multiple similarities with theta phase precession (Maurer and McNaughton, 2007; O’Keefe and Recce, 1993), especially in the non-spatial context (Pastalkova et al., 2008; Qasim et al., 2021; Robinson et al., 2017; Takahashi et al., 2014; Terada et al., 2017). This includes the negative circular-linear correlation between the gamma activity and theta phase across the task stages (Figure 4) as well as the decodability of current and future states based on the HGA during different theta phases (Figure 5). Simultaneous single-neuron recordings from the AMY, HPC, and OFC during contextual modulation behavioral paradigms could inform the extent of similarity between these phenomena. As pointed out by Jensen (2001), the computational advantage of the systematic change in theta-gamma phase is an increased information transfer capacity. Specifically, if the temporal context of the information encoded by spiking activity or HGA (a proxy measure of spiking) (Ray and Maunsell, 2011) is determined by the underlying theta phase, then the temporal context could be decoded by the downstream reader.

We reveal neural computations within the AMY-HPC-OFC network that support context-dependent modulation of facial perception. Specifically, contextual information is processed within the MTL, projects to the OFC via an HPC-to-OFC route, and integrates with facial sensory information through AMY-OFC bidirectional interactions. We identify a critical role of theta phase shifts in mediating long-range communications between the MTL and OFC to incorporate past and future information. By unveiling the functional anatomy and temporal dynamics of these neural representations, we show how discrete emotional events are sequentially coded to shape our perception and guide our social behavior.

Limitations of the study

First, we demonstrate context-specific modulation of facial perception, with emotional valence from context biasing participants’ interpretation of the subsequent facial expression. Because time constraints in the clinical setting, we did not collect participants’ valence ratings of the context images. Instead, when computing the contextual modulation strength, we used valence ratings from the IAPS rated by an independent group of subjects. Although previous literature has reported high reliability of the IAPS in different age groups (Grühn and Scheibe, 2008) and diverse cultural backgrounds (Lohani et al., 2013; Soares et al., 2015), individual participants’ intrinsic (unmeasured) valence ratings for a context image might differ from crowd-sourced ratings. Second, the contextual processing we studied here is a temporal contextual modulation, with context and face presented sequentially because this task structure disentangles neural dynamics at each stage. It will be interesting to extend the study under a real-life setup, with context and face presented simultaneously. Third, we also demonstrate the existence of phase shifts with high gamma bursts at later time points during the trial occurring at earlier theta phases. The phenomenon of phase shifts observed in this study is similar to the phase precession reported in previous studies using single-cell recordings (Qasim et al., 2021). However, without simultaneous single-cell recordings, we can hardly study the relationship between the phase shifts and phase precession. With HGA as a proxy measure of spiking, future work studying both mechanisms could reveal insights into phase precession at different neural scale levels. All behaviors and neural signals in this study were recorded in individuals with drug-resistant epilepsy. Although we constrained our analyses to brain regions contralateral to the seizure onset zones and removed epileptic discharges under the guidance of an experienced epileptologist, it is possible, but unlikely, that the epileptic brain might alter the observed neural patterns.

STAR★METHODS

RESOURCE AVAILABILITY

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the Lead contact, Jie Zheng (zhengjie.may@me.com).

Materials availability

This study did not generate new unique reagents.

Data and code availability

Preprocessed data have been deposited at GitHub (https://doi.org/10.5281/zenodo.6970426) and are publicly available as of the date of publication. Analyses codes for key results of the paper were custom-written in MATLAB and are publicly available on GitHub (https://doi.org/10.5281/zenodo.6970426). Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Eight patients with drug resistant epilepsy (3 Male, 5 Female, Table S1) performed the task, while they were hospitalized at the University of California, Irvine Medical Center for surgical evaluation and were implanted with intracranial depth electrodes (Integra or Ad-Tech, 5-mm inter-electrode spacing). The electrode placements were exclusively guided by the clinical needs of localizing the seizure onset zone for possible surgical resections. Informed consent was obtained from each subject before the task and the research protocol was approved by the Institutional Review Board (IRB) of the University of California, Irvine. Subject selection was based on the following inclusion criteria: 1) normal IQ; 2) electrodes implanted in either the amygdala, the hippocampus, or the orbitofrontal cortex contralateral to or outside of the epileptogenic region were included for the analysis of neurophysiology signals.

METHOD DETAILS

We developed a decision-making task to study the processing of emotional sequences in humans (Figure 1A). Each trial started with a fixation period (Fixation; black screen with a white cross in the middle for 0.5s), followed by a context image display (Context; 1s) and a maintenance period with a blank screen (Maintenance; 0.5s). After the maintenance period, a neutral face image appeared (Face; 1s) and subjects were instructed to rate the valence of the face by mouse clicking to move the blue triangle on top of the color bar (Decision; self-paced). Warmer colors on the color bar denote more positive valence while colder colors represent more negative valence. The starting position of the blue triangle is always at the center of the color bar, that was informed to the subjects as a neutral rating position. Subjects were allowed to adjust their choices before clicking the “accept” button to lock in their decisions. Context images were selected from the International Affective Picture System (Lang et al., 2008) based on their valence scores (ranging from 1 to 9, unpleasant to pleasant), with 35 images each for different valence groups (valence scores, mean ± s.t.d.: negative = 3.87 ± 0.81; neutral = 5.23 ± 0.45; positive = 7.62 ± 0.63). The face images were all neutral faces, with 35 identities (female: 18; male: 17) selected from the NimStim Face Stimulus Set (Tottenham et al., 2009). Each neutral face was paired with 3 context images (1 negative, 1 neutral and 1 positive valence). The task was programmed in python using PsychoPy2 (version 1.8.02) software (Peirce et al., 2019). For each trial, subjects’ ratings and reaction times were recorded for further analysis.

Data collection

The task was presented on an Apple Macbook pro, which was placed on an overbed table at a comfortable distance in front of subjects. An external Apple keyboard was connected to the recording laptop to capture subjects’ responses. Trial onsets and offsets were detected by a photodiode attached at the bottom corner of the recording laptop’s screen (without blocking experimental stimuli). The photodiode signals (i.e., synchronization channel) and intracranial EEG signals were acquired simultaneously using a Nihon Kohden recording system (256 channel amplifier, model JE120A) with an analog-filter above 0.01Hz and a 5000Hz sampling rate.

QUANTIFICATION AND STATISTICAL ANALYSIS

Behavioral data processing

All the subjects completed the task (i.e., 105 trials). Due to the difference in subjects’ rating preferences (i.e., some subjects’ rating scale were more intensive while others’ rating range were more extensive), for each subject, we normalized the originally recorded ratings Rorig relative to the common neutral point (0.5), which subjects were informed as the middle of the color bar during the Rate stage:

Rmax / Rmin was the biggest/smallest rating score within each subject. The normalized ratings Rnorm were within the range of [0, 1] and maintained subjects’ rating polarity (negative vs. positive). For example, for a subject with face ratings ranging from 0.4 to 0.9, normalized against the common neutral point (0.5) will assign all the faces with ratings above 0.5 as positive. However, directly rescaling from [0.4, 0.9] to [0, 1] will falsely assign faces with original ratings between [0.5 0.65] as negative. We also rescaled the valence of context images from [1, 9] into [0, 1]. We then calculated the correlation between the normalized face valence and the normalized context valence using Pearson correlations (function corr.m from MATLAB Statistics and Machine learning Toolbox, Figure 1B).

Due to recording time constraints in the hospital setting, we didn’t independently capture subjects’ valence ratings of each face (i.e., without context image presented ahead). To ensure that the correlation was not driven by the potential bias toward the neutral face itself (i.e., subjects might perceive the neutral faces as positive or negative), we calibrated subjects’ face ratings paired with positive or negative context by subtracting their ratings of the same face paired with the neutral context:

The Pearson correlation between the calibrated face ratings R and context valence was also computed (Figure S1).

To quantify the strength of emotional bias attributed from the context valence (Figure 1C), we computed the contextual modulation strength for each trial based on the fraction of the face valence deviation (i.e., normalized face valence ratings subtract the neutral score, 0.5) over the context valence deviation (i.e., normalized context valence ratings subtract the neutral score, 0.5):

The absolute value of Equation (3) represents the strength of modulation, with bigger numbers indicating stronger valence influence carried from context images. The positive sign denotes consistent modulation (i.e., a face is rated as negative/positive when a negative/positive context image is presented ahead) while the negative sign represents opposite modulation (i.e., a face is rated as negative despite a positive context image presented ahead). Zero means no modulation.

Electrode localization

The electrode localization was performed using pre- and post-implantation structural T1-weighted 1mm isotropic MRI scans as well as post-implantation CT scans. For each subject, the post-implantation MRI and CT scans were registered to the pre-implantation scans using a 6-parameter rigid body transformation implemented with Advanced Normalization Tools (Avants et al., 2011). A high-resolution anatomical template (0.55mm) with labels of the amygdala, hippocampus and orbitofrontal cortex were aligned to the pre-implantation MRI scans using ANTs Symmetric Normalization. The electrode localization was guided by this co-registered anatomical template and was determined by examining the anatomical labels of the electrode artifacts from the co-registered pre- and post-implantation MRI scans. For visualizing electrodes across all the subjects (Figure 4D), we obtained MNI coordinates for all the electrodes by aligning each subject’s preoperative scan to the Colin27 brain template (Holmes et al., 1998) using a concatenation of an affine transformation followed by a symmetric image normalization (SyN) diffeomorphic transform (Avants et al., 2008).

Neurophysiological data preprocessing

The preprocessing of raw neurophysiology data was conducted using customized MATLAB scripts based on the open-source Toolbox Fieldtrip (Oostenveld et al., 2011). The neural recordings were first downsampled to 2000Hz and band-pass filtered between 1 and 250Hz using a zero-phase delay finite impulse response (FIR) filter with Hanning window. After demeaning the signal, we used Welch’s method to estimate its power spectral density (PSD). The line noises (usually 60Hz and its harmonics) inspected from PSD plots were removed using a discrete Fourier transform. The filtered signals were re-referenced to the nearest white matter electrodes on the same depth electrode probe based on the electrode localization results. The epileptiform discharges were manually inspected under the supervision of a neurologist (J.J.L.), who was blinded to electrode locations and trial information (e.g., stimuli onsets and subject’s performance). There were no seizures recorded in any subject while performing the task, and only the electrodes contralateral to or outside of the seizure onset zone were included in the analyses.

Task-induced power

The pre-processed data was FIR band-pass filtered between 1 and 250Hz through 35 logarithmically spaced frequencies. Then the Hilbert transform was applied to extract the analytic amplitude and instantaneous phase from each filtered time series. Next, we epoched the continuous data of analytic amplitude and instantaneous phase into trials. The trials overlapped with epochs that had been marked during the epileptiform inspection were removed from further analysis. For each trial, the task-induced power was calculated by squaring the analytic amplitudes and normalizing to the averaged pre-trial baseline power using the relative change in decibel conversion (dB). As shown in Figures S2A and S2B, the time-frequency representations were averaged across all the trials and all the electrodes within the same region of interest (amygdala, hippocampus, or orbitofrontal cortex).

Cross-regional phase synchrony

We computed phase synchrony between signals across each electrode pair within the same hemisphere (i.e., AMY-HPC, HPC-OFC, AMY-OFC). The phase synchrony was quantified using weighted phase lag index (WPLI) (Vinck et al., 2011) (window size = 500 ms, step size = 100 ms) sliding through the whole trial. Since WPLI is solely based on the imaginary component of the cross-spectrum between two signals, which is less sensitive to volume conduction or correlated noise source induced phase synchrony. Phase synchrony (i.e., WPLI) ranges from 0 to 1, with bigger numbers indicate stronger phase synchrony. For each frequency bin (10 bins between 2 and 16 Hz), we computed weighted phase lag index using the function implemented in the Fieldtrip (function ft_connectivityanalysisi.m). To statistically determine the significance of the observed phase synchrony, we used a non-parametric statistical approach (i.e., permutation test; Function ft_timelockstatistics.m from Fieldtrip Toolbox) that randomly permuted trial labels 1000 times (i.e., disrupt the correspondence between theta oscillations from two regions) to correct for multiple comparisons on a cluster level (i.e., across electrodes, across timepoints and across frequency bins). The phase synchrony was considered significant when it exceeded the 95% percentile of values in the null distribution.

Granger causality

To assess the directional influence between each cross-regional electrode pair (e.g., A and B), we computed spectral Granger causality (Geweke, 1982), which quantifies the prediction error of the signal in the frequency domain by introducing another time series. Before fitting to the multivariate autoregressive model to compute the spectral Granger causality, the time series data were low-pass filtered at 85 Hz, down-sampled to 250 Hz and normalized within each trial (e.g., subtracting the temporal mean and cross-trial mean). Then, we defined the model order using the Multivariate Granger Causality (MVGC) Toolbox based on the Akaike information criterion (Barnett and Seth, 2014). The model order for each subject varied from 11 to 14. The Granger Causality index was computed separately for each task stage (fixation, context, maintenance, face, and decision) for both directions (A to B, B to A). To access the significance of Granger Causality, we implemented the same method as described in the section ‘Cross-regional phase synchrony’ by randomly shuffling the trial labels 1000 times and the significance threshold was the 99th percentile of the surrogated data. As shown in Figure 2C, electrodes demonstrating significant Granger causality influence are connected, with the color denoting the origin of the directional influence and the thickness indicating the significant level of the directional influence.

Phase amplitude coupling (PAC)

We first computed the theta-gamma PAC (averaged within each trial across all task stages) across all the cross-regional electrode pairs. To avoid registering spurious manifestation of PAC, for both electrodes in each electrode pair, we computed the PSD of its original signal and the PSD of the envelope of filtered high gamma activity at the electrode-level within the low frequency range (1– 32 Hz). We then utilized an iterative algorithm to subtract the aperiodic background and looked for maximum (peak) of the residual signals. A peak is defined as point above the noise threshold (calculated from the standard deviation of the residual signal) and the putative peak frequency is determined by fitting the residual signal to a Gaussian (Donoghue et al., 2020). Electrode pairs with peaks detected within the theta range (4–12Hz) in all four PSD plots were included in the PAC analyses. For electrode pairs that met our “peak” criteria, we then computed theta-gamma PAC (averaged within each trial across all task stages) and a surrogate distribution of theta-gamma PAC by shuffling trial labels for 100 times. For electrodes with significant theta-gamma PAC (i.e., real theta-gamma PAC exceeds 95% of the null distribution), we calculated the PAC between each cross-regional electrode pair (window size = 500ms, step size = 100ms) sliding through the entire trial (Samiee and Baillet, 2017). PAC was calculated as the Euclidean norm of the summed polar vectors ( as the real component, Ø as the imaginary component) and was normalized to power to minimize the influence of signal magnitude in the measurement of coupling strength. The PAC strength was then Z score normalized to the baseline period and was averaged across all the frequency bins within the high gamma frequency range (70–250Hz). To statistically determine the significance of the observed PAC, we used similar non-parametric permutation test that randomly permuted trial labels 1000 times (i.e., disrupting the correspondence between high gamma activity and theta oscillations) to correct for multiple comparisons on a cluster level (i.e., across electrodes and across timepoints). PAC was considered significant when it exceeded the 95% percentile of values in the null distribution. We further averaged PAC across each task stage (fixation, context, maintenance, face, and decision) and over all pairs within the same region. In other words, if one subject has multiple HPC-OFC electrode pairs, we average PAC across all the HPC-OFC electrode pairs within each trial. We then correlate these averaged PAC with the contextual modulation strength (Figures 3C–3F, S3A, and S3B). To evaluate the significance of correlation, correlation coefficients were recomputed 1000 times after randomly assigning the contextual modulation strength to trials (i.e., disrupting the correspondence between PACs and the contextual modulation strength) to create null distributions. Correlation coefficient was considered significant if it exceeded 95% percentile of values in the null distribution.

High gamma activity burst

We used similar methods described in the previous study (Lundqvist et al., 2018) to detect bursting behavior within the high gamma range (70 –250 Hz). First, the preprocessed data within each trial period was frequency decomposed with the 1 Hz frequency resolution (using multi-taper approach adopted with frequency-dependent window lengths, with 4–8 oscillatory cycles and frequency smoothing as 0.2 of the central frequency f0, i.e., f0 ± 0.2 f0) and 1 ms temporal resolution to obtain smooth temporal profiles of power within the high gamma range (70–250 Hz). For each frequency bin within the high gamma range, high gamma activity bursts were defined as the time windows when the instantaneous power exceeded the two standard deviations above the mean power over the past 10 trials and lasted for at least three oscillatory cycles of the central frequency. Consistent with (Lundqvist et al., 2018), we created a local time-frequency map around the neighborhood of each high gamma activity burst with a two-dimensional Gaussian function to estimate its temporal duration (i.e., time interval where power was higher than half of the maximum power) and frequency bandwidth (i.e., frequency range where power was higher than half of the maximum power). The occurring phase of high gamma activity was defined as the time when the maximum power present in the detected high gamma activity bursts relative to theta oscillations.

Phase shift

We assessed phase shift in all HPC-OFC and OFC-AMY electrode pairs (see pair numbers in Table S1) when transitioning from the maintenance period to face presentation (i.e., [1, 2.5] seconds relative to the trial onset). For each HPC-OFC electrode pair, we quantify the relationship between the timing of OFC high gamma burst relative to HPC theta phase (i.e., burst phase) and the time by computing circular-linear correlation coefficients (Function circ_corrcl.m from MATLAB Circular Statistics Toolbox (Berens, 2009)). Notably, theta oscillations in humans are transient and often non-sinusoidal, which could lead to phase shifts when using Hilbert-transform on a narrow-band pass signal to extract phase information. Therefore, as a complementary approach, we also used waveform-based approach (Belluscio et al., 2012; Cole and Voytek, 2019; Dvorak and Fenton, 2014; Eliav et al., 2018) for phase estimation, which ensures the real peaks and troughs and can capture the non-sinusoidal characteristics of the signal. Similar results were observed regardless of the phase estimation methods (Figures S4E–H and Table S3). To access the significance of phase shift, for each electrode pair, we randomly assigned phases from the distribution of all burst phases to each high gamma burst and recompute the circular-linear correlation coefficients between the scrambled burst phase and the time. The same procedure was repeated 1000 times to form a null distribution of circular-linear correlation coefficients, which controlled for any effect of spurious phase estimates by disrupting the correspondence between burst phases and time. An electrode was considered exhibiting significant phase shift only if its circular-linear correlation value (i.e., negative correlation) is less than 95% of the null distribution. Similar analyses were performed on OFC-AMY electrode pairs as well except computing circular-linear correlation coefficients between the timing of AMY high gamma burst relative to OFC theta and the real time.

Bayesian decoding

For each electrode, we formed a three-dimensional binary temporospectral event matrix (n_frequency × n_timebin × n_trials), with 1 indicating a high gamma activity burst and 0 assigned to the rest of the matrix values, as the inputs for Bayesian decoding analyses. A memoryless Bayesian decoding algorithm was used to estimate the information of time and trial information (i.e., context valence and face valence). To decode the time information, the posterior probability of time relative to the trial onset was estimated given high gamma activity burst information (HGA, occurring time of high gamma activity for Figure 5A and occurring theta phase of high gamma activity bursts for Figures 5B and 5E) as:

The prior probability P(time) is uniform distribution by design. P( | time) is the distribution of , which is a normal distribution (Scherberger et al., 2005). Then the posterior probability of the time can be estimated as: where C(HGA) is a normalization factor that ensures the sum of all probabilities to be one. The same method was applied for the decoding of the trial information (i.e., context valence and face valence) with the occurring phase of HGA in the ROI highlighted in Figures 5B and 5E as inputs. To evaluate the decoding significance of trial information, the same decoding procedure was performed 1000 times using the surrogated data by shuffling the trial labels (context valence or face valence). The significant threshold (dashed horizontal line in Figures 5C, 5D, 5F, and 5G) was the 99th percentile value of the null distribution.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Pre-processed data	This paper	https://doi.org/10.5281/zenodo.6970426
Software and algorithms
MATLAB 2019b	MathWorks	RRID: SCR_001622
Python	Python.org	RRID: SCR_008394
Psychopy2	Peirce, 2009	https://www.psychopy.org
ANTs	Led by Brian Avants; Github	http://stnava.github.io/ANTs
Fieldtrip 20,220,206	Oostenveld et al. (2011)	http://www.fieldtriptoolbox.org/
Custom MATLAB Code	This paper	https://doi.org/10.5281/zenodo.6970426