Introduction

The brain continuously integrates information from different sensory systems enhancing its ability to filter noisy signals, to construct coherent percepts and guide appropriate actions (Stein and Meredith, 1993; Alais et al., 2010; Green and Angelaki, 2010; van Atteveldt et al., 2014). Such processes require the flexible and coordinated orchestration of distributed neuronal populations within and across brain areas (Singer and Gray, 1995; Engel et al., 2001; Fries, 2005; Senkowski et al., 2008; Siegel et al., 2012). Numerous studies have shown that this orchestration relies mostly on three modes of intrinsically generated coupling: synchronization of neural phase, the correlation of amplitude envelopes and the coupling between phase and amplitude (Siegel et al., 2012; Raichle, 2010; Canolty and Knight, 2010; Engel et al., 2013; Galindo-Leon et al., 2019). In addition, these coupling modes operate in multiple frequency bands and scales, mediating different cognitive functions (Engel et al., 2001; Fries, 2005; Siegel et al., 2012; Buzsaki, 2006; Jensen and Mazaheri, 2010; Varela et al., 2001; Engel and Fries, 2010), and providing different channels for long-range communication (Bastos et al., 2015; Michalareas et al., 2016). Phase coupling, in particular, has been suggested to serve the routing of information through cortical networks and to promote selective communication between distant brain areas (Engel et al., 2001; Fries, 2005; Womelsdorf and Fries, 2007), as well as to mediate stimulus selection and prediction (Engel et al., 2001; Fries, 2005; Engel and Fries, 2010; Michalareas et al., 2016; Lima et al., 2011). However, how phase coupling supports the brain communication across scales, and how it relates to behavior, is not well understood.

A natural scenario to test the role and mechanisms of phase coupling in brain dynamics at multiple scales are multisensory tasks. In the human brain, studies of multi-scale network dynamics have remained challenging although substantial advances have been made using non-invasive methods such as EEG or MEG. In particular, only a few studies have addressed the relation between phase coupling and multisensory processing (Keil et al., 2014; Giordano et al., 2017). In animal studies, the potential role of functional connectivity involved in multisensory processing has so far been investigated only by multielectrode recordings from the same cortical areas (Lakatos et al., 2007; Kayser et al., 2008) or by simultaneous recordings from at most two different regions (Ghazanfar et al., 2008). Short-range neural interactions are well established at the level of microcircuits (Coulter et al., 2011; Bastos et al., 2012; Feldmeyer et al., 2013). Invasive studies on functional connectivity between distant population in behaving animals typically have recorded simultaneously only from relatively small numbers of sites (Womelsdorf and Fries, 2007; Gregoriou et al., 2009; Miller and Buschman, 2013; Liebe et al., 2012; Salazar et al., 2012).

Electrocorticographic (ECoG) arrays offer an option to study simultaneously the relation between multiple-scales of network dynamics and behavior (Bastos et al., 2015; Bockhorst et al., 2018; Rubehn et al., 2009; Stitt et al., 2017). Because of their bio-compatibility and recording stability, they have been used for the study of functional connectivity in networks underlying cognitive and sensorimotor functions (Crone et al., 2006; Bosman et al., 2012; Keller et al., 2014; Fukushima et al., 2015; Lewis et al., 2015) in humans (Canolty and Knight, 2010; Fox et al., 2018), monkeys (Bastos et al., 2015; Bosman et al., 2012) and rodents (Toda et al., 2018). Furthermore, small-sized ECoG electrodes can provide local information at the level of cortical columns (Bockhorst et al., 2018; Baratham et al., 2022) which, when distributed across distant regions in an array, provide a unique opportunity to study the interaction between local and global cortical dynamics.

Here, we investigated how changes in functional connectivity and power, indicative of global and local dynamics, respectively, determine the detection of multisensory stimuli. We hypothesized that in a multisensory detection task the functional coupling between modalities, rather than local power, determines the performance. We tested this hypothesis in a large-scale cortical network involving visual, auditory, somatosensory and parietal areas using a custom-made 64-channel ECoG array in behaving ferrets. Animals were trained in a 2-alternative-forced-choice paradigm to detect brief auditory and visual stimuli presented either left or right from the midline. A key finding of our study is that phase coupling of visual and auditory regions with parietal cortex was predictive of task performance. Our results suggest that fluctuations in the networks state, particular with respect to long-range connectivity, are decisive for the response of the upcoming task during multisensory processing.

Materials and methods Animals

