The participants in the experiment were 12 healthy adults (aged 25.4–41.9 years, mean 33.7 years; six women). None of the participants reported a history of head injury, neurological disease, hearing problems, severe medical illness, or drug abuse. The experiment was approved by the Ethics Committee of Kyushu University.

The paradigm consisted of auditory stimulus sequences composed of standard stimuli with a probability of 80% and deviant stimuli with a probability of 20%, which were delivered in random order until at least 150 deviant stimuli were presented. Tone bursts of 500 Hz for standard stimuli and 550 Hz for deviant stimuli (10-ms rise and 20-ms fall) with a 100-ms duration were delivered monaurally through plastic tubes (length, 6 m; inner diameter, 8 mm). The hearing threshold was determined for each ear of each subject, and stimuli generated by a tone-burst-generator (Kyushu-Keisokuki, Fukuoka, Japan) were delivered at intensities of 50 dB above the threshold (46). The stimulus onset asynchrony (SOA) was 500 ms, and the presentation rate of the stimuli represented by f_p was 2 Hz. Stimuli were delivered to each ear in separate runs, with masking noises delivered to the contralateral ear to avoid cross-hearing (47). Inversed stimuli (550 Hz for standard and 500 Hz for deviant) were presented monaurally in separate runs. These stimuli were counterbalanced. In the current study, only data from right-ear stimulation and using 500-Hz standard/550-Hz deviant stimuli were analyzed. Subjects were instructed to ignore the auditory stimuli while they lay on the bed and watched a silent movie (16).

MEG was acquired using a 306-channel (204 planar gradiometers and 102 magnetometers) whole-head system (Elekta-Neuromag, Helsinki, Finland) in a magnetically shielded room. The sampling rate was 1,000 Hz, with a band-pass filter of 0.03–330 Hz. EEG was simultaneously recorded using 19 scalp electrodes according to the international 10–20 system, although the sparse EEG data were not analyzed in the current study.

To compare the component distribution patterns between the four methods, a two-way repeated-measures analysis of variance (rmANOVA) was used to analyze the center (z-scored M_max and z-scored C_max, respectively), and slope of the first PCA component with within-subjects factors of APPROACH (subtraction vs. weighted) and DECOMPOSITION (BSS_T/k vs. infomax). For the post hoc tests, multiple comparisons were performed using paired t-tests with Bonferroni correction. The significance level was set at p < 0.05.

We counted c, where the non-linear approximation reached the 5% threshold. It was assumed that components above the 5% threshold significantly contributed to the MMR and were defined as “dominant components,” whereas those that did not meet the threshold did not contribute to the MMR.

The results of the spatio-temporal cluster permutation are shown in Figure 3A. Among the six clusters, Clusters #1 and #2 (20 left temporal sensors at 96–276 ms with the alpha level of p = 0.004, 24 right front-temporal sensors at 105–266 ms with the alpha level of p = 0.008) were consistent with the empirical findings. On the other hand, Clusters #3, #4, and #6 contained the late latency and with lower alpha levels of p = 0.01. The Cluster #5 was at 102–212 ms mainly from parietal sensors with the alpha level of p = 0.008. Therefore, Clusters #1 and #2 were selected to define the reference standard (i.e., a selection of MMR sensors from 44 left temporal and right front-temporal sensors at an MMR time range of 96–276 ms).

