Imaging of neural oscillations with embedded inferential and group prevalence statistics

Magnetoencephalography and electroencephalography (MEG, EEG) are essential techniques for studying distributed signal dynamics in the human brain. In particular, the functional role of neural oscillations remains to be clarified. For that reason, imaging methods need to identify distinct brain regions that concurrently generate oscillatory activity, with adequate separation in space and time. Yet, spatial smearing and inhomogeneous signal-to-noise are challenging factors to source reconstruction from external sensor data. The detection of weak sources in the presence of stronger regional activity nearby is a typical complication of MEG/EEG source imaging. We propose a novel, hypothesis-driven source reconstruction approach to address these methodological challenges. The imaging with embedded statistics (iES) method is a subspace scanning technique that constrains the mapping problem to the actual experimental design. A major benefit is that, regardless of signal strength, the contributions from all oscillatory sources, which activity is consistent with the tested hypothesis, are equalized in the statistical maps produced. We present extensive evaluations of iES on group MEG data, for mapping 1) induced oscillations using experimental contrasts, 2) ongoing narrow-band oscillations in the resting-state, 3) co-modulation of brain-wide oscillatory power with a seed region, and 4) co-modulation of oscillatory power with peripheral signals (pupil dilation). Along the way, we demonstrate several advantages of iES over standard source imaging approaches. These include the detection of oscillatory coupling without rejection of zero-phase coupling, and detection of ongoing oscillations in deeper brain regions, where signal-to-noise conditions are unfavorable. We also show that iES provides a separate evaluation of oscillatory synchronization and desynchronization in experimental contrasts, which has important statistical advantages. The flexibility of iES allows it to be adjusted to many experimental questions in systems neuroscience.


Introduction
The role of neural oscillations in population codes of brain functions, and the possible mechanisms of inter-regional communication between brain regions are not entirely understood. Source imaging techniques with magnetoencephalography (MEG) or electroencephalography (EEG) are time-resolved, non-invasive tools used to test a great diversity of neurophysiological hypotheses [1]. In principle, MEG/EEG imaging can map multiple regional sources of oscillatory activity from external sensor data. However, spatial smearing and heterogeneous signal strength across brain locations limits the performance of current source imaging methods. Consequently, if nearby brain regions express an effect of interest, the area of stronger magnitude will mask the detection of weaker sources, as illustrated in Fig 1. The detection of multiple oscillatory sources therefore remains challenging to MEG/EEG imaging. This limits the insight about distributed brain dynamics that can be gained from the technique.
MEG/EEG localization of oscillatory generators typically relies on a procedure that is non optimal in terms of signal detection. Source time-series are first reconstructed using imaging or beamforming approaches [2]. Second, inferential statistics based on the experimental hypothesis are tested at each voxel of the source space-e.g., using the ratio of oscillatory power between two experimental conditions. Significant and spatially-distinct regional clusters are then interpreted as distinct sources of oscillations. This approach hinders the detection of weaker or deeper sources in the presence of stronger regional activity. We refer to this methodology as the standard approach.
We introduce a novel methodology to alleviate this problem. The technique performs imaging with embedded inferential and group prevalence statistics (iES) altogether. With iES, the experimental hypothesis is not deferred to the stage of statistical inference on the estimated source values. Rather, it explicitly constrains the solution to the hypothesis tested. In essence, iES reduces the (spatial) dimensions of the data, to detect and equalize the contribution of source components that are consistent with the tested hypothesis. The iES methodological apparatus is based on generalized eigen decompositions [see e.g. 3], nonparametric statistics [4,5] and subspace scanning [MUSIC, see 6]. The paper is organized as follows: Results: We first give a high-level description of the iES approach, starting with the basic principles and then illustrating a full group analysis with an MEG dataset. We then illustrate the advantages of iES using experimental data and simulations, and show that iES 1) has key statistical advantages, yielding improved detection sensitivity, 2) can be used in conjunction with the standard approach, for complementary estimation of source strengths, 3) improves the detection of functionally connected regions, and that 4) iES can implement a wide range of experimental hypotheses. Discussion: We then put iES in the context of previous work and discuss limitations. Materials and Methods: Finally, we provide all the experimental details and the full mathematical formulation of the iES approach.

Overview of the approach
We describe the basic principles of iES and illustrate the steps involved using a MEG data example. The method per se is detailed in Materials and Methods.
Basic principles. We propose to transcribe the experimental hypothesis into a quality function f(s) over a signal s. f(s) is defined such that it returns larger values if s is consistent with the hypothesis. In Fig 2a we show the quality functions featured as examples in this article, with multiple possible variations, as discussed below. For example, in the case where stimulusinduced responses in the gamma band  are expected, f(s) would be designed to return the ratio of gamma power between time segments when the stimulus is presented vs. when the stimulus is absent. Note that the tested design is a directed one: in the latter example, testing for gamma power increases vs. decreases are two different hypotheses that are evaluated separately. Fig 2a shows several use cases and signals s that are either consistent or inconsistent with the hypothesis quantified by f(s) values.
Let x[t] denote the MEG/EEG time-series recorded from an array of channels (Fig 2b) and X the data matrix (channels × time samples t). The quality function f is used in an optimization problem to identify spatial filters w j and spatially-filtered signals s j ½t ¼ w T j x½t in the data such that the quality function is maximized as argmax w f ðw T XÞ ð1Þ and spatial patterns p j , j = 1, . . ., D describing the spatial patterns that contribute to the recorded time-series in the decomposition as x½t ¼ p 1 Á s 1 ½t þ ::: þ p D Á s D ½t þ ð2Þ (Fig 2c). The combination of spatial filters, patterns and corresponding signals is comparable to the notion of 'components' in independent component analysis [ICA,7], which yields mixing (spatial patterns p j ) and unmixing (spatial filters w j ) matrices as well as ICA time-series (spatially filtered signals s j [t]). The spatial patterns of signals that conform to the tested hypothesis (i.e. have high quality function values f(s)) represent a subspace of the MEG/EEG channel space [the signal subspace, see 6].
The MEG/EEG forward model is a three-column matrix G(ρ) that describes the signal produced at the channel array by an elementary current source at location ρ with a dipole moment m, a vector that has three values corresponding to the three directions in space (Fig 2d). Therefore the forward fields of current sources with different dipole orientations also form a subspace of the MEG/EEG channel space. The iES method proceeds by scanning each elementary brain location of the source space. The source space can be a uniform 3-D grid of the brain volume, or restricted to the cortical surface (Fig 2e). At each tested brain location, the correspondence between the forward fields from this location and the data spatial patterns identified by the quality function is evaluated using the measure of subspace correlation (subcorr). This latter quantifies the smallest principle angle between two subspaces [6]. Intuitively, the data and the physical forward models are compared at each brain location, with respect to the experimental hypothesis. The procedure generates a map of possible sources, which activity accounts for the experimental effect of interest.