All experiments were approved by the independent Hamburg state authority for animal welfare (BUG-Hamburg) and were conducted in accordance with the guidelines of the German Animal Protection Law. Data were recorded in four adult female ferrets (Mustela putorius furo; Euroferret, Dybbølsgade, Denmark) aged 2 and 4 years (n = 2 each). Animals were individually housed in standard ferret cages equipped with an enriched environment, under controlled conditions (21°C, 12-h light/dark cycle, with lights on at 8:00 a.m.). The animals had ad libitum access to food pellets, while access to tap water was restricted for 8 h prior to the experiments and training procedures. All behavioral testing was conducted during the light cycle, specifically between 10:00 a.m. and 2:00 p.m. The animals were treated with a deslorelin acetate (4.7 mg implant) chip to control their ovarian activity (Proháczik et al., 2010). Animals were trained on a two-alternative forced-choice task, described in detail below, involving lateralized detection of visual and auditory stimuli that we have established previously (Hollensteiner et al., 2015). After completion of all recording sessions (for details, see below), animals were euthanized in deep anesthesia with a lethal dose of pentobarbital (400 mg/kg i.p.). Subsequently, animals were transcardially perfused with fixative (4% paraformaldehyde in phosphate-buffered saline) and brains were removed for histological processing.

Sensory stimuli

All experiments were carried out in a dark sound-attenuated chamber (Acoustair, Moerkapelle, Netherlands). Visual and auditory stimuli were generated with the Psychophysics Toolbox (Brainard, 1997). Visual stimuli were presented on an LCD monitor (Samsung SyncMaster 2233, frame rate 100 Hz) placed 20 cm in front of the animal. Auditory stimuli were generated digitally with a sample rate of 96 kHz and delivered through two Beyerdynamic T1 speakers. Auditory stimuli consisted of white noise bursts of 100 ms duration, including 5 ms cosine rise/fall, presented with intensities between 2 and 62 dB SPL. Visual stimuli consisted of circular gratings (22.5°, 0.2 cycles/°, 5 Hz) with Michelson contrast (Cm) between 0 and 0.38. A static random noise pattern located in the center of the screen was presented to indicate trial onset by a decrease in contrast (not overall luminance). Multisensory stimuli consisted of both stimulation types presented simultaneously.

Training

Initially, ferrets were handled and accustomed to a housemade flat-bottomed tube, which served as the animals’ enclosure during the experiment. Figure 1A shows a schematic of the experimental setup. Subsequently, animals were trained in the spatial detection task, during which the unimodal as well as crossmodal psychophysical thresholds were estimated. In the task, we used a crossmodal approach that bears on established unimodal training paradigms (Nodal et al., 2008). Animals were restricted from access to water for a period of 8 h before the measurements, and conditioned by using water rewards during the task. In the first phase of the study, the animals were trained to detect unimodal auditory and visual stimuli presented in a randomized fashion. Auditory and visual unimodal detection thresholds were determined using 20 different stimulus amplitudes (auditory: 2–62 dB SPL; visual: 0–0.38 Cm) in a 1 down/3 up staircase procedure (Kaernbach, 1991). Next, unimodal and bimodal thresholds were assessed in a combined approach, using the previously determined unimodal thresholds to adjust the test parameters (for details of the procedure see Hollensteiner et al., 2015). Subsequently, the ferrets were accustomed to electrophysiological recordings during the detection task with a reduced set of stimulus amplitudes (eight per modality), adjusted to the individual psychometric functions to acquire a higher number of trials in the performance range of interest around 75% accuracy.

Figure 1. Experimental setup and lateralized detection task. (A) Schematic of the experimental setup in a top view: (a) LCD-screen, (b) speakers, (c) three light-barrier-waterspout combinations (left, center, right; the red dot indicates a broken light-beam), (d) pedestal. (B) Sequence of events in a single trial: (I) inter-trial window, (II) baseline window, (III) delay window, (IV) stimulus window, and (V) response window. The three circles below each frame represent the state of the light-barriers (white = unbroken, red = broken). The center of the screen displays a static circular random noise pattern.

Detection task