Figure 3B represents the grand-averaged ERFs elicited by standard and deviant stimuli and the MMR difference sensor waveforms in sensor space (x¯sub). As indicated by the results of the cluster permutation (Figure 3A topographical map in Clusters #1 and #2), prominent activity occurred in the bilateral temporal and right frontal sensors at the peak latency (Figure 3B topographical map). Note that the offset of the M100 was included in the MMR time for both standard and deviant ERFs (textured arrows in Figure 3B).

Figure 4 represents the results of the decomposition together with the sensor space analysis of a representative subject (Subject 2). The resemblance of the reference standard (Figure 4A) from this subject was compared with each component from four methods (subtraction-BSS_T/k, weighted-BSS_T/k, subtraction-infomax, and weighted-infomax, in Figure 4B(i-iv), respectively). One component in the weighted-BSS_T/k was discriminable with similar morphology [Figure 4B(ii) upper panel, red line, BSS107] and topographical map [Figure 4B(ii) right] to the reference standard of the peak time (Figure 4A, 140 ms). When this component was back-projected into the sensor space (Figure 5, red lines), the left and right temporal sensors (dotted areas) within MMR sensors at MMR time range closely represented the reference standard (blue lines). Accordingly, the corresponding component had discriminable MmaxDC_φ and CmaxDC_φ among other components in the scatter plot [Figure 4B(ii) lower panel, red arrow]. Moreover, this component showed a minor additional topographical representation in the left temporal sensors, which corresponded with the reference standard of 260 ms. No components were discriminable using the infomax methods [Figure 4B(iii, iv)]. The subtraction-BSS_T/k [Figure 4B(i)] provided two components (red and green arrows) that had a moderate value of MmaxDC and CmaxDC.

Scatter plots of M_max and C_max for each component are depicted for the four different methods for all subjects (Figure 6A). While most of components had lower M_max and C_max values in the four methods, in the weighted-BSS_T/k [Figure 6A(ii)], one or a few components represented high MmaxDC_φ and CmaxDC_φ values individually. The z-scored plot shows the distribution of salient (MMR-related) and inconsequential (non-MMR-related) components according to the quadrants based on a 90% z-score (Figure 6B). The salient components were mostly located in the RU quadrant (major) in the weighted-BSS_T/k [Figure 6B(ii)], whereas in the subtraction-BSS_T/k [Figure 6B(i)], they were equally distributed between the RU (major) and RL (pseudo) quadrants. The two infomax methods [Figure 6B(iii and iv)] had salient components mostly in the LU (minor) or RL (pseudo) quadrants. Most components are inconsequential components in all four methods (the numbers of salient components are shown in Table 1).

The component distribution patterns of these salient components (major, minor and pseudo) were further investigated using PCA. The averaged center of the distribution of the salient components and the first PCA component are superimposed on the z-scored plots of salient components in Figure 7 (individual plots are shown in Supplementary Figure 3). The average z-scored M_max, C_max, and slope for the four methods were 3.1 ± 1.4, 0.6 ± 0.4, and 0.29 ± 0.2, respectively, for subtraction-BSS_T/k [Figure 7(i)], 3.7 ± 1.4, 1.6 ± 0.5, and 0.14 ± 0.1, respectively, for weighted-BSS_T/k [Figure 7(ii)], 0.8 ± 0.5, 1.6 ± 0.7, and −0.60 ± 0.4, respectively, for subtraction-infomax [Figure 7(iii)], and 0.7 ± 0.5, 1.6 ± 0.6, and −0.60 ± 0.4, respectively, for weighted-infomax [Figure 7(iv); Table 1].

The rmANOVA results of the z-scored M_max of the center (Figure 8A) revealed a significant main effect of DECOMPOSITION [F_(1,9) = 76.9, p < 0.001], which indicated that the z-scored M_max in both BSS_T/k methods was significantly larger than that in both infomax methods. There was no significant interaction between APPROACH and DECOMPOSITION [F_(1,9) = 2.7, p = 0.1] or main effect of APPROACH [F_(1,9) = 0.9, p = 0.4]. These results suggested that the salient components of both BSS_T/k methods were located in the right quadrant, whereas those of both infomax methods were located in the left quadrant.

The rmANOVA results of the z-scored C_max of the center (Figure 8A) revealed a significant interaction between APPROACH and DECOMPOSITION [F_(1,9) = 60.4, p < 0.001] and significant main effects of APPROACH [F_(1,9) = 40.2, p < 0.001] and DECOMPOSITION [F_(1,9) = 9.7, p < 0.01]. The post hoc analysis revealed that the z-scored C_max of subtraction-BSS_T/k was significantly lower than that of weighted-BSS_T/k and those of both infomax methods (weighted-BSS_T/k, p < 0.0001; weighted-infomax, p < 0.0005; subtraction-infomax, p < 0.001). These results suggested that the salient components of the weighted-BSS_T/k and both infomax methods were located at the border between the upper and lower quadrants, whereas those of the subtraction-BSS_T/k were located in the lower quadrant.

The rmANOVA results of the slope of the first PCA component (Figure 8B) revealed a significant main effect of DECOMPOSITION [F_(1,9) = 65.5, p < 0.001], which indicated that the slope in both BSS_T/k methods was significantly larger than that in both infomax methods. There was no significant interaction between APPROACH and DECOMPOSITION [F_(1,9) = 0.2, p = 0.6] and no main effect of APPROACH [F_(1,9) = 2.4, p = 0.2]. These results indicated that the locations of the salient components in both BSS_T/k methods had positive spatio-temporal correlations regarding the MMR (i.e., the slope had a positive value), whereas those of both infomax methods had negative correlations (i.e., the slope had a negative value).

In conclusion, the distribution of the salient components was mostly in the RU quadrant (major) with weighted-BSS_T/k [Figures 6B(ii), 7(ii)], the RL quadrant (pseudo) with subtraction-BSS_T/k [Figures 6B(i), 7(i)], and the LU (minor) or RL (pseudo) quadrants with the two infomax methods [Figures 6B(iii and iv), 7(iii and iv)]. Both BSS_T/k methods [Figures 6B(i and ii), 7(i and ii)] showed positive spatio-temporal correlations while both infomax methods showed negative correlations [Figures 6B(iii and iv), 7(iii and iv)].

Figure 9 shows the results of M_ave after cumulative back-projection in a representative subject (Subject 2). The curvature of the weighted-BSS_T/k [Figure 9(ii)] was steep in the first component (c = 1, red arrow, corresponding to BSS107), which suggested that in the weighted-BSS_T/k, only one component contributed highly to the MMR. Note that this component was a major component localized on the RU quadrant [Figure 6B(ii)]. On a contrary, other components represent a minimal increase in M_ave [e.g., green and yellow arrows from weighted-BSS_T/k in Figure 9(ii) or all three arrows from two infomax methods in Figure 9(iii and iv)]. These were either pseudo- or minor components (Figure 6B). In addition, it is notable that the third component of subtraction-infomax [Figure 9(iii), yellow arrow, corresponding to pseudo-component in Figure 6B(iii)] negatively contributed to the MMR. Moreover, the first component of subtraction-BSS_T/k [Figure 9(i), red arrow] showed a mild increment in M_ave, which corresponds to this component being classified as a major component in Figure 6B(i).

Supplementary Figure 4 shows the RC lines (upper panels) together with their approximate lines (lower panels) in individual subjects. In Subject 2, the approximate lines of the weighted-BSS_T/k show that the first component (red arrow) represented a contribution as high as 30%, whereas later components (green and yellow arrows) provided much lower contributions. We counted c, where the non-linear approximation reached the 5% threshold (gray dotted lines; i.e., the dominant components). The number of dominant components is shown in Table 2. In the weighted-BSS_T/k, 1–3 components significantly contributed to the MMR, except for one subject (Subject 5). In the subtraction-BSS_T/k, 2–6 components contributed to the MMR. The two infomax methods had few components that significantly contributed to the MMR. These results indicated that one or a few dominant components contributed to the MMR in weighted-BSS_T/k, whereas no components represented the MMR in infomax.

The MMR has been widely used in many fields of human neuroscience (10, 15, 16). Conventionally, the subtraction approach was needed to extract the MMR from other auditory ERP/ERF. However, there are several problems with such a method, which include increased noise and the inability to exclude neural adaptation. Several approaches have been proposed to avoid the neural adaptation problem (53, 57); however, all such approaches depend on subtraction. Our novel approach avoids subtraction. In general, the MMR is a relative component because a common response is included in standard and deviant ERFs, and the MMR is then defined as the difference waveform based on the original theory underlying the MMR (i.e., the memory-comparison process). The MMR should be present in deviant epochs but not in standard epochs. Thus, only deviant epochs are needed for its decomposition.

We made two assumptions underlying the successful decomposition of the weighted-BSS_T/k: (1) The MMR occurs periodically within a specific time range (i.e., the MMR time range) and in the deviant epochs; (2) Exogenous/obligatory ERFs highly correlate with themselves in the non-MMR time. Originally, BSS_T/k was expected to highlight periodic signals using T (35–38). The MMR time range (96–276 ms) was defined according to the spatio-temporal cluster permutation analysis, which was assumed to reveal the statistically significant time range in which the MMR occurs. Both the offset response of the M100 and the MMR fall into this time range, whereas the onset response of the M100 occurs outside of the time range (Figure 3B). Assigning a weight to this time range may minimize the joint M100 and MMR effect. The weighting emphasizes the target response (i.e., the MMR) within the window, whereas the response outside the window (i.e., the onset of the M100) takes away the response (i.e., the offset of the M100) if they are highly correlated. Analogous to the subtraction approach (as subtraction separates such responses by subtracting deviant responses from standard responses), the weight assignment on a specific time range may differentiate the MMR from other responses. Supplementary Data (2-1) and Supplementary Figure 5 support our assumptions; the assignment of the weight outside the M100 in the standard epochs did not result in the extraction of remarkable components that represent the M100 (Supplementary Figure 5B).

We obtained four main findings. First, the weighted-BSS_T/k decomposed one or a few components (<3) that manifested the MMR among the many components obtained from multi-channel data (Figures 4, 6, 9, Supplementary Figure 4, and Table 2). We refer to this decomposition result as specification. Multi-channel recordings of electromagnetic fields emerging from neural currents in the brain generate large amounts of data (28). Thus, this specification makes interpretation and comparisons among groups easier. Our primary aim was to extract MMR in a few dominant components. The dominant component was the component that had the most discriminable M_max and C_max, and thus, it contributed most highly to the MMR (Figure 9, Supplementary Figure 4, and Table 2). We do not assume that the dominant component manifests a single MMR source; instead, it may represent the network or a series of MMR sources (Figure 4). Other irrelevant activities were redistributed among the remaining components. Since our method (BSS_T/k) depends on the theory that utilized correlations between components instead of strong independence (i.e., ICA), it would result in extracting components with keeping physiological correlation that may represent several generators or network of MMR. If bitemporal and frontal MMR sources are highly correlated, with a certain delay, these sources should be extracted in a few components using our time-delayed correlation method. Indeed, it is known that these sources have separate temporal dynamics (58) but interact with each other (59). In contrast, it is difficult to identify any dominant components using ICA, where each extracted component represents one or two dipolar sources. This is discussed in the following section.

Second, the decomposed components revealed positive spatio-temporal correlations regarding the MMR, and the center of the distribution of the salient components was in the RU (major) quadrant (Figures 6–8, Supplementary Figure 3, and Table 1). According to Eq. (2), the decomposed component contains the mixing matrix (spatial) and signal source (temporal). A positive spatio-temporal correlation in the decomposed component suggests that the component is physiologically meaningful (9). In turn, with a positive spatio-temporal correlation, a component that shows the most similar morphology regarding the MMR also has the most similar topography regarding the MMR. This relationship is particularly important when targeting the response with an unknown generator source. The temporal information can be mutually applicable to the detection of the target, without a priori knowledge of its precise generator. For example, in Figure 4B(ii), if the MMR topography is unknown, BSS107 can be selected as the MMR component based on its discriminable amplitude.

Third, each component was obtained from individual data and the results were statistically significant. This indicated that weighted-BSS_T/k is generally applicable to individual subjects, unlike group-ICA.

Fourth, a new cohort from subjects with low SNR in the sensor-space analysis regarding MMR (subsection Spatio-Temporal Cluster Permutation to Define the MMR Time Ranges and Sensors or the Reference Standard) demonstrated a few MMR-related components in weighted-BSS_T/k when the same MMR time range was used for the weight assignment [Supplementary Data (2-2) and Supplementary Figure 6]. This MMR time range was independently determined in this cohort. These results may indicate that the generous setting of the weight time range can be available as long as the crucial time range is covered.

Based on these results, the application of our approach provides potential benefits that the sensor-space subtraction method does not, despite its status as the current gold standard for revealing MMR. Our single-trial, contrast-free approach would minimize the effect of refractoriness and maximize the temporal information underlying the neural mechanism of MMR. Our approach would therefore provide a new approach toward investigating further insights into the physiology of MMR.

Both ICA methods (subtraction- and weighted-infomax) consisted of a collection of minor or pseudo components (Figures 4, 6, 9, Supplementary Figure 4, and Table 2). Most components were located in the LU (minor) or RL (pseudo) quadrants (Figures 6–8 and Table 1). The slope of the first PCA component showed a negative spatio-temporal correlation (Figures 7, 8, Supplementary Figure 3, and Table 1). There were no dominant components that manifested the MMR in either of the ICA methods (Table 2). The decomposition method in ICA is based on stochastic properties and does not depend on the time structure; thus, spatio-temporal dissociations may occur (34). Several papers have reported successful decomposition of the MMR using ICA (9, 25, 60–64); however, most results were derived from oligo-channel recordings. When the number of sensors/channels sensing the MMR is relatively small, the MMR can be extracted by one or a few components. However, such specification in multi-channel data is rarely shown in ICA studies because a greater number of channels results in poorer estimation accuracy of the components (25). If we assume fewer numbers of sources (e.g., tens) but use larger numbers (204) of sensors for ICA decomposition, the components of interest will likely be (i) split into sub-components and (ii) located where the SNR of each component is reduced. This is consistent with our previous work where ICA decomposition showed fragments of interictal epileptiform discharges from a single epileptogenic zone (42). Furthermore, most ICA studies are based on cluster analysis (e.g., group-ICA), not individual analysis. Generalization of the application of ICA to the MMR was not demonstrated in these studies.

Lastly, although subtraction-BSS_T/k follows the conventional subtraction approach, it performed better than the two ICA methods, especially the subtraction-infomax. The center of the distribution of the salient components was in the RL quadrant (pseudo), yet it maintained a positive spatio-temporal correlation (Figures 6–8, Supplementary Figure 3, and Table 1). A possible interpretation of these findings is that these components may represent the partial generators of MMR sensors. The difference between subtraction-BSS_T/k and subtraction-infomax may explain the theoretical difference between BSS_T/k and infomax (time-delayed correlation vs. strong independence). The decomposition of the subtraction-BSS_T/k was less successful than that of the weighted-BSS_T/k. There were more dominant components (<6; Table 2) in the subtraction-BSS_T/k than there were in the weighted-BSS_T/k. From the viewpoint of specification, fewer dominant components are desired. In conclusion, both BSS_T/k methods, which use time structure, performed well in extracting the MMR; however, the weighted approach was the most sensitive.

The current study aimed to extract the MMR as a distinct component using a combination of the periodical arrangement and assignment of a weight. The specific effect of each technique should be investigated in a future study, which may help achieve a better understanding of the physiology of the MMR.

Because there was no confidence in terms of source localization of extracted components, although there are several ICA and SOBI studies (9, 32, 41), this view may provide potential benefits given that components may encompass several sources or networks of MMR. This should be investigated in future studies.

Our method is not dependent on the number of components. Our motivation was not to apply dimension reduction to maximize the multi-channel MEG data. However, the application of our method to different numbers of sensors, different MEG systems, or another type of sensor (magnetometer) is an interesting but open question. Theoretically, our weighted method can possibly be applied to any clinical neurophysiology data to investigate ERFs, which include higher cognitive functions where the elicitation of the target requires subtraction, and the target is subject to a specific assumption about the time window in which it occurs in multi-channel data. In a paradigm where stimuli are jittered and thus are not periodic, our weighted method will also be applicable by concatenating the epochs.