Fig 2. Basic principles of iES. a) Examples of designs:
The experimental design (shown as a black trace) determines the quality function f(s), so that this latter takes high values for signals consistent with the hypothesis (in orange; the signals that do not correspond to the tested hypothesis are shown in blue). b) MEG data: the multichannel MEG recordings are captured in the matrix X = {x[t = 1], . . ., x[t = T]}. c) Computing the signal subspace: spatial patterns P = {p 1 , . . ., p D } are extracted from the MEG data by optimizing the quality function with respect to spatial filters W = {w 1 , . . ., w D }. Whereas W is used to extract the signals of interest from the multichannel MEG data, P are the forward fields of these signals as they contribute to the measured MEG data. d) Computing the forward model: shown are the MEG spatial patterns G(ρ) generated by two tangential dipoles at location ρ in a single subject. e) Subspace correlation as a scanning metric: The spatial patterns from c) and d) span a subspace of the MEG sensor space. A grid of source locations is scanned with a subspace correlation metric [6], quantifying the smallest possible angle between the data and source subspaces. This yields a distributed map of scores, which highlights possible source locations consistent with the hypothesis. Since the effect strength is entirely captured in the quality function f(s), it does not directly influence the subcorr values. Therefore, the contribution of the hypothesis-consistent sources are equalized in the resulting maps, which contributes to their detection, regardless of their respective strengths.
Extension to group-level analysis. So far, subspace scanning techniques have been mostly used to identify equivalent dipole sources in single-subject MEG/EEG data [6,8]. We describe a principled approach to conduct group level analyses with the proposed method, using an example data set. The data was obtained from a variation of the visual attention experiment in [9], where a contracting circular grating was presented to participants in MEG. In each trial, after 3-5 seconds following stimulus onset, a change in the contraction speed occurs and participants had to indicate their perception of the change with a button press. To illustrate the methodology, iES was used to identify the regions where gamma-band (50-85Hz) power was stronger during visual stimulus presentation, with respect to baseline, prestimulus periods. Gamma oscillations generated in occipital visual regions are expected to be reliably enhanced by this paradigm [9]. The quality function f induced was defined as the ratio between the gamma power during the interval [1,3] seconds post-stimulus onset, and the gamma power during the baseline interval [-2, 0] seconds (0s corresponds to stimulus onset). Fig 3 shows the iES group analysis workflow in detail. First, the signal subspace estimation described above yields components from each participant that can be interpreted similarly to those of a principal component analysis (PCA) decomposition (Panel a). The subset of spatial patterns retained for source-space scanning corresponds to the iES components with highest quality function scores, exceeding a threshold f Ã induced . This threshold is determined by a permutation procedure under the null hypothesis of exchangeability of baseline and stimulus data segments. A permutation histogram of f induced values is obtained, and f Ã induced is set to the value that is higher than n perm − p crit /2 Ã n perm values of the permutation distribution, with p crit = .05 and the number of permutations n perm = 1000 in the example presented. In our example, this procedure yields five spatial components to be included in the definition of the signal subspace (Fig 3a). The sensor spatial patterns suggest occipital signal origins of stronger gamma-band activity during stimulus presentation. The iES decomposition, akin to PCA, produces orthogonal signal components s j . Thus the corresponding spatial patterns do not necessarily represent anatomically distinct sources: the spatial localization of corresponding sources is subsequently obtained via scanning of the source space.
The computation of this subspace is performed at the individual level, resulting in different numbers of subspace dimensions retained per subject. At the group level, our approach acknowledges that the effect being tested may be absent in some participants. Concretely, their data may not contain a spatial pattern whose f induced exceeds the critical value f Ã induced . Rather than pretending otherwise and averaging across all participants, as in the standard approach, we put forward the concept of population prevalence γ to account effectively for the variability of the tested effect in the group (Fig 3b) [see also 10, for similar discussions]. This notion enables to form a prediction on how many subjects in the sample are expected to show an effect. A prevalence null hypothesis, H 0 : γ γ 0 , can be tested using a simple binomial distribution. The null hypothesis can be e.g., that the effect is absent from the population (γ 0 = 0, global null hypothesis) or that it is present in less than half of the population (γ 0 .5, majority null hypothesis). The null hypothesis is rejected if observing the number of subjects presenting the effect has a probability lower than a critical value (here p crit = .05). In the present example of induced gamma oscillations, all subjects in the sample show the effect of interest. This means we can reject the majority null hypothesis (γ 0 .5), and the highest γ 0 that can be rejected at maximizing the f induced quality function), these signals present the largest ratio of gamma power between stimulus and baseline periods. Strong, tonic gamma oscillations are clearly visible after stimulus onset, along with reductions in alpha/beta power [9]. A 3.2-Hz oscillatory component is also found: it corresponds to the entrainment of lower-frequency neural components at the pattern-repetition frequency of the contracting circular grating.
We used the subcorr metric to produce source-level maps for each effect-prone subject, i.e. whose data features a non-null signal subspace. We show an example of an individual subcorr map in Fig 3e, which as expected, indicates a spatial peak in occipital visual regions. The Fisher-z transform arctanh(subcorr) can be applied to obtain a sharper map (referred to as subcorr-z). We then performed statistical inference at the group level, using group averaged subcorr maps in a permutation procedure (Fig 3f). A permutation distribution of the maximumstatistic is computed under the null hypothesis of exchangeability of signal subspace with a dimension-matched subspace drawn from the opposite end of the decomposition spectrum in Fig 3a. The null hypothesis thus reflects the assumption that the effects were not localized and spatially consistent across the tested cohort. This procedure yielded a statistically thresholded map of average subcorr values, highlighting the brain regions spatially consistent across the group, with an activity profile responding to the experimental question of interest. Here, the resulting map pointed to the visual cortex as the source of the gamma oscillations induced by the visual stimulus. This result was expected from published reports, and therefore further strengthens the validity of the proposed approach.

Distinct evaluation of positive and negative effects improves statistical power
The iES source maps highlight sources whose signals are consistent with a directed hypothesis across a group of subjects. When two experimental conditions are contrasted, this implies that two distinct source maps can be produced: for instance in the previous case example, one map corresponding to increased oscillatory power in one condition over the other; the other map corresponding to decreased oscillatory power. The benefit resulting from this is that mutual interference in the detection and statistical evaluation of the two sets of sources is avoided.
We demonstrate these methodological assets using the same experimental MEG data as above. We analyzed task-induced oscillations in the beta band (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30), with the hypothesis that they were strongly suppressed during attention-demanding tasks in the occipital visual cortex [9,11]. We also wished to test whether other brain regions would reveal a selective increase in beta power during stimulus presentation. This contrast thus serves to illustrate how a strong power effect (decreased beta power) can challenge the detection of weaker opposite responses (increased beta power) with the standard approach, but not with iES. Fig 4 shows results for the hypothesis of increased beta band power during stimulus presentation. The data from an example subject (Panel a) contained one spatial component consistent with that hypothesis. At the group level, only eight subjects out of 17 showed the effect of interest. Here, we shall emphasize the importance of the notion of effect prevalence, since the majority null hypothesis could not be rejected (Fig 4b). However, the prevalence null hypothesis can be rejected up to γ 0 = 0.22, which indicates there is a subgroup of the population from which our subjects were drawn, which shows the hypothesized effect. To better illustrate the significance of this notion, let us first assume the effect is not present in the population. With a probability of 0.95, one may still observe out of chance an effect in up to 3 out of the 17 subjects. The prevalence test therefore indicates that the observed data are unlikely under the assumption that prevalence is 22% or less (at a false positive rate of p < .05). Thus we pursued further the analysis of the subgroup of 8 participants (see Fig 4b), bearing in mind that the results may not generalize to the majority of the population. The validity of such a decision depends on whether the scientific question is pertinent to generic vs. restricted effects among participants. For instance, it can be particularly valuable for identifying effects that are more specific of a sub-type of participants in terms of behaviour or clinical condition. Fig 4c and 4d shows typical signal traces in a subject from the subgroup presenting stronger beta and alpha oscillations building up during stimulus presentation. The sharp waveforms and the combined alpha/beta spectral pattern were typical of the somatosensory mu rhythm [12,13]. The effect was localized to right postcentral regions, as shown in the example subject and the group subcorr maps (Fig 4e and 4f). This result replicates previous observations of lateralized beta oscillations during an attention-demanding task [14]. , only one significant spatial dimension was retained for the signal subspace contributing to stronger power in the beta band. b) Subspace computation, group level: γ = .22 was the highest population prevalence that could be rejected at a p = 0.05, thus the majority null hypothesis could not be rejected. The analysis was pursued with the subgroup (n = 8) of participants that showed the hypothesized effect. The purpose was to appreciate the spatial concordance across subjects and compare iES to standard source imaging approaches. c) Spatially filtered signals, example subject: induced power changes in the band of interest (beta, but also in alpha band) are clearly visible in 3 example trials. d) Spatially filtered signals, group level: induced power changes in the band of interest were found in the participant subgroup (n = 8). e) Subspace correlation maps, example subject: the hypothesized effect localized to the right post-central/ parietal cortex. f) Subspace correlation maps, subgroup level: the effect localized to the right post-central gyrus. Note that this effect cannot be generalized to the majority of the population that the subjects were drawn from (see b) but only to a subset, which may present interesting capacity for identifying subtypes in participants. https://doi.org/10.1371/journal.pcbi.1005990.g004