The trial schedule of the detection task is shown in Figure 1B. To initialize a trial, the ferret had to maintain a central head position to break the central light-barrier for 500 ms. This caused the static random noise pattern in the center of the screen to decrease in contrast informing the animal that the window for stimulus presentation (from 0 to 1,000 ms after trial onset) had started. During this interval, the animal had to further maintain its central head position. A stimulus was presented for 100 ms on either the left or the right side. Stimulus presentation could be unimodal visual (“V”), unimodal auditory (“A”) or a temporally congruent bimodal combination of an auditory and a visual stimulus presented on the same side. This combination either consisted of a variable auditory stimulus accompanied by a visual stimulus of constant contrast (“Av”) or, conversely, a visual stimulus of varying contrast supported by an auditory stimulus of constant amplitude (“Va”). The intensity value of the accompanying second stimulus with constant amplitude was set at a level of 75% accuracy. After stimulus offset, the animal had to respond within 700 ms by moving its head to the stimulated side; otherwise, the trial was considered as a miss (no response). If the response was correct the animal received a water reward (~80 μL) from a spout at the stimulus position and could immediately start the next trial. If the response occurred prematurely (before stimulus onset or within 100 ms after stimulus onset), was incorrect (wrong side) or omitted, the trial was immediately terminated, followed by a 2,000 ms interval during which no new trial start could be initialized.

ECoG array implantation

Micromachining technology (Rubehn et al., 2009) was used to design and implement an ECoG array that matched the anatomy of the ferret brain (Figures 2A,B). The thin-film (10 μm) polymide-foil ECoG contained 64 platinum electrodes with a diameter of 250 μm each, arranged in a hexagonal formation at an inter-electrode distance of 1.5 mm (Figures 2C,D).

Figure 2. ECoG recordings from the ferret brain. (A) Lateral view of the left hemisphere. Lines depict sulci in the posterior part of the hemisphere. (B) Anatomical organization of the posterior part of the ferret brain (adopted from Bizley et al., 2007). (C) Schematic of the implanted ECoG (contact spacing: 1.5 mm; ø: 250 μm). (D) In-situ picture of the ECoG array over the left posterior hemisphere during implantation. LS, lateral sulcus; SSS, suprasylvian sulcus; S2, secondary somatosensory cortex; PPr/c, rostral/caudal posterior parietal cortex; 17, 18, 19, 20, 21, early and higher-order visual cortex; SSY, suprasylvian field; A1, primary auditory cortex; AAF, anterior auditory field; ADF, anterior dorsal field; PPF, posterior pseudosylvian field; PSF, posterior suprasylvian field. (B–D) Adapted from Stitt et al. (2017).

The surgical procedure for the implantation of the ECoG array started with the induction of general anesthesia by using a ketamine-medetomidine mixture (10 mg/kg and 1 mg/kg, respectively). During surgery, anesthesia was maintained by ventilating the animal with isoflurane (1–1.5%) in combination with 70% N20 and 30% O2. To monitor the state of anesthesia, physiological parameters such as the electrocardiogram (ECG) and rectal temperature were monitored throughout the surgery. All surgical procedures were performed under sterile conditions. After the operating field was prepared, a craniotomy was performed using a saline-cooled ultrasonic microsaw (Mectron) in order to expose the posterior half of the left cerebral hemisphere. The dura was carefully removed and the ECoG array was gently placed on the surface of the cortex such that it covered occipital, temporal and parietal areas (Figure 2D). The dura was then folded back over the ECoG array and an absorbable artificial dura was placed above the ECoG, covering the trepanation to full extent. The excised piece of bone was fixed back in place with titanium plates and screws and subsequently the gaps were filled with fast set absorbable bone putty. Finally, the ECoG’s connector was placed on top of the skull and fixed in place with titanium screws and dental acrylic. After the surgery, the animals received analgesic (carprofen, 2.5 mg/kg) and antibiotic (enrofloxacin, 5 mg/kg) medication. The animal was allowed a recovery period of a minimum of 8 weeks post-surgery. Throughout this post-operative phase, body weight and behavioral patterns were closely monitored. Regular handling sessions were maintained to ensure the animal remained habituated to the experimental setup. To align electrode positions across animals, the ECoG array placement and the cortical parcellation introduced by Bizley et al. (2007) was utilized (Figure 2B). For each animal, the exact ECoG array position over the posterior cortex was photographically documented during surgery. The locations of all 64 ECoG electrodes were then projected onto a scaled illustration of a ferret model brain. Each electrode was then mapped to the underlying cortical area.

Electrophysiological recording and preprocessing

Local field potentials (LFPs) were digitized at 1.4 kHz and sampled simultaneously with a 64-channel AlphaLab SnR recording system (Alpha Omega Engineering, Israel). The high pass filter was set at 0.1 Hz and the low pass filter at 357 Hz.

To ensure comparability across animals and modalities, only recording blocks with an accuracy of 75 ± 10% were considered in all electrophysiological analyses. Furthermore, this selection of trials ensured better comparability between all four modalities because the value that was fixed during the bimodal stimulation was set at a level of 75% accuracy for all stimulation amplitudes. To assess a sufficient amount of trials with stimulus presentation contralateral to the implanted hemisphere per modality in the 75% accuracy range, trials collected on different days were pooled; (see Hollensteiner et al., 2015) for non-stationarity effects across sessions.

All offline data analysis was performed using custom scripts in Matlab (The Mathworks Inc., MA). The hardware-filtered LFP signals were band-pass filtered with a phase-preserving fourth-order Butterworth filter between 2 and 200 Hz. Next, band-stop Butterworth filters at 49–51, 99–101 and 149–151 Hz were applied to notch the line noise. Subsequently, the continuous LFP signals were cut into epochs aligned to trial, stimulus and response onset, respectively. In each of these analysis time windows data were cut from 500 ms pre- to 500 ms post-onset, with 500 ms real data padding on each side to prevent edge artifacts in frequency space. Afterwards, we applied independent component analysis (ICA) (Hyvärinen and Oja, 2000) to the concatenated analysis window data in order to detect and remove ECG, muscle and eye blink artifacts. Next, we re-referenced the LFP time series of each ECoG contact with the LFP time series of its next posterior neighbor. This processing step created 55 virtual electrodes from the 64 recorded contacts. The new virtual electrode was located in the midpoint between both real electrode positions. Subsequently, the oscillatory signal components of unimodal and bimodal trials were analyzed using spectral decomposition.

Spectral power analysis

Channel-wise spectral power was computed by taking the square of the absolute value of the Fourier transform. Spectra were computed for all three analysis time windows using a Hanning window approach (2–200 Hz, 2 Hz steps, 500 ms window; for the LFP spectrogram the window was shifted from –500 ms to 500 ms around the window of interest) in 1 ms steps.

Quantification of phase locking

To estimate phase synchronization between ECoG signals, we calculated the phase locking value (PLV) (Lachaux et al., 1999) across all frequencies (2–200 Hz). In general, the instantaneous phase θ was extracted from the analytic signals that were produced by the Fourier transform of the convolution of the ECoG time series with the Hanning window. The PLV between channels A and B at carrier frequency f is defined as follows (Equation 1):

PLVABf=1N|∑n=1NeiθnAf−θnBf|    (1) For an event related paradigm the circular average of phases is in general non-vanishing, and a non-vanishing result using the classical definition of PLV does not imply statistical dependence between phases in two electrodes. Therefore, we corrected for evoked effects and defined the PLV here as Equation 2 PLVABf=eiθAf−θBf−eiθAfeiθBf    (2)

where 〈〉 denotes average over trials. For an infinite number of trials, the PLV defined here vanishes exactly if the phase in electrode A and B are statistically independent. To contrast PLV between modalities, we first computed the global mean PLV for the same time windows selected in the spectral power analysis. Subsequently, the PLV between anatomical regions, as defined in Figure 2B, was calculated by averaging across all virtual electrodes overlaying the same area. To normalize and compare PLV values from hit and miss trials we computed the sensitivity index (d′), which is defined by Equation 3:

d′=μhit−μmissσhit2+σmiss2    (3)

where μ denotes the mean and σ the standard deviation for hit and miss trials, respectively.

Principal component analysis of functional connectivity

We used a similar approach to that of Leonardi et al. (2013) to identify hidden sub-networks during the baseline period that may be associated with performance by applying principal component analysis (PCA) to trial-based connectivity matrices. To facilitate the contrast in the analysis we separated matrices in groups of hits and misses. Before applying PCA, the global mean connectivity matrix was subtracted from matrices of all trials within each group before concatenating all recording sessions. To avoid redundancy we took only the lower triangular part of each nxn connectivity matrix (excluding elements of the diagonal) and then vectorized it to obtain a vector of n(n − 1)/2 elements. The vectors of all animals were concatenated along the trial dimension into an array of (n − 1)/2 × T elements, with T being the total number of trials. PCA was then performed on the resultant matrix. Eigenvectors or principal components (PCs) can be understood as features or subnetworks that characterize the variation across the collection of correlation matrices and represent network patterns. Because these capture connectivity patterns they are also called eigenconnectivities (Leonardi et al., 2013). Finally, to reduce dimensionality and obtain the main eigenconnectivities, matrices were thresholded by the lowest and strongest 10% with respect to their mean.

Directionality analysis

To estimate the direction of causal influences we calculated the phase slope index PSI (Ψ) (Nolte et al., 2008). PSI is highly robust against false estimates caused by confounding factors of very general nature, and is defined as Equation 4:

Ψ˜ij=imag∑f∈FCijfCijf+δf    (4) where Cijf=Sijf/SiiSjj is the complex coherency, S the cross-spectral matrix, δf the frequency resolution. F is the set of frequencies over which the slope is summed. The complex coherency C is a Hermitian matrix and consequently PSI is antisymmetric and, hence, directional. Intuiti vely, PSI is the slope (as a function of frequency) of the phase difference between two signals averaged within a frequency band. The slope is larger, when the time lag between the signal and the coupling increases. The sign depends on which signal is earlier. It is constructed such that it vanishes (apart from stastical fluctuations) if the observed signal is an instantaneous superposition of independent sources, and the average across frequencies is weighted with coupling strength. To normalize Ψ˜ we divided it by its standard deviation (Equation 5) Ψij=Ψ˜stdΨ˜    (5) estimated by the jackknife method. The acquired PSI matrices were thresholded equivalently to the eigenconnectivity matrices and their sign was extracted. Results

To investigate neural responses and task-related fluctuations of functional connectivity, we recorded local field potentials (LFPs) via 64-channel ECoG arrays chronically implanted in four adult female ferrets.

Multisensory effects on behavior

Before the electrophysiological recordings were conducted, a psychometric evaluation was performed on each animal to determine the stimulus parameters (see Methods). Figures 3A,B show examples of psychometric functions fitted to the behavioral data for the A and Av, and for the V and Va conditions, respectively. For further analysis, data were taken only from trials with stimulus intensities at threshold to allow comparisons across conditions and animals (Figures 3A,B).

Figure 3. Trial selection for electrophysiological analysis and corresponding RTs. (A) Example data for task performance from one ferret. The plot shows psychometric functions for the uni- and bimodal condition. Stimulus response function fitted to performance data for the eight different stimulus intensities for the unimodal auditory condition (A) and the condition where the variable auditory stimulus was accompanied by a visual stimulus of constant contrast (Av). The dot diameter indicates the number of trials at a given stimulus intensity. The green ellipse indicates the stimulus intensity range around 75% performance from which trials were selected for subsequent analysis. The unmasked parts of the graphs indicate the range of the actually tested stimulus amplitudes. (B) Data from the same animal as shown in A for the unimodal visual condition (V) and the condition where the variable visual stimulus was accompanied by an auditory stimulus of constant intensity (Va). (C) Reaction times (mean ± standard deviation) for the four different conditions for the selected trials. The four different shapes (circle, diamond, cross, star) represent the single trial RT of each animal. Error bars indicates the mean RT ± STD across animals. Black asterisks indicate significant differences between modalities or pooled data from bimodal (Va and Av) and unimodal (A and V) stimulus pairs (brackets).

Reaction times (RT) have been broadly used to quantify differences between unimodal and multisensory stimulation. A t-test was run to investigate the differences between bi- and unimodal conditions. Responses to bimodal stimulation were, on average, faster than those to unimodal stimulation (p < 0.05). To evaluate whether this effect was driven by the auditory or the visual condition, we computed a one-way analysis of variance (ANOVA) on the mean RT of each animal and condition, with sensory modality as a main factor. This revealed a main effect of condition [F (3, 12) = 5.3, p < 0.05], but post hoc t-tests only revealed significant differences between the A-trials and both bimodal conditions. The population averages of RTs for Av, Va and V were not significantly different (Figure 3C). Nevertheless, the population mean RT of V (0.207 s) was slightly higher than both, Va (0.194 s) and Av (0.197 s). The significant difference between the RTs in the A and the Av condition suggest a multisensory interaction effect leading to reduction in RT for the bimodal stimulation. However, our data do not suggest a substantial improvement of visual detection by a concurrent auditory stimulus.

LFP power reflects stimulus processing and response preparation

A major aim of our analysis of the electrophysiological data was to test the relation of power and functional connectivity to task performance by contrasting hit and miss trials. In the hit trial group, we included all trials with stimulus intensities around 75% threshold and stimulation contralateral to the implanted ECoG (265 ± 14 hit trials per animal). To match the number of hit trials we used all miss trials (236 ± 45 miss trials per animal) throughout all sessions (15 ± 2 sessions per animal). Miss trials were defined as those with sensory stimulation that did not elicited a behavioral response, i.e., the ferret maintained the centered head position throughout the response window. The number of false responses, i.e., orientation of the ferret to the side contralateral to the sensory stimulation, were too few and highly variable across animals (104 ± 58 false responses per animal) and, thus, were not considered for analysis of the electrophysiological data.