Imaging of neural oscillations
To compare these findings with those from the standard approach, we obtained source maps of log-power ratios using minimum-norm imaging kernels. We used the MNE implementation of Brainstorm, with default parameters [15, see also descriptions in Materials and Methods]. The resulting maps were statistically thresholded following the same permutation procedure based on the maximum statistic. Note that with this procedure, distinct maps of positive and negative effects cannot be produced. For comparison purposes, we used the data from the subgroup (n = 8) that showed the desired effect of higher beta power during stimulus presentation.
Fig 5a shows the complete beta band iES results (i.e. increases and decreases). In addition to the increased stimulus-induced beta power over right postcentral regions, we observed beta suppression in the visual cortex. The prevalence assessment revealed that this latter effect was observed in the entire group, and thus may generalize to the majority of the population. In the minimum-norm map (Panel b), the suppression of beta oscillations in the visual cortex was also readily observed, with similar spatial extension. However the increased, stimulus-induced beta oscillations over the right central regions were absent from the minimum-norm map (left) significant average subcorr map (p < .05, see text for procedure). Note that results were obtained from the subgroup of participants that presented the hypothesized effect (n = 8, see Fig 4). (right) histogram from observed data and permutation tests to derive a subcorr threshold corresponding to p < .05. b) Minimum-norm imaging results: average maps of log-transformed power ratios (stimulus/baseline, p < .05). Note that the distinction between positive and negative effects is not possible. The results were derived from the same subgroup (n = 8) to allow comparison with a), the results obtained with the full group (n = 17) are shown as an outline. Contrary to iES, no increase in beta power could be detected over the right post-central gyrus region, with the same subgroup of subjects. Unthresholded maps are shown in the supplementary material. (right) histograms of observed data and permutation tests to determine significance of minimum-norm maps at p < .05. Note how the strong negative effects inflated the permutation distribution and prevented the detection of the smaller positive effects. As shown using iES, positive and negative effects could be evaluated separately and specifically. produced from the 8 subjects presenting the effect in iES. The non-thresholded maps are shown in Supplementary Material, and confirm that a positive peak was indeed present in the minimum-norm maps, but was not deemed statistically significant. The reason for the observed discrepancy between methods can be understood from the permutation and data histograms ( Fig 5, right column): By definition, the permutation histograms of the log-power ratios were mirror images for the evaluation of positive and negative effects respectively. This was the case because we drew exhaustive permutations from the data from the 8 subjects (2 8 = 256). Thus every unique permutation of labels had a corresponding opposite permutation. The consequence was that the variance and spread of the resulting distribution were determined by the strongest effect in magnitude-here the negative effect of beta suppression. The histogram of the observed data indicated that the right tail of the histogram indeed did not reach the statistical threshold. The iES allowed to test two directed hypotheses separately. Hence the permutation distributions were distinct and adapted to each respective hypothesis, revealing the positive effect in the iES statistical source map that was absent in the standard approach.

Supplementary insight gained compared to standard approaches
We detail in Methods that iES requires the estimation of cross-spectral or covariance matrices, and their decomposition in the generalized eigenvalue framework. This means that, in addition to the subcorr statistical maps produced, a corresponding map of the standard approach can be obtained by applying a minimum-norm imaging kernel to those matrices, which allows plotting the value of the quality function f at each location of the source grid. Fig 5 shows an example of this approach to obtain a map of log-power ratios. We emphasize that the combined use of subcorr and minimum-norm source maps enabled by the proposed method provides complementary information with respect to the experimental hypothesis of interest.
We demonstrate such benefit using the same visual-attention MEG data, to detect the origins of narrow-band oscillations (Fig 2a). The corresponding quality function f narrow quantifies the ratio of signal power in a frequency range of interest with respect to the total power of the broadband signal. Such a quality function highlights signals with a peaky spectral profile [16], which is of specific interest when studying stimulus-independent ongoing oscillations. We used the data of the ongoing visual stimulus period ([1, 3] s after stimulus onset) to investigate the anatomical origins of three frequency bands of interest: theta (4-8 Hz), alpha (8)(9)(10)(11)(12)(13) and beta (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30). The reference broadband signal against which to contrast possible effects in the narrow frequency bands of interest was taken between 2 and 100 Hz.
We compared the subcorr statistical maps with the minimum-norm maps of f narrow (Fig 6). The log-transform of the ratios was not applied because negative effects were of no interest to the question, thus a symmetric measure was not required. A threshold 0 < f Ã narrow < 1 for selecting relevant signal subspace patterns was computed with the bootstrap procedure described in Materials and Methods. In the alpha and beta bands the results were similar between our approach and standard imaging. Commonly observed brain regions as strong sources of these ongoing rhythms were found [see e.g., 17]. Alpha activity was prominent in medial occipital-parietal regions; beta activity was stronger over bilateral sensory-motor regions. Alpha band oscillations were also found prominently over the right postcentral region, which parallels the finding of enhanced alpha and beta power during the stimulus period in the same brain area, as shown in the previous section.
We found differences between iES and minimum-norm maps in the theta band. The subcorr statistical map revealed involvement of the medial temporal lobes (MTL) bilaterally, and of medial frontal/anterior cingulate regions. Theta oscillations in MTL, including the hippocampus and parahippocampal regions, have been extensively described [18]. Due to their relative depth and therefore lower MEG signal-to-noise ratios (SNR), they have been considered more challenging to detect [19,20,21]. The MNE power-ratio maps though showed a lateralized distribution of theta activity in the right MTL. We argue that both results are not mutually exclusive: they indicate that both the left and right MTL were consistent sources of theta oscillations in the tested group. However, the effect strength in the right MTL was higher in the average power ratios of theta. Such insight could not be gained with either approach taken separately and required the direct comparison of the iES and MNE statistical maps. Fig 7 shows simulation results to illustrate and underline further the difference in sensitivity between the iES and standard approach. For each simulation run (300 iterations) we generated five minutes of data. Source time-series with a 1/f spectral profile were generated for 68 source locations distributed evenly across the brain according to the Desikan-Killiany atlas from Freesurfer [22]. For two of these locations (precentral left and right), we selectively amplified power in the frequency of interest (8)(9)(10)(11)(12)(13) to obtain a specified ratio f narrow between narrowband and broadband power (1-100 Hz). Whereas an f narrow of 0.6 was targeted for source 1, the targeted f narrow for source 2 was varied between 0.2 and 0.6. After generating MEG data from this simulation setup, we applied iES (with a f Ã narrow threshold of 0.22) and the standard approach to detect narrow-band oscillations in the frequency band of interest and computed a metric that quantified the probability of detecting both sources of narrow-band oscillations.  The sources of narrow-band signals were mapped for the theta, alpha and beta frequency bands using a) iES subspace scanning and b) power ratios from minimum norm imaging (narrow-band over broadband 2-100 Hz). iES allows for statistical thresholding across the group using permutation procedures that are equivalent for all use cases. The theta band results showed marked differences between the two approaches in deeper, medial temporal regions. iES revealed bilateral sources whereas MNE power ratio maps pointed at predominant source activity in the right hemisphere. https://doi.org/10.1371/journal.pcbi.1005990.g006 Imaging of neural oscillations approach scales with the differences in f narrow between the two sources, whereas iES' sensitivity is not influenced by uneven source activity and detects sources above the chosen threshold with a constant probability. This encourages using the different sensitivity profiles of the two methods in conjunction, to obtain complementary information as shown in the data example above. We show the results of a related simulation in S4 Fig where source 2 was a deep source (left parahippocampal) and source 1 a cortical source (left precentral). In this challenging scenario iES outperforms MNE in detecting both sources.