We checked whether the spectral characteristics of LFPs immediately before task execution predicted the animal’s performance. To assess how LFP spectral characteristics evolve during hit, miss and bi- and unimodal stimulation, we computed the power across time windows aligned to baseline, stimulus presentation and response onset for each animal individually (Figure 4). Results of all animals were pooled after correction for different trial numbers. Power differences between hit and miss trials were analyzed in the theta (4–8 Hz), alpha (8–16 Hz), beta (18–24 Hz), gamma (30–80 Hz) and high-gamma (80–160 Hz) frequency band across the analysis time windows (Supplementary Table S1). In addition, one-way ANOVAs were calculated on power values within each frequency band with condition as the main factor to examine differences related to stimulus type (Supplementary Figure S1 and Supplementary Table S1).

Figure 4. Spectral power analysis of pooled hit and miss trials. (A) Grand average power spectra for all hit (green) and miss (red) trials ± standard error of the mean (SEM). For the baseline window absolute power is shown. Note that for the time windows around stimulus and response onset, the spectral changes shown are relative to the baseline window. Asterisks indicate significant differences between hits and misses within the specified frequency band (FDR corrected). Labels indicate theta, alpha, beta, gamma and high-gamma band, respectively. (B) Time-frequency representation of power in the three analysis time windows, expressed as change relative to baseline before trial onset for all hit and miss trials. The vertical line represents trial, stimulus and response onset, respectively.

Figure 4A illustrates the grand average spectra, while Figure 4B shows the time-frequency representation of power changes during each time window (baseline, stimulus and response onset) for both hit and miss trials. During the baseline, only the high-gamma frequency-band (80–160 Hz) showed moderated but significant differences, with higher power for miss trials (p < 0.05; FDR corrected). Comparison of the average spectra in the stimulus onset analysis window showed significant reductions for hit and miss trials in the gamma frequency-range (30–160 Hz) compared to baseline (Supplementary Table S1). However, hit trials showed significantly higher power compared to miss trials after stimulus onset in both gamma (p < 0.001) and high-gamma (p < 0.001) frequency bands (Figure 4A). Around response onset, hit trials consistently showed significantly higher power compared to miss trials in the theta- (p < 0.01), alpha- (p < 0.01), gamma- (p < 0.001) and high-gamma band (p < 0.001). Statistical comparison across time windows within hit- and miss-trials revealed significant decrease and increase in the gamma frequency range compared to baseline and stimulus onset, respectively (Supplementary Table S1; all p-values, within and across windows of interest, FDR corrected).

Supplementary Figure S2 shows the topographies of the respective spectral differences between hit and miss trials. In the baseline window, power topographies showed high spatial uniformity across all frequency bands. Around stimulus onset, notable differences occurred mainly in the beta band where, in hit trials, higher beta power differences occurred in parietal areas. The highest regional differences occurred in the response window. Occipital regions exhibited higher power across all frequency bands during response onset in miss trials compared to hit trials. In contrast, auditory and parietal areas showed increased power in lower frequency bands (theta and alpha) for hit trials. An increase of gamma band power occurred in the response window for hit trials in regions around the lateral sulcus (Supplementary Figure S2).

To investigate stimulus-type dependent effects within hit-trials a one-way ANOVA between the four conditions (Av, Va, A and V; Supplementary Figure S1 and Supplementary Table S1) was calculated. It showed no main effect in the factor condition across frequency bands during baseline, stimulus or response onset. This suggests that power changes related to stimulus processing and response preparation did not depend on the modality of the presented stimulus or on crossmodal interactions.

Functional connectivity predicts performance

To determine whether functional connectivity predicts the animals’ performance we computed the PLV between all pairs of ECoG channels for the different frequency bands. PLV analysis corrected for evoked effects was performed for the same time windows as the power analysis: aligned to baseline, stimulus and response onset. Figure 5 displays the population average PLV spectra for all hit and miss trials as well as the difference between hits and misses (Figure 5B) and the relative spectral changes in the stimulus and response onset window relative to baseline (Figure 5B). Paired t-tests were applied between hit and miss trials within and across all time windows for each frequency band separately. Within each analysis window, significantly higher PLV was observed for hits compared to miss trials in all frequency bands, with exception of the theta band during response onset (p < 0.05, FDR corrected; Figure 5A).

Figure 5. Grand average of functional connectivity for hit and miss trials. (A) Grand average phase locking value (PLV) for all hit (green) and miss (red) trials for the analysis time windows (baseline, stimulus and response onset). (B) Differences between hits and misses for the three analysis windows and changes of PLV around stimulus and response onset, expressed as changes relative to the baseline window. Asterisks indicate significant differences between hits and misses within the specified frequency band (p < 0.05, FDR corrected). Labels in A denote theta, alpha, beta, gamma and high-gamma band, respectively.

When compared to the baseline time window, PLV was significantly reduced in the stimulus and response onset windows in nearly all frequency bands. Exceptions were in the theta frequency band for miss trials around stimulus as well as response onset.

A one-way ANOVA within analysis time windows and frequency bands with stimulus condition as the main factor revealed significant differences in PLV between stimulation conditions (Av, Va, A and V; Supplementary Figure S3). During the baseline period there were no significant effects in the main factor condition in the theta frequency band [F (3, 44) = 2.46, p > 0.05]. However, starting with the alpha [F (3, 76) = 4.98, p = 0.003], beta [F (3, 60) = 3.71, p = 0.016], gamma [F (3, 412) = 31.89, p < 0.001] and high-gamma [F (3, 652) = 56.99, p < 0.001] frequency bands, there was a main effect of conditions. This effect persisted also for the windows around stimulus onset [theta: F (3, 44) = 1.23, p > 0.05; alpha: F (3, 76) = 3.74, p = 0.015; beta: F (3, 60) = 3.5, p = 0.021; gamma: F (3, 412) = 26.12, p < 0.001; high-gamma: F (3, 652) = 55.67, p < 0.001] and response onset [theta: F (3, 44) = 2.31, p > 0.05; alpha: F (3, 76) = 3.63, p = 0.017; beta: F (3, 60) = 3.67, p = 0.017; gamma: F (3, 412) = 25.12, p < 0.001; high-gamma: F (3, 652) = 51.24, p < 0.001]. Furthermore, post-hoc t-test confirmed a constant pattern of significant differences between conditions (Supplementary Table S2). In particular, for hit trials connectivity was higher in the A and Av conditions than in the V condition, and there was a trend for connectivity in bimodal Va trials to be higher than in unimodal V trials (Supplementary Figure S3).

Taken together, in contrast to the power, the PLV revealed significant differences between hit and miss trials. These differences occurred over a broad range of frequencies across all analysis time windows. These results suggest that large-scale phase-fluctuations in network state are determinant to the animals’ performance.

Large-scale coupling shows specific changes during the task

To examine functional connectivity in relation to the topography of cortical areas, we assigned the data of each ECoG contact to a distinct region based on the functional map of cortical areas from Bizley et al. (2007) (see Methods) and constructed functional connectivity matrices accordingly (Figure 6). Cortical areas were grouped in three functional systems, comprising visual (areas 17, 18, 19, 20, 21), auditory (areas A1, AAF, ADF, PPF and PSF), and parietal areas (SSY, PPc, PPr, S2). Differences between hit and miss trials were expressed using the sensitivity index (d′).

Figure 6. Matrices of functional connectivity across anatomical areas. Cells show population average phase locking values (PLV) for the different frequency bands, expressed as d′ for hits vs. misses, for all analysis time windows. Abbreviations as in Figure 2.

In line with the grand average PLV results, the connectivity matrices generally displayed positive d′ values and, thus, higher functional connectivity for hit compared to miss trials (Figure 6). This effect was most prominent during the baseline window and decreased during the stimulus presentation. For the alpha, gamma and high-gamma band, the PLV difference between hit and miss trials increased again in the response window. Generally, connections between early sensory cortices and parietal cortex displayed the highest differences between hit and miss trials, in particular during the baseline window (Figure 6; Supplementary Figure S5).

Eigenconnectivity reveals subnetworks related to performance

We asked whether the trial-to-trial fluctuations in functional connectivity during the baseline period could affect the animal’s performance. We first identified the subnetworks that most strongly contributed to fluctuations during the baseline across trials. To find the intrinsic structure in the patterns of connectivity, we applied principal component analysis (PCA) to functional connectivity matrices for hits and misses along the trial dimension (see Methods). We focused on a subset of areas (18, 19, 21, PPc, PPr, A1, and AAF) from the recorded functional systems for the analysis of task-related subnetworks to reduce dimensionality.

Figure 7 shows eigenconnectivity patterns in the beta and gamma bands for hit and miss trials in the first and second PCs. Only those connections which had the smallest and the largest fluctuations of functional connectivity in the baseline window are displayed (complete matrices and explained variance for PCs 1–6 in all frequency bands are shown in Supplementary Figure S4). PC1 did not show substantial differences between hits and misses (Figure 7A; Supplementary Figure S4). For both hits and misses, the first PC revealed that the largest fluctuations occurred between auditory and visual areas, and the smallest between visual and parietal cortices. Importantly, PC2 revealed substantial differences between hit and miss trials. These differences were most prominent in the theta, alpha, and beta bands, but also occurred in the gamma band (Figure 7B; Supplementary Figure S4). In the beta band, a higher fluctuation in functional connectivity between AAF and visual areas was a common pattern in miss trials for PC2. A similar pattern was shown in the gamma band (Figure 7). In contrast, the highest connectivity fluctuations in the beta band for hits were observed within the visual system. From the third to the sixth components, the connectivity fluctuations revealed more complex patterns, with increases as well as decreases in both hit and miss trials (Supplementary Figure S4). The directionality in Figure 7 and Supplementary Figure S4, which was extracted from PSI matrices, represents the directionality of the coupling in subnetworks and does not necessarily carry information about directionality of the fluctuation.