Assessment of functional connectivity
Because of spatial smearing, the study of functional connectivity is a challenging problem for MEG and EEG source imaging (see Fig 1). Since the seed region is maximally correlated with itself and neighbouring regions, with correlated time series due to field spread of the MEG/ EEG inverse operator, functional connectivity maps tend to be biased towards artificially inflated values of connectivity measures. This issue is discussed in [23] and generally addressed Simulation results comparing sensitivity of iES and standard approach. a) Examples of simulated time-series that follow a 1/f spectral distribution (grey trace) or target a pre-specified f narrow , which is the ratio between narrowband and broadband power (blue traces). b) Simulation setup: Two sources of interest in blue targeting pre-specified f narrow (blue traces) are embedded in background brain noise composed of 1/f signals evenly distributed across 66 locations. c) Metric of detection probability: We quantified the probability that the two sources of interest were detected in a source map by using a range of different thresholds: the two sources were detected, if they were contained in two separate clusters after thresholding. Here we show 4 different thresholds in two simulation scenarios using the standard imaging approach. In the first scenario, sources were detected with 2 out of 4 (detection probability: 0.5) threshold values. In the second scenario, sources were detected only with 1 out 4 (detection probability: 0.25) threshold values. This configuration illustrates the issue of concurrent sources with different strengths on the detection of separate clusters of activity. d) Comparison of methods: the maps from each simulation run were thresholded using 50 different values to estimate a detection probability as in c). Since the range of data values for both MNE and iES were different, we normalized the detection probability by the maximum value obtained in each method. Thus we did not compare the absolute detection probability between the two methods, but rather how it varied with respect to the difference in f narrow , between the sources of interest. https://doi.org/10.1371/journal.pcbi.1005990.g007 with methods that discard all contributions of zero phase-lag time series, either by orthogonalizing signals [23] or via measures of the imaginary part of coherence [24]. However, zero-lag coherence between distant regions is plausible theoretically [25] and was observed physiologically [26]. We demonstrate the relevance of iES to address this issue, by studying amplitude correlations in the alpha band (8-13 Hz) with respect to an anatomical seed placed in the sensorimotor cortex. The tested hypothesis was to reveal amplitude correlations with homologous contralateral brain regions [23,27,28]. Fig 8 shows results from resting-state data obtained during the same recording session as the visual stimulus experiment. a) shows example time series of co-occurring oscillatory bursts, which form the basis of amplitude correlations between two distant brain regions. The time series were extracted from bilateral central regions. We show occurrences of 270˚/90˚phase differences during alpha bursts-note that the phase estimation of MEG source signals has a 180˚ambiguity due to arbitrary conventions on source direction [29]-and of 180˚/0˚phase differences, which would be discarded by other methods [23,24]. We argue that the zero-lag correlations shown here are not spurious, as evidenced by their differences in waveform and amplitude dynamics. This data example provides a proof of principle that studying zero-lag connectivity using MEG is achievable. We next proceeded to map significant inter-regional amplitude correlations in the presence of field spread.
We extracted the source time series y ref from the left central sulcus location that was the closest to the activation peak (MNI coordinates [-39, -27, 55] mm) corresponding to the search term 'finger' in the Neurosynth meta-analysis tool [30]. We defined the iES quality function Correlation of alpha amplitudes between the seed region (circle) and the rest of the brain, using minimum-norm imaging in an example subject. c) (left) iES subcorr map showing source locations whose amplitudes correlated with the seed region's at r > .4. The homologous contralateral region is emphasized in this map. (right) the same map with the data projected away from the spatial pattern of the seed region. d) Same as b) but averaged over the group e) same as c) but averaged over the group and statistically thresholded using the permutation approach explained above. https://doi.org/10.1371/journal.pcbi.1005990.g008 Imaging of neural oscillations f ampcorr (s, y ref ) as the correlation between the amplitude of other source time series s and y ref .
The outcome to this optimization process is the third use case of iES and is illustrated in Fig  2a. The optimization is done using a solution described in [31] and involved splitting the data in epochs of one second length. We set a correlation threshold of r > 0.4 for spatial components to be included in the signal subspace for the subcorr analysis. Fig 8b and 8c shows-in a single subject example-that iES was able to reveal the contralateral anatomically-homologous region as the primary distant connected region with the reference brain location. The conventional minimum-norm based map of correlation values was dominated by spurious crosstalk correlation surrounding the seed region. The performance of iES is explained by the equalized contribution of spatial components that are consistent with the embedded hypothesis (r > 0.4). To further limit the contribution of the seed region to the data, it is possible to project the signal subspace and forward fields away from the spatial forward field of the seed region, as illustrated in Fig 8c. The group analysis further reveals that connectivity maps were dominated by crosstalk effects from the seed reference signal, both in the minimum-norm based maps and in the raw subcorr map (Fig 8d and 8e). Projecting away the seed's contribution before computing every subject's maps was necessary to confirm the hypothesized contralateral coupling. Note that with iES and in contrast with other approaches [23], the temporal dynamics of the seed region are not projected away from the data; only the spatial topography of the seed region is subtracted from the sensor data. Thus iES does not exclude the detection of physiological zero-lag coupling.
We provide a simulated example in S2 and S3 Figs that illustrates this point. If there is a source that is correlated in amplitude with the seed and the underlying oscillations have a zero phase difference, the two sources will be captured by only one subspace pattern. In this case the coupled region will be picked up by the iES procedure only after the signal subspace and forward fields are projected away from the contribution of the seed topography.

Applicability to a wide range of experimental questions
As summarized in Fig 2a, iES can be used for a greater variety of experimental designs: whenever a reference signal y ref defined 1) on a trial-by-trial basis or 2) as a continuous signal is considered, f ampcorr is used to obtain subcorr maps of sources, whose source dynamics correlate with y ref . We illustrate such case in Fig 9, using simultaneous MEG and pupil diameter recordings. We first formed the iES hypothesis based on recent demonstrations in mice [32] that continuous pupil diameter fluctuations correlated with alpha power at rest. Measures of pupil diameter were extracted from continuous video eye-tracking recordings by fitting an ellipse to the pupil on a frame-by-frame basis. We demonstrate the case of trial-by-trial correlations by analyzing pupil diameter changes prior to visual stimulus presentation from data presented in Fig 3. The signal subspace was defined with spatial components who were deemed significant below a p-value of 0.05 computed by a shuffling procedure across trials. The iES maps indicated brain regions in the occipital cortex, consistent with the upcoming onset of a visual stimulus.
We emphasize that a specific strength of the iES approach is its versatility: it can be extended to a great variety of experimental designs and research hypotheses, since the experimental question is formulated as an optimization problem. We derive in Methods the mathematical formulation for iES coherence with a reference signal, as an additional experimental use case. The experimental hypotheses discussed here all have corresponding quality functions that can be solved analytically. An identical framework can be used for hypotheses that require numerical optimization of the corresponding spatial filters. We foresee that the introduction of the iES approach will establish a generic framework for an increasing number of experimental contexts related to a growing diversity of research questions.