Figure 7. Eigenconnectivity (EC) of hit and miss trials for the baseline window. Depicted are EC networks for beta and gamma of the first two PCs, which together explain 42 and 32% of the variance in the PCA for hits and misses, respectively. Green and red arrows represent hit and miss networks, respectively. Thick solid lines indicate the upper 10%, and thin dashed lines the lowest 10% of connections in the EC distribution. Arrowheads indicate the direction of information flow extracted by PSI. Colored circle segments represent the three functional systems (visual = magenta; parietal = blue; auditory = orange). Note that EC matrices for the first six PCs with all connection are shown in Supplementary Figure S4. Abbreviations as in Figure 2.

It should be noted that the eigenconnectivity patterns cannot be directly compared to the global connectivity differences between hits and misses described in the preceding sections because they reflect changes in subnetworks, relative to the average network coupling strength, in the two groups of trials. Importantly, our analysis shows that changes in connectivity do not occur arbitrarily across the network, but are distributed in subnetworks that show specific increases or decreases in phase coupling.

Between- and within-system interactions differ in their relation to behavior

The topography of functional connectivity changes was further analyzed by contrasting between-system and within-system interactions (Figure 8; Supplementary Figure S5). A previous study from our group showed that this is a suitable classification to determine global/local effects associated to changes in brain state (Stitt et al., 2017). For hit trials, interactions between the three functional systems were, on average, stronger than functional connections within the respective system. Regarding the within-system connectivity, the auditory system showed the lowest d′. In the period from baseline to stimulus onset, the strength of d′ decreased across all frequency bands. With response onset the spectral profile of d′ became more heterogeneous. In the theta-band d′ decreased further, whereas it remained unaltered in the beta-band, and increased in the alpha-, gamma- and high-gamma band.

Figure 8. Comparison of long-range (between-systems) and short-range (within-system) functional connectivity. (Left) Cortical regions were grouped into visual (magenta), parietal (blue), or auditory (orange) systems. PLV was averaged for all pairs of regions that were located within the same system (grey arrows) and pairs between systems (red arrows). Subsequently, d′ was calculated for hit and miss trial PLVs. (Right) Differences of between-region and within-region PLV (expressed as d′) are plotted as a function of carrier frequencies. Large values indicate higher between-system PLV during hit compare to miss trials. Symbols indicate significant differences between baseline and stimulus onset (*), baseline and response onset (o), and stimulus onset to response onset (+), respectively.

For statistical analysis, we averaged PLV d′ values and derived the difference of d′ values between long-range (between-system) and short-range (within-system) connectivity and calculated statistics across analysis time windows for each frequency band (Figure 8; for statistics across systems and within analysis time windows and frequency see Supplementary Figure S5; for condition specific statistics see Supplementary Figure S6). To reveal dynamic effects within frequency bands of interest a one-way ANOVA with time window as the main factor was computed. There was no main effect in the theta band [F (2, 6) = 4.39, p > 0.05]. The ANOVA exposed significant effects for the main factor time window alpha [F (2, 12) = 17.17, p < 0.001]; beta [F (2, 9) = 14.43, p = 0.002]; gamma [F (2, 75) = 3.71, p < 0.001]; high-gamma [F (2, 120) = 17.94, p < 0.001]. Post hoc t-tests revealed significant differences between baseline and stimulus onset for all frequency bands with the main effect (alpha: p = 0.011; beta: p = 0.002; gamma: p < 0.001; high-gamma: p < 0.001), between stimulus and response onset for the gamma frequency range (gamma: p = 0.009; high-gamma: p < 0.001) and between baseline and response onset in the alpha- (p < 0.001) and beta- (p = 0.006) frequency bands (Bonferroni corrected).

Overall, comparison of functional connectivity between and within cortical systems across the different analysis time windows revealed that between- compared to within-system connectivity was stronger during the task period compared to the pre-stimulus baseline. In hit-trials, phase coupling was higher between systems in the gamma frequency range during stimulus onset compared to baseline and response onset. Furthermore, low frequency (alpha and beta) coupling between systems increas