Discussion
We developed the imaging with embedded statistics (iES) method to produce image maps of the sources of neural oscillations from MEG/EEG sensor data. In this article, we showed with ground-truth simulations and experimental data, that iES identifies source patterns that are challenging to standard approaches, especially when masked by field-spread from other sources in the volume. Specifically: 1) iES generates separate maps for event-related increases and decreases in oscillatory power, which facilitates statistical inference; 2) iES can be used in pair with standard approaches such as wMNE, with iES identifying source regions and wMNE extracting their respective amplitudes. We also showed that 3) detection of functionally connected sources presents an extreme case of source strength imbalance that can be solved using iES, without rejecting the possibility of zero-phase delay coupling between regions. And finally: 4) iES can be flexibly extended to a greater variety of experimental designs, by defining quality functions suitable with the neuroscience hypothesis to be tested with the collected data.
The iES method builds on previous methodological work. In particular, we acknowledge inspiration from previous subspace scanning approaches adapted to MEG/EEG, such as the multiple signal classification (MUSIC) [6] and RAP-MUSIC [8] methods. These previous approaches produced subspace correlation maps, which were subsequently pruned to sets of discrete equivalent dipoles. These methods have inspired an abundant literature, with multiple extensions covering different experimental questions [16,31] and the estimation of time series interdependencies in functional connectivity [33,34,35,36]. We wish to emphasize however that MUSIC-type methods do not provide a clear path to determine the dimensionality of the signal subspace for each subject. With iES, we propose novel solutions using nonparametric statistics for each use case described. Previous suggestions for component selection [37,38] are applicable mostly for the original MUSIC solution, which performed principal component analysis of the event-related fields. It is also unclear how previous methods could indicate how many elementary dipole sources could be adjusted to the data, e.g., how many iterations of RAP-MUSIC were required. With iES, subspace correlation maps are used to perform statistical inference at the group level, thus circumventing the necessity of registering heterogeneous discrete elementary dipole models across participants. The goal of iES thus is to produce a distributed and statistically thresholded map. Depending on the data, experimental context and hypothesis embedded in iES, such maps can reveal one spatial peak-e.g., in the visual cortex in the example illustrated in Fig 3 or more source regions, as for example in Fig 6, which highlights sources of beta-band oscillatory activity over bilateral sensorimotor regions. The iES methodology also explicitly addresses the recurrent issue of heterogeneity in the expression of an effect of interest across individuals in a tested group of participants. We propose to use an innovative approach in the field, based on prevalence statistics. We acknowledge that prevalence testing was used in multivariate decoding analyses [10], but the strategy used in iES is quite distinct.
In sum, iES has broad, practical value to research using M/EEG imaging, where group-level inferences are very common. In addition, iES features single-subject analysis and group-based evaluation through the concept of effect prevalence, which can lead to new insights on subject stratification. We consider this feature has great potential value in identifying subgroups of participants e.g., responding differently to a specific paradigm manipulation in experimental psychology and cognitive science, or a therapeutic intervention in clinical trials. A further interesting property of iES is that we separate two steps of the analysis, namely 1) evaluating the presence of a hypothesized effect in single subjects (as well as its prevalence across the group), e.g. increased oscillatory power in a frequency band of interest during stimulus presentation; and 2) the spatial consistency of this effect across the group/subgroup. Standard MNE imaging would not be able to detect an effect that is present in all the subjects of a given study, but is not spatially consistent across the group. In contrast, with iES we would detect that the effect is present in all subjects, even if the group spatial maps do not show peaks that exceed the statistical threshold. Again, we see this as a strong asset of iES, since it could prompt a more in-depth analysis of spatial patterns in individuals. This, in turn, could lead to the discovery of subgroups that differ in the spatial patterns they produce in a given experimental paradigm.
Limitations of iES in its present form are essentially in the definition of hyper-parameters in an ad-hoc manner. This is not specific of iES, as even the standard approaches do not provide, in practice, hyper-parameter estimation from the data or based on strong theoretical background. For instance, the estimated covariance matrices in the first step of iES are regularized using a fixed regularization parameter (see Fig 10 for the two variants used in this paper). Recently, automatic methods [39] for selecting those regularization parameters have been proposed, via maximization of the likelihood of unseen data (under the cross-validation principle). We acknowledge that these methods can be easily used for hyper-parameter selection within the iES framework presented. Another iES hyper-parameter is the p-value threshold used when selecting the spatial patterns to define the signal subspace. We used one of the common thresholds (e.g., 0.001 in the bootstrap procedure for the results in Fig 6), and this certainly lacks theoretical foundations. To improve the rationale on this selection, a more established theoretical framework on the physiological mechanisms explaining the relative power ratios observed between the typical frequency bands of electrophysiology would be required. Some recent work is going in that direction and could inform the selection of iES hyper-parameters, for example by modelling power spectra as a mixture of 1/f spectral noise, as well as narrowband oscillations [see e.g., the discussion in 40,41]. This open scientific question is beyond the scope of the methodological advances we propose, and models embedded in iES can be improved as advances are made.
There is a fast growing literature on supervised learning methods for M/EEG that test specific questions regarding encoding or decoding of e.g., stimulus features [42,43]. For the most part, these latter are framed as regression or classification problems and, as purely datadriven methods, they don't have a model of how a specific spatial pattern was obtained. Our contribution with iES is strongly related to these aspects: iES reframes an experimental question into an optimization problem, which can be seen as a supervised-learning objective to identify spatial patterns of interest from the data. We also propose with iES to interpret these patterns using the physical forward model of MEG/EEG, which yields group-level maps in source space. In principle, iES can be used in studies using supervised learning to find physiological sources of a discriminative pattern, provided the learning procedure is aimed at optimizing linear channel weights. This would not be possible if non-linear methods are used, e.g. with support vector machines with radial basis functions. Finally, we described iES with use cases that have quality functions yielding closed-form solutions. Future work should reveal how iES can be used with non-convex quality functions (i.e. presenting multiple local minima) solved with the typical apparatus of numerical optimization.

MEG data
Participants and ethics statement. 17 healthy participants were recruited (21-45 years; 5 female). The study was approved by the Montreal Neurological Institute's ethics committee (NEU-11-036), in accordance with the Declaration of Helsinki. All participants gave written informed consent and were compensated for their participation.
Stimuli. Subjects were presented with a variation of the visual stimulation paradigm in [9]: A circular sine wave grating (diameter of 5 with 100% contrast) contracts towards the fixation point (velocity: 1.6 deg/s). The contraction accelerated (velocity step to 2.2 deg/s) at an unpredictable moment between 3-5 seconds after stimulus onset. Subjects had to indicate with an index-finger button press that they detected the velocity change. The button press ended one trial and the stimulus was turned off. Inter-trial intervals were 5 seconds long with a jitter.

Imaging of neural oscillations
During the inter-trial interval subjects were presented with a central fixation cross. Stimuli were generated using the Psychophysics Toolbox [44].
Experimental procedure. Participants received both oral and written instructions on the experimental procedure and the task. The recording session started with a 5-minute resting state run with eyes open. The participants were presented with 10 test trials, to become familiar with the task. They were then presented with a total of 240 stimulus sequences (trials). Participants performed 60 trials per acquisition block. After each block, they received a feedback on the accuracy of their responses. A trial was considered correct if subjects responded within 500 ms after the actual velocity change occurred. Between each block, the participants were given a break of self-determined length. After completion of the 240 trials, subjects were given a 15-minute break. A further 5-minute resting-state recording concluded the session. Data acquisition. The participants were measured in a seated position using a 275-channel VSM/CTF MEG system with a sampling rate of 2400 Hz (no high-pass filter, 660 Hz antialiasing online low-pass filter). Three head positioning coils were attached to fiducial anatomical locations (nasion, left/right pre-auricular points) to track head movements during recordings. Head shape and the locations of head position coils were digitized (Polhemus Isotrak, Polhemus Inc., VT, USA) prior to MEG data collection, for co-registration of MEG channel locations with anatomical T1-weighted MRI. Eye movements and blinks were recorded using 2 bipolar electro-oculographic (EOG) channels. EOG leads were placed above and below one eye (vertical channel); the second channel was placed laterally to the two eyes (horizontal channel). Heart activity was recorded with one channel (ECG), with electrical reference at the opposite clavicle.
A T1-weighted MRI of the brain (1.5 T, 240 x 240 mm field of view, 1 mm isotropic, sagittal orientation) was obtained from each participant, either at least one month before the MEG session or after the session. For subsequent source analyses, the nasion and the left and right preauricular points were first marked manually in each participant's MRI volume. These were used as an initial starting point for registration of the MEG activity to the structural T1 image. An iterative closest point rigid-body registration method implemented in Brainstorm [15] improved the anatomical alignment using the additional scalp points. The registration was visually verified and adjusted manually, if necessary.

MEG data preprocessing
All MEG data analysis steps were performed with Brainstorm [15], with the novel approaches described in this paper implemented as a Brainstorm plug-in written in MATLAB (available through: https://github.com/pwdonh/ies_toolbox).
Artifact removal and rejection. Eye-blink and heart-beat artifacts were removed from MEG data using a PCA-based signal source projection (SSP) method, using recommended procedures [45]. The ECG and EOG channels were used to automatically detect artifact events. Noisy MEG channels were identified by visually inspecting their power spectrum and removing those who showed excessive power across a broad band of frequencies. The raw data were further visually inspected to detect time segments with excessive noise e.g., from jaw clenching or eye saccades. When epoching the data, we automatically excluded all trials that overlapped with these noisy time segments. Sinusoid removal at the power line frequency and harmonics (60, 120, 180 Hz) was applied to the continuous data. A high-pass filter above 1 Hz was also applied to reduce slow sensor drifts. The MEG data were centered around the baseline mean after epoching. All the filters used in the current study are zero phase shift non-causal finite impulse filters coded and documented in Brainstorm.
MEG data were epoched to the interval [-2, 3] seconds around the visual stimulus onset. We refer here to the stimulus period as the interval [1,3] seconds, and the baseline period as [-2, 0] seconds with respect to visual stimulus onset. The intervals were chosen to focus on the steady-state part of the oscillatory response, as opposed to the transient at stimulus onset [9].
Intra-subject coregistration. Prior to the computation of signal subspaces, we performed a between-run coregistration of the MEG data based on recorded head positions, using the movement correction method similar to [46] available in Brainstorm. Briefly, we computed forward models G k based on the head positions of different runs k and G avg using the average head position. Then we computed coregistration operators O k to project the MEG data from different runs into the same space as O k X k . O k was computed as O k ¼ G avg V n S À 1 n U T n , using the singular value decomposition (SVD) G k = USV T truncated corresponding to the largest n singular values. The index n is set so as to preserve 99.99% of the squared singular value spectrum. We additionally took into account that different runs had slightly different SSP projectors applied (see above). We thus apply these projectors to the forward fields of individual runs before computing the coregistration operators.
Source models. We defined a volumetric source grid on the MNI152 2009c nonlinear anatomical template [47], using an adaptive procedure: an outer layer of 4000 grid points was produced based on a brain envelope covering cortical and subcortical structures. This outer layer was then shrunk and downsampled by a factor of 2.2. This procedure was repeated to result in a total number of 20 layers containing 25,740 grid points. For each subject, we computed a linear transform of individual anatomy to MNI coordinates using affine coregistration as implemented in SPM12 [function spm_maff , 48]. We applied the inverse linear transform to project the default source grid onto each subjects' individual anatomy. Thus a source grid consisting of 25,740 points was produced for each subject, with one-to-one correspondence between points across subjects, as produced by Brainstorm volume source grids.
Forward modeling of neural magnetic fields was performed using the overlapping-sphere model [49]. Conventional MEG source imaging was obtained by linearly applying the weighted-minimum norm operator [2]. The wMNE was obtained as follows (using default hyperparameters in Brainstorm): First, a regularized version of the noise covariance matrix, estimated from same-day empty-room MEG recordings, was produced by adding a diagonal matrix to the original noise covariance array. The additive diagonal elements were channelspecific regularization weights set to 10% of the average of the noise variance observed across MEG channels from the empty-room recording. Second, the regularized noise covariance array was subsequently whitened, following its eigendecomposition. Third, the lead-field matrix obtained from forward-field modelling of the elementary current dipoles in the distributed source model was whitened, by applying left multiplication with the whitened and regularized noise-covariance estimate. Fourth, compensation for uneven source depth was applied by computing source-specific depth-compensation weights as the minimum between the inverse of the norm of each source's forward field and 10 times the smallest source forward field norm in the model (for regularization purposes). Fifth, the source covariance matrix was assembled from the cross-product between all depth-weighted source forward fields. Sixth, the wMNE linear kernel was obtained via a regularized version of the source covariance matrix, with a hyperparameter set to 3. Note that the code of this implementation of wMNE is open-access (bst_wmne.m) in the distribution of Brainstorm and is consistent with that of the MNE software.
For the extraction of a seed time-series for functional connectivity analysis in Fig 8 we found the optimal dipole direction at the seed location using SVD of the filtered 3 component time-series extracted with the MNE kernel K(ρ).

iES formulation
The iES method described in this paper is based on subspace scanning, which processes the entire spatio-temporal MEG data matrix X, instead of reconstructing neural activity independently at each time point for the whole source space [see e.g. MUSIC, 6]. The method features two steps, as shown in Fig 2: 1) extraction of the relevant spatial patterns from the data (signal subspace identification), and 2) scanning of the source space for contributions that explain the identified spatial patterns (subspace scanning step per se).
Subspace scanning. Here we describe how subspace scanning was used in the original MUSIC approach [6]. First a set of MEG topographies is identified that captures the signal components of the MEG data matrix. We define the notion of signal in the following section. In [6], this first step was equivalent to performing a PCA of the event-related average MEG signals. The D components corresponding to the largest PCA eigenvalues spanned the signal subspace: span(P s ). P s is a full column-rank M × D matrix where M is the number of sensors and D is the dimensionality of the signal subspace. A particular MEG topography v lies within the signal subspace, if there exists a linear combination t of the columns of P s such that v = t T P s . A geometric measure quantifies how close a particular MEG topography lies to the signal subspace. For instance, the cosine of the angle between v and the projection of v onto P s is a suitable measure [6].
A dipolar source at location ρ is described by an orientation (θ) and an amplitude (a) parameter. All the MEG topographies that can be generated by this source are described by a linear combination of the forward fields of dipoles along the three orthogonal spatial directions [2]. The resulting 3-column forward field matrix G(ρ) thus also spans a subspace. A MEG topography produced by this dipole lies within the signal subspace if there exist linear combinations u and t such that u T G(ρ) = t T P s . Due to noise in measurements and inevitable approximations in the forward model metric, a perfect match cannot be expected. Thus we use the subspace correlation metric as the cosine of the smallest principal angle between subspaces subcorrðGðrÞ; P s Þ ð3Þ as defined in [6]. This metric quantifies how close the two subspaces lie to each other, and thus how well a dipole source at the scanned location fits the signal subspace. This metric is applied at each possible location across the anatomical volume.
Computing the signal subspace. We now describe how finding the signal subspace can be seen as the solution to an optimization problem, which opens to a wide range of new possible applications. In the standard MUSIC case, the first column of the signal subspace p 1 is a vector/topography that, when applied to the event-related average of the MEG data X, results in a signal that has maximum variance (broadband power): it is a solution to the optimization problem subject to a norm constraint on p. The next subspace column p 2 is the solution to the optimization problem where P ? is the orthogonal projector away from the first subspace column P ? ¼ I À p 1 p T 1 . This corresponds essentially to a PCA of the event-related average X. More generally, the signal subspace could be constructed by the solution of an optimization problem where the function f is chosen according to the experimental question of interest. The timeseries s is the signal obtained by applying the spatial filter to the data as s = w T X. In standard MUSIC, the experimental question of interest is to find the sources that have the strongest contribution to the event-related responses, thus the quality function is where " s is the event-related trial average of s. However, many other experimental questions can be expressed as an optimization problem. For example we might be interested in finding sources whose power is correlated with a reference signal y (such as an EMG recording or an audio stimulus envelope). In that case we would set the quality function as where Corr(a, b) is the correlation of signals a and b. While, in principle, it is possible to use any quality function and proceed with numerical optimization, the subspace method is specifically attractive for quality functions that can be solved analytically for computationally efficient implementations. Here we focus on a set of quality functions that can be solved using the generalized eigenvalue problem (GEP). We show solutions for four different experimental use cases, three of which are illustrated in Fig 2a. Subspace computation using the GEP.
The generalized eigenvalue problem [GEP, see e.g. 50] Aw = λBw, for symmetric matrices A and B, arises in optimization situations like or equivalently To show how the GEP can be used to define a subspace, we focus on a) induced responses as a first use case (see Fig 2a). Here we are interested in finding sources whose power in a frequency band of interest [f 1 , f 2 ] differs between two conditions or time periods, e.g. stimulus and baseline periods. Thus the quality function becomes where Pow(s a,b )[f] is the power of a signal s at frequency f in time periods a and b. The power of a signal in a given frequency band [f 1 , f 2 ] can be approximated by the variance of the signal filtered in that frequency band. The quality function can thus also be written as where superscript filt indicates that the signal was filtered in the frequency band of interest. For readability, we will drop this superscript in the following. Since the bandpass-filtered signal is zero mean, we can compute the variance using the dot product 1 LÀ 1 s T s where L is the number of time samples. The quality function thus becomes where C is the estimated covariance matrix of the filtered MEG signals X and w is a spatial filter topography. This is now in the form of the GEP shown above and has been used in the field of brain-computer interfaces as Common Spatial Patterns [CSP, e.g. 3,51]. Alternatively, one can define the quality function directly in the frequency domain, and compute C (a,b) as the average of the real part of the estimated cross-spectral density matrices in the frequency band of interest where C XX [f] is the estimated M × M MEG cross-spectral density matrix at frequency f. The GEP can now be solved by defining a whitening projector from the SVD: USV T = C b , which equalizes the variance along the principal axes of C b , as required in the constraint of Eq 10. We then solve the ordinary eigenvalue problem where the eigenvector ϕ is now a spatial filter in the whitened data space. The eigenvalue λ provides the ratio of power in the two conditions, thus is equal to f induced . This means that the signals of interest, which maximize the quality function in Eq 12, can be estimated from the MEG data asŝ where w combines the two steps of whitening (P ? ) and filtering in whitened space (ϕ T ) to obtain a spatial filter in the data space as in Eq 13. The data generated by a specific source signal can in turn be estimated byX where p is the spatial pattern vector, or forward field, of the source signal in sensor space, since an inverse whitening step (US 1/2 ) is applied to the forward pattern in whitened space (ϕ). [see 52, for further discussion on the distinction between spatial patterns and filters]. Solving the GEP this way, one obtains M spatial patterns p j that can be ordered according to their corresponding quality function scores f induced = λ. The columns of the signal subspace matrix P s are then defined by the spatial patterns that exceed a threshold f Ã induced , yielding an M by D subspace matrix, where D is the number of spatial patterns exceeding the threshold. We discuss the estimation of this threshold below. As described above, the anatomical source space can then be scanned by computing subspace correlations with the forward fields at each source location, using where we obtain the left singular vectors of the two matrices such that GðrÞ ¼ U G S G V T G and P s ¼ U P S P V T P and the squared subspace correlation corresponds to the maximum eigenvalue obtained as Additional use cases based on the GEP. We describe three other use cases that can be solved using an appropriate quality function in combination with the GEP (see Fig 2a). The solutions are convenient in that they only require to change the definition of C a and C b in Eq 13.
We define a signal showing b) narrowband oscillations as a signal that has increased relative power in a frequency band of interest with respect to broadband power. We thus define a signal frequency band of interest ½f s 1 ; f s 2 and a broad noise frequency band ½f n 1 ; f n 2 . The quality function then becomes Analogously to Eq 13, this quality function can be expressed in the form of the GEP as where Solving the GEP, we obtain spatially filtered signals that are ordered according to their ratios of power f narrow in the signal and noise frequency bands. This approach is similar to what has been described in [16] as spatio-spectral decomposition.
As next use case, we consider the case of c) amplitude modulation using a solution described in [31]. Here we wish to find sources whose amplitude fluctuations in a frequency band of interest covary with the value of a reference variable. This might be a slow time-varying signal y ref , or a variable that is defined on a trial-by-trial basis such as reaction time or task difficulty. Here we describe the former case, but the latter follows easily [31].
The data are split into epochs denoted by the index e. Epoch length needs to be short enough to allow capturing fluctuations in the reference signal y ref , and long enough to estimate the power of data signals X filtered in the frequency band of interest. The results in where X(e) denotes the data matrix in epoch e. We thus maximize the covariance between the power of s = w T X(e) and the value of y ref ðeÞ normalized by the power of s = w T X = 1. Assuming that y ref is a zero-mean and unit-variance signal, this can be solved in the GEP framework by setting where C X (e) is the estimated covariance matrix of the filtered MEG signals in epoch e. C a thus represents a weighted (by y ref ðeÞ) average, and C b the unweighted average of the single epoch covariance matrices. Please refer to [31] for a derivation of these results. We can obtain the correlation values from the above as In the analysis examples we used f ampmod to compute the spatial filter basis using the GEP in a computationally efficient manner. Ordering and selecting the components to be included in the signal subspace was then based on f ampcorr . When the research hypothesis requires testing for source dynamics that are d) coherent with a reference signal y ref at a specific frequency, the quality function becomes The reference signal y ref can be an external stimulus such as the envelope of an audio signal, a simultaneously measured peripheral signal such as EMG, or a neural time-series extracted using source imaging. Magnitude squared coherence is computed as the ratio of cross-spectral to auto-spectral densities as where C sy [f] is the estimated cross-spectral density between signals s and y at frequency f, and C yy [f] is the auto-spectral density of signal y. Since C yy [f] is constant, we can leave it out of the quality function and remain with Now setting s = w T X we get where C Xy [f] is the column vector containing the estimated cross-spectral densities between the reference signal and the MEG signals, and the H superscript stands for conjugate transpose. The optimization problem can be solved by invoking the GEP as in Eq 22 and setting Note that, as the matrix C Xy ½f C H Xy ½f is Hermitian, w T C Xy ½f C H Xy ½f w will be a real number and thus be equal to w T ReðC Xy ½f C H Xy ½f Þw. Hence we only need to keep the real part for C a . Covariance regularization. As described above, at the core of our approach lies the computation of a subspace using the GEP. This requires the inversion of matrix C b , which can be numerically unstable if some of its singular values are small (see Eq 15 for its influence on the whitening projector). We thus regularize matrices C a and C b by adding values to the diagonal as where α is a regularization parameter. We refer to this technique as diagonal loading. where s (a,b)j are the singular values of the narrow-band covariance matrices C (a,b) computed from stimulus (a) and baseline (b) periods, respectively. This can be understood as the expected f based on the overall power across sensors in both time periods. Because these components do not carry physiological information, we make sure to never include them in the signal subspace. They were detected using a simple bootstrap procedure: a confidence interval on the mean power of each component (s aj and s bj ) is estimated from the unregularized data by sampling with replacement from epochs e = 1, . . ., E. We then verify whether s aj and s bj computed from the regularized data lie within the 99.9% confidence interval, and discard the components where this is not the case. We used this approach in the analysis examples on induced responses and narrow-band oscillations (Figs 3, 4, 5, 6 and 7) and set α = 0.05.
We also describe an alternative regularization approach, referred to as truncated SVD. This entails removing the columns associated with the smallest singular values from U during the computation of the whitening projector in Eq 15. We define the regularization parameter and keep the singular values making up 100(1 − )% of the cumulated singular value spectrum. We see in Fig 10b that this results in a smaller number of components extracted from the GEP, however yielding similar f values at both ends of the spectrum. We used this approach in the analysis examples on amplitude modulation (Figs 8 and 9) and set = 0.001.
Estimating the dimensionality of the signal subspace.
All the methods described above result in a matrix P where each column p j is a spatial pattern associated with a quality function score f j . We now need to determine which of these spatial patterns to include in the signal subspace P s that will be used for scanning. This can be done by setting a threshold on the quality function scores f j in a hypothesis-driven way. These scores are readily interpretable as e.g., the power ratio between conditions (induced responses) or the correlation between a reference signal and neural amplitude time-series (amplitude modulation). In the following we show how we can set a threshold based on approaches from non-parametric statistical testing, including permutation and bootstrap procedures. The result is a M × D signal subspace matrix P s , with dimensionality D being the number of significant spatial patterns. In the case of permutation tests, spatial patterns p j are included in the signal subspace matrix P s if the associated quality function score f j is inconsistent with the null hypothesis. In the case of the bootstrap, spatial patterns are included in P s if the associated bootstrap distribution of f j is consistent with the alternative hypothesis. We describe these procedures in the following.
For a) induced responses we compare oscillatory power between two conditions using permutation testing. We describe the case of stimulus-baseline contrast where, for each epoch e, a data matrix for baseline and stimulus periods is available to perform a paired test. Other cases can be derived easily using standard approaches in non-parametric statistics [see e.g. 4,5]. Under the null hypothesis of no difference, the condition labels are exchangeable with respect to the statistic of interest f (here f induced ), which is the ratio of power between the two conditions which is defined for each of the potential columns j of the signal subspace matrix P s . The data are divided in e = 1, . . ., E epochs, from which we compute empirical covariances C a (e) and C b (e). We run O permutations, where at each iteration, a binary permutation vector ω of length E is drawn. At each permutation we solve the GEP based on the permuted condition labels and compute a maximum statistic as and the f max values are logged at each iteration. We then obtain a null distribution of f max (assuming exchangeability of the condition labels) against which to test the observed f j 's to obtain a permutation p-value. In this paper we use O = 1000 permutations.
In the second use case b) narrowband oscillations the power ratios f j (f narrowband ) will differ depending on the frequency band of interest. Due to 1/f in electrophysiology power spectra, low-frequency bands have higher f j 's than high-frequency bands. To find spatial patterns p j whose relative power stands out from the rest of the activity, we define an expectedf , as the ratio of overall power in the narrow-and broad frequency bandŝ where s (a,b)j are the singular values of the estimated narrow-and broadband cross-spectral densities C (a,b) as defined in Eq 23. We use a bootstrap procedure to find the dimensions j that reliably lie above this expected power ratio. A large number of bootstrap samples can be obtained by sampling with replacement from the epoched data, and logging the mean values over the f j (e) of each selected epoch. A confidence interval based on the obtained bootstrap distribution is obtained and a dimensionality D up to which the confidence interval does not contain the expected power ratio is therefore defined. In Fig 6, we used a confidence interval of 99.9% to define the signal subspace. Use cases c) amplitude correlation and d) coherence compute a measure of temporal association between time-series. The temporal ordering between the reference signal y and the data X can be scrambled under the null hypothesis of no association. Since the computation of f ampmod involves splitting the data into epochs e (see Eq 24), we can compute a null distribution of the respective f values by randomly re-assigning epochs between the reference signal y(e) and the MEG data X(e). At each iteration the GEP is solved resulting in a single null distribution of f Ã against which to test all the observed f j .
Projection of the seed topography for functional connectivity analyses.
We can project out the topographic contribution of the seed location ρ s in a functional connectivity analysis as shown in Fig 8. Using subspace correlation we define the topography g at location ρ s with orientation θ that maximizes the fit with the signal subspace as g ¼ Gðr s ; yÞ; where y ¼ argmax y subcorrðGðr s ; yÞ; P s Þ ð38Þ Then we find the orthogonal projector to be applied to both the signal subspace and the leadfield matrices, so that we can scan the source space as Group analysis. Testing the effect prevalence.
We have derived the dimensionality of each subject's subspace using tests described in the previous section. The subspace matrix P i of subjects i = 1, . . ., N has an estimated dimensionality ofD i . IfD i > 0, one can claim that subject i shows the effect of interest, i.e. there is a spatial dimension in which the null hypothesis can be rejected. For example, the effect of interest could be that the power of gamma oscillations in one spatial dimension of the subject's sensor data is stronger during presentation of a stimulus than during rest. As a first step for grouplevel analyses, we test if the mere presence of the effect is generalizable to the population. If the effect is deemed generalizable, we run a procedure to test if there exist consistent source spatial locations across the group, where the effect originates from (see next section).
The first step requires formulating a prevalence hypothesis [see e.g. 10,53,54]. In this framework, a true effect is assumed to be present in a proportion γ of the population. Hence if a subject i is randomly selected from the population We then specify a prevalence null hypothesis that γ is smaller than or equal to a certain proportion γ 0 . In order to claim that the effect is generalizable to the population, an intuitive value for γ 0 is 0.5, i.e. the effect would be present in the majority of the population. If we observed that K out of N subjects showed an effect ðD i > 0Þ, we can define a p-value for the likelihood of K or more out of N subjects showing an effect, if the prevalence across the population is smaller than or equal to γ 0 : If this p-value is below a specified significance level, the effect is deemed generalizable to the population. A certain subject i can show an effect both if the effect is actually present (with sensitivity β), or because of a false positive (at the specified α for the single-subject tests). The probability to pick a subject i from the population that shows an effect, assuming a population prevalence of γ, is thus The probability to pick a subject from the population that shows no effect is Thus the probability to observe K out of N subjects with an effect (see top panel) is pðKjg; NÞ ¼ N K ðgb þ ð1 À gÞaÞ K ðgð1 À bÞ þ ð1 À gÞð1 À aÞÞ The sensitivity β is usually not known, and therefore is fixed at 1, to remain conservative. Computing the p-value as in Eq 42 to test the prevalence null hypothesis requires to sum over these values for K and higher and then to maximize over the range of γ values covered by the null hypothesis pðk ! Kjg g 0 ; NÞ ¼ As discussed in [10], one can also report the largest γ 0 value under which the null hypothesis can be rejected at the given significance level. This can be interpreted as the lower bound of a one-sided confidence interval about the true population prevalence γ, which can be of interest to the research question.
Statistical thresholding of subcorr maps across subjects.
In the previous section we tested the null hypothesis that an effect is present e.g. in only half of the population or less. We did not test the spatial consistency of the effect. This might be sufficient for a given research question: e.g. whether beta-band oscillations can be used to discriminate between two different conditions in a majority of subjects. The inference procedure on the spatial consistency of the effect is described in the following.
In order to derive spatial inferences across the group, we specify a new null hypothesis with respect to the subcorr values as the statistic of interest used to localize effects in individual subjects. The average subcorr value across the group at a source location ρ is computed as where P s i is the signal subspace of subject i in a group of i = 1, . . ., N subjects. Note that P s i consists of the D i first columns (the significant spatial patterns) from the M × M matrix P i . If there is no spatially consistent effect across subjects, we can randomly flip the ordering of spatial patterns in P i and select the D i spatial patterns associated with the smallest quality function values: We refer to this subspace as the noise subspace P n i . Under the null hypothesis, signal and noise subspaces are exchangeable with respect to the average subcorr statistic. We now run O permutations, where at each iteration we draw a binary permutation vector ω of length N. Then we compute the average subcorr value based on shuffled subspaces, where for each subject we use and keep the maximum subcorr Ã (ρ) over the volume at each iteration to obtain a null distribution against which to test the observed subcorrðrÞ values across the volume. Note that there can be subjects that show no effect, i.e. whose signal subspace P s has dimensionality of D = 0 (as discussed in the previous section). Since their subcorr will be zero everywhere (whether using the signal or noise subspace) they will not change the result of the statistical inference. Hence they are omitted in this part of the analysis. We want to emphasize that this does not constitute a case of circular analysis [55], since we use two independent statistics for the two analysis steps. Inference on the quality function value f is used to determine whether a subject shows an effect. This does not tell us anything about the spatial consistency of the effect in the subgroup, or whether we expect higher subcorr values from scanning the signal vs. scanning the noise subspace (which would violate the null hypothesis tested here). These were set to occur at the same time in left/right homologous regions in 50% of the time, such that the resulting amplitude correlation in the alpha band (8-13 Hz) was above 0.4. Panel a) shows such a signal pair. The phase delay between these oscillatory bursts was uniformly random (see enlarged signal parts). We then performed an iES analysis as shown in Fig 8 of the main text (functional connectivity) with a left precentral seed. We can see in panel b) the obtained subspace patterns for a threshold of r > .4. Note that we obtained two spatial patterns that corresponded to the leadfields of the two regions correlated to the seed (precentral left/ right). Panel c) shows that the iES maps revealed the correlated contralateral source with or without the projection step explained in the main text. We can also see that only the correlated sources were present in the map (precentral left/right) and not the sources that had oscillatory bursts in this frequency band (rostralmiddlefrontal left/right) but were not correlated to the seed. (TIF) S3 Fig. Simulation of sources with correlated amplitudes (zero lag). We simulated the same setup as in S2 Fig with the difference that the co-occuring oscillatory bursts in left/right homologous regions had zero phase delay (see enlarged signal parts in panel a). Even though the amplitude correlations were still above 0.4, a linear combination of the leadfields of both sources was captured in one subspace pattern (panel b). This is because the decomposition produces subspace patterns whose corresponding signals are orthogonal. We see in panel c) that the iES map without projection applied did not reveal the contralateral correlated region. However when projecting out the seed topography, the peak in the contralateral region was revealed, as expected. This illustrates the difference between iES and other approaches that orthogonalize signals, which would remove the correlated oscillatory bursts in this scenario. (TIF) S4 Fig. Simulation results analogous to Fig 7, with a superficial and a deep source. a) Two sources of interest (precentral left and parahippocampal left, according to the Desikan-Killiany atlas, [22]) targeting pre-specified power ratios are embedded in background brain noise composed of 1/f signals evenly distributed across 66 locations. The superficial source is simulated at a fixed power ratio (narrowband vs. broadband power) of 0.6, the power ratio of the deep source is varied between 0.2 and 0.6. b) Normalized detection probability is calculated as in Fig 7 and shows that, while the scenario is more challenging for iES (higher power ratio of the deep source is needed to detect both sources), it outperforms MNE, which did not detect both sources in this scenario. (TIF)