{'fraction_non_alphanumeric': 0.04371966277782509, 'fraction_numerical': 0.01599250617389083, 'mean_word_length': 4.809439002671416, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Weighted Blind Source Separation Can Decompose the Frequency Mismatch Response by Deviant Concatenation: An MEG Study

The mismatch response (MMR) is thought to be a neurophysiological measure of novel auditory detection that could serve as a translational biomarker of various neurological diseases. When recorded with electroencephalography (EEG) or magnetoencephalography (MEG), the MMR is traditionally extracted by subtracting the event-related potential/field (ERP/ERF) elicited in response to “deviant” sounds that occur randomly within a train of repetitive “standard” sounds. However, there are several problems with such a subtraction, which include increased noise and the neural adaptation problem. On the basis of the original theory underlying MMR (i.e., the memory-comparison process), the MMR should be present only in deviant epochs. Therefore, we proposed a novel method called weighted-BSST/k, which uses only the deviant response to derive the MMR. Deviant concatenation and weight assignment are the primary procedures of weighted-BSST/k, which maximize the benefits of time-delayed correlation. We hypothesized that this novel weighted-BSST/k method highlights responses related to the detection of the deviant stimulus and is more sensitive than independent component analysis (ICA). To test this hypothesis and the validity and efficacy of the weighted-BSST/k in comparison with ICA (infomax), we evaluated the methods in 12 healthy adults. Auditory stimuli were presented at a constant rate of 2 Hz. Frequency MMRs at a sensor level were obtained from the bilateral temporal lobes with the subtraction approach at 96–276 ms (the MMR time range), defined based on spatio-temporal cluster permutation analysis. In the application of the weighted-BSST/k, the deviant responses were given a constant weight using a rectangular window on the MMR time range. The ERF elicited by the weighted deviant responses demonstrated one or a few dominant components representing the MMR that fitted well with that of the sensor space analysis using the conventional subtraction approach. In contrast, infomax or weighted-infomax revealed many minor or pseudo components as constituents of the MMR. Our single-trial, contrast-free approach may assist in using the MMR in basic and clinical research, and it opens a new and potentially useful way to analyze event-related MEG/EEG data.

The mismatch response (MMR) is thought to be a neurophysiological measure of novel auditory detection that could serve as a translational biomarker of various neurological diseases. When recorded with electroencephalography (EEG) or magnetoencephalography (MEG), the MMR is traditionally extracted by subtracting the event-related potential/field (ERP/ERF) elicited in response to "deviant" sounds that occur randomly within a train of repetitive "standard" sounds. However, there are several problems with such a subtraction, which include increased noise and the neural adaptation problem. On the basis of the original theory underlying MMR (i.e., the memory-comparison process), the MMR should be present only in deviant epochs. Therefore, we proposed a novel method called weighted-BSS T/k , which uses only the deviant response to derive the MMR. Deviant concatenation and weight assignment are the primary procedures of weighted-BSS T/k , which maximize the benefits of time-delayed correlation. We hypothesized that this novel weighted-BSS T/k method highlights responses related to the detection of the deviant stimulus and is more sensitive than independent component analysis (ICA). To test this hypothesis and the validity and efficacy of the weighted-BSS T/k in comparison with ICA (infomax), we evaluated the methods in 12 healthy adults. Auditory stimuli were presented at a constant rate of 2 Hz. Frequency MMRs at a sensor level were obtained from the bilateral temporal lobes with the subtraction approach at 96-276 ms (the MMR time range), defined based on spatio-temporal cluster permutation analysis. In the application of the weighted-BSS T/k , the deviant responses were given a constant weight using a rectangular window on the MMR time range. The ERF elicited by the weighted deviant responses demonstrated one or a few dominant components representing the MMR that fitted well with that of the sensor space analysis using the conventional subtraction approach. In contrast, infomax or weighted-infomax revealed many minor or pseudo components as constituents of the MMR. Our single-trial, contrast-free approach may assist in using the MMR in basic and

INTRODUCTION
The mismatch negativity component in electroencephalography (EEG), and its magnetoencephalographic (MEG) counterpart the mismatch field (or mismatch response, MMR), are eventrelated responses (EPRs/ERFs) widely used to measure auditory processing in cognitive neuroscience (1)(2)(3)(4)(5)(6). The MMR is recorded using an oddball paradigm, where the repeated presentation of a stimulus (standard) is occasionally replaced by a different stimulus (deviant). The MMR is then computed as the difference between the deviant and standard responses. This difference representing the MMR is typically found around 100-250 ms after the onset of the deviant stimulus (7). Previous studies have revealed a cortical network consisting mainly of the bilateral temporal regions, but also the frontal and parietal regions, which is involved in the generation of the MMR (8)(9)(10). The prevailing view is that the MMR reflects the detection of change in the auditory system that can be measured without attention, although alternative interpretations exist (11)(12)(13)(14). The MMR has therefore been widely used to assess auditory processing in children and clinical groups (10,15,16).
Originally, it was suggested that the occurrence of the MMR relates to the presence of a short-term memory trace where the memory-comparison process detects a discrepancy between the neural representation of the regularity inherent in the recent stimulation and the representation of the current deviant stimulus (17). On the basis of this hypothesis, obtaining a difference waveform by subtracting the standard response from the deviant response is the only way to identify the MMR. However, there are several problems associated with the subtraction approach. First, the subtraction reduces the signalto-noise ratio (SNR) because the noise present in the standard responses is added to the noise in the deviant responses. Second, the neural adaptation process, especially with frequency MMR, can affect the difference waveform. The auditory system has a tonotopic organization from the cochlea through to the cortex (18). Stimulus repetition leads to repeated initiation of patterns of neural activity (e.g., the M100) that habituates as a function of the repetition rate (19,20). In the classic oddball protocol, the neural response to standard stimuli is attenuated by these repetition suppression effects. This suppression is greater for the standard stimuli than for the less frequent deviant stimuli. The adapted and non-adapted neural activity presents not only different amplitudes, but also different temporal dynamics. Thus, the subtraction approach does not simply reflect the MMR (i.e., a memory-based comparison) but also the differential adaptation of neurons (13). Therefore, the study of the temporal dynamics of the MMR might convey critical information regarding the nature of the underlying neural generators. Hence, to effectively reveal the MMR, another approach considering the temporal information, instead of the subtraction approach, is desirable.
Each EEG electrode or MEG sensor records a linear combination of signals from several sources (21). Multichannel EEG/MEG, which typically involves hundreds of sensors, provides detailed spatio-temporal distribution patterns, which obviously complicate the interpretation of signals and topographies. Independent component analysis (ICA), which is a blind source separation (BSS) method, is a stochastic method that can be used to decompose such complex data into a set of spatiotemporal components, each of which comprises a fixed spatial distribution and an associated signal (22,23). Each component signal is a weighted sum of the sensor or electrode signals, which in turn are weighted sums of the dynamics of the neural sources (24). ICA/BSS can provide signal sources without any a priori information about their occurrence in biological signals. In general, the single-trial approach of ICA/BSS can utilize temporal information, because the contraction of information occurs during the averaging process of the ERP/ERF. A single trial may contain all kinds of non-brain artifacts and spontaneous EEG/MEG processes, whereas decomposing an average of all trials not only minimizes the contributions of those neural and artifactual processes that are not reliably time-and phaselocked to experimental events but also removes event-related brain dynamics among trials (25). As artifacts often exhibit stereotypical patterns that differ from those of brain activity, ICA/BSS can mostly be used to separate artifactual patterns (26)(27)(28). In fact, ICA/BSS has been used to extract event-related activities in only a handful of previous studies (29)(30)(31)(32)(33). Owing to the components being computed based purely based on their statistical independence, physiological perspectives are not taken into account (28,34). Considering that regional brain activities substantially correlate with each other, an approach requiring strong independence may not be the most fruitful (35)(36)(37)(38).
An approach for refining ICA/BSS using time-delayed correlation, or the decorrelation method (DC) has also been considered (39). Time-delayed correlation takes account of the characteristic time structure of the signals of interest, including the periodicity and/or morphology. Thus, time-delayed correlation measures the correlation between two signals, then maximizes the correlation between components. For example, several studies applied DC to second-order blind identification (SOBI) to separate periodic signals, such as cardiac and oscillatory brain activity, because periodic signals are wellcorrelated with delayed signals and non-delayed original signals (40,41). As a result, well-correlated signals were extracted in one (or a few) components. However, in SOBI, most approaches examine the time structure of the target signals subjectively. When parameters are highly specified, featured components are more independent, and therefore target signals collapse because of strong independence and the SOBI method becomes equivalent to ICA (39)(40)(41).
In an attempt to develop solutions to address the limitations of ICA and SOBI, we proposed a novel method of BSS called the T/k (fractional) type of DC (BSS T/k ) (35)(36)(37)(38). This method shares the fundamental concept underlying DC such as SOBI but is more focused on the periodicity of the target signal. The BSS T/k method is based on extracting time points (i.e., time-delayed parameters) determined by the parameters T and k, which represent periodicity concerning a fundamental and harmonics (Supplementary Data (1) and Supplementary Figure 1). BSS T/k allows weak independence among the components. Setting timedelayed parameters in this way results in highlighting the characteristics of target ERFs that are periodically presented. Previously, we demonstrated that somatosensory-evoked fields in response to periodic electrical stimuli can be decomposed into a few components using the BSS T/k algorithm in 64 channel magnetometers of CTF (35)(36)(37)(38). Using a generalization of BSS T/k , non-periodic interictal epileptiform discharges that were assumed to originate in a single epileptogenic zone were decomposed into one dominant component (42).
For the MMR paradigm, where deviant stimuli are presented in random order, we proposed to use a modification of BSS T/k , which we termed weighted-BSS T/k (43). In weighted-BSS T/k , we only used deviant responses that were concatenated into a periodical arrangement. Then, deviant responses were assigned a constant weight (rectangular window) on the specific time interval that represents MMR (i.e., the MMR time range). This is known as a window function in the time domain. The MMR time range was defined in a data-driven manner using sensor space subtraction (i.e., the reference standard). Through these procedures, the correlation between MMR and the responses outside of the MMR time range (e.g., the M100) can be minimized; thus, weighted-BSS T/k , which underlies time-delayed correlation, can effectively extract the MMR. We hypothesized that weighted-BSS T/k would extract one or a few dominant components that can discriminate the MMR from background brain noise and other artifacts or other irrelevant ERFs. As the first application in the cognitive neuroscience of weighted-BSS T/k using only deviant epochs, we aimed to extract components that resemble the reference standard because subtraction is currently the gold standard for identifying MMR. We applied both BSS T/k and infomax (ICA) separately to the same weighted multi-channel MEG data (weighted-BSS T/k and weighted-infomax, respectively), and used the subtraction approach (subtraction-BSS T/k and subtraction-infomax) as a more general approach to investigate how the single-trial approach works, and then, statistically compared the similarity of each component to the reference standard to test a further hypothesis that BSS T/k is more sensitive than infomax.
It was not our aim to use the subtraction-BSS T/k /weighted-BSS T/k to separate independent MMR sources. Typically, statistically independent components separated by preprocessing with ICA are expected to be associated with one or two dipolar sources (9,23,44,45). We instead made a more general assumption that a component extracted by subtraction-BSS T/k /weighted-BSS T/k will relate to multiple sources or a network of activity generating the MMR. In this sense, few decomposed components are better than many, as long as they represent the reference standard. Thus, the extraction of MMR in a few components would simplify the interpretation of MMR in regard to clinical and research applications.

Participants
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.

Stimuli and Procedures
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).

Data Acquisition
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.

Preliminary Process
The temporal signal space separation method (TSSS) using MaxFilter 2.2.13 (Elekta-Neuromag, Helsinki, Finland) was applied to the sensor level data with the default setting of an inside expansion order of 8, an outside expansion order of 3, automatic optimization of both inside and outside bases, a subspace correlation limit of 0.980, and a raw data buffer length of 10 s (48,49). Notch filters were applied to suppress power line frequency and its harmonics (60,100,120,180,200,240, and 300 Hz). Data from the 204 planar gradiometers were used for all subsequent analyses. Hereafter, all analysis steps are shown in Figure 1 and summarized in Figure 2.

The Subtraction Approach
The conventional subtraction approach for sensor space analysis was used as a reference (Figures 1A, 2A). Before averaging across epochs, the data were low pass filtered at 30 Hz and, epochs exceeding 4,000 fT/cm on any planar gradiometer channel were excluded from the average. Based on our experience, some ocular artifacts leak into the good epochs. Therefore we took extra care and visually inspected the data to remove eye movements. However, the impact of this procedure was minimal because the number of epochs removed for each subject was 0-1. Each epoch contained a 600-ms time window ranging from 100 ms pre-stimulus to 500 ms post-stimulus onset, with the stimuli being periodically presented (SOA = 500 ms or f p = 2 Hz). The MMR difference sensor waveform (i.e., x sub ) was calculated by subtracting the averaged deviant ERFs from the averaged standard ERFs for each subject (Figures 1A, 2A); where x (n) represents the MEG sensor data at the discrete time, n. x (n) reflects the averaged sensor waveform of x (n) across epochs. x std (n) and x dev (n) are averaged standard and deviant responses, respectively.

Decomposition Process
The decomposition methods of BSS T/k and infomax were applied separately to each subject's sensor dataset, which contained 204 sensors. The sensor data were originally decomposed into a set of spatio-temporal components; where A is a mixing matrix, and s is a signal source. BSS T/k was applied; where A DC is a mixing matrix of BSS T/k , and s DC is a signal source of BSS T/k . Hereafter, we refer to BSSq (q = 1, 2, 3, . . . , 204) as a specific component obtained after the application of BSS T/k . We briefly describe the BSS T/k method here; full details are provided in previous studies (36,38). As a preliminary step, we conducted a sphering procedure to orthogonalize and normalize the time-series data for input sensors. We then conducted an iterative Givens rotation to minimize the absolute sum of offdiagonal elements of the normalized correlation matrices at the parameters. Specifically, the Jacobi-like algorithm proposed by Cardoso and Souloumiac (50,51) was used in the BSS T/k method to approximately solve the simultaneous diagonalization problem at specific times. Regarding the period T = 1/f p with sampling frequency f s , the time-delayed parameters τ can be defined by: For ICA, we used the infomax algorithm (25,49), which was implemented in MNE-python (52) using the default setting; where A ICA is the mixing matrix of infomax, and s ICA is the signal source of infomax. Hereafter, we refer to ICAq (q = 1, 2, 3, . . . , 204) as a specific component obtained after application of infomax. The number of principal components from the pre-whitening step that was passed to the ICA algorithm was 204, which corresponded with the number of sensor inputs. Accordingly, we obtained 204 components with associated time courses and spatial distributions.

Two Different Approaches (Subtraction and Weighted)
After applying the decomposition methods (BSS T/k and infomax) to the sensor space data, we obtained the MMR difference source waveform (i.e., s sub ; Figures where s std (n) and s dev (n) are the averaged source waveforms across epochs (i.e., ERFs) elicited by the standard and deviant stimulus, respectively, obtained from each decomposition method. The novel method, the weighted-BSS T/k , is expected to be a more sensitive approach of extracting the MMR. The basics of the method lie in the periodical arrangements and assignments of weights on the MMR time range. Although our BSS T/k method is expected to highlight periodic signals, the deviant epochs occur randomly, not periodically. To obtain periodical The conventional subtraction approach for sensor space analysis. Sensor data (x) consist of standard (std) and deviant (dev) epochs. The MMR difference sensor waveform (x sub ) is calculated by subtracting the event-related field (ERF) to the standard (x std ) stimulus from the ERF to the deviant (x dev ) stimulus. BSS T/k , fractional type of decorrelation method; MMR, mismatch response; SOA, stimulus onset asynchrony; f p , presentation rate of stimuli. (B) The subtraction approach with two decomposition methods (BSS T/k or infomax). Data from 204 sensors (x) are decomposed into 204 components (s). The subtraction approach is followed for each component; subtraction of the deviant ERF (s dev ) from the standard ERF (s std ) makes the MMR difference source waveform (s sub ). Note that deviant epochs occur randomly, not periodically. (C) Deviant epochs of sensor data are concatenated and new raw data are made (x dev ). This process makes periodical arrangements of deviant epochs. (D) Sensor data (x ϕ ) are assigned with a weight on the MMR time (from n 1 to n 2 , highlighted in purple shadows) on x dev using a window function. The window function is shown in the inset figure. (E) Data with 204 sensors assigned with a weight (x ϕ ) are decomposed into 204 components (s ϕ ) with two decomposition methods. In each component, the ERF of the deviant epochs assigned with a weight is obtained (s ϕ ). The inverted black arrows between A and B and between D and E represent the back-projection process in a group of several components. arrangements, we concatenated the deviant epochs to form new raw data (x dev (n); Figure 1C). To highlight the MMR that was included in the deviant epochs, we then weighted the MMR time range (around 100-250 ms, from n 1 to n 2 ) defined by the spatio-temporal cluster permutation (subsection Spatio-Temporal Cluster Permutation to Define the MMR Time Ranges and Sensors or the Reference Standard), with the weight described by the window function of the rectangular window (inset between Figures 1C,D); where ϕ describes a window function and the * reflects its repeat operation. The segmentation of data (epoch number, mean 174.3 ± 19.6 [standard deviation]) was multiplied by the window function values. Equation (8) indicates, where n = Index of deviant epoch − 1 · SOA + m. Here, m is the given time point within every deviant epoch. Equation (8) indicates that this window function, Eq. (9), was applied repeatedly ( Figure 1D, purple shadow) to the concatenated sensor data (x dev (n)). We then applied the BSS T/k and infomax methods separately to the weighted data (weighted-BSS T/k and weighted-infomax; Figure 1E); Finally, after lowpass filtering (30 Hz), we obtained the ERFs (i.e., s ϕ ; Figures 1E, 2C). That is, s ϕ DC (n) and s ϕ ICA (n), elicited by the weighted deviant stimulus, instead of subtraction.
Two assumptions underlie the successful decomposition of the weighted-BSS T/k . First, the MMR occurs in the MMR time (n 1 ≤ n ≤ n 2 ) only in deviant epochs. Second, exogenous/obligatory ERFs (e.g., the M100) highly correlate with themselves in the non-MMR time (n < n 1 , n 2 < n). The offset response of the M100 often intrudes on the MMR within the MMR time, which is one of the reasons why the subtraction approach is necessary (53). To minimize the joint M100 and MMR effect, a rectangular window in the non-MMR time is used to keep the correlation of the offset and onset of the M100 and extract these as distinct components from an MMR component using weighted-BSS T/k , which underlies time-delayed correlation. However, it is expected that weightedinfomax, in contrast to weighted-BSS T/k , does not decompose the MMR effectively because infomax does not depend on time structure.

Spatio-Temporal Cluster Permutation to Define the MMR Time Ranges and Sensors or the Reference Standard
Currently, the only way to identify MMR is via sensor-space subtraction. We therefore used sensor-space subtraction as a reference standard. A data-driven approach was used to find significant MMR time ranges and sensors in all subjects. Among the 12 subjects, two did not exhibit a prominent MMR during the initial screening of the visual inspection of sensor space subtraction (confirmed by three independent inspectors, TMat, SK, and KK.) and were thus excluded from further analysis.
Individual MMR difference sensor waveforms, x sub , were tested if they were different from 0 across the 10 subjects, with the multiple comparison problem being addressed using a clusterlevel permutation test across space and time (54). We used 1,024 permutations, and the cluster-defining threshold was set at p = 0.01. Selected samples were clustered based on both spatial and temporal adjacency (i.e., spatio-temporal cluster permutation). Our motivation to use the spatio-temporal cluster permutation method was to verify the empirical knowledge that MMR occurs around 100-250 ms in the bilateral front-temporal sensors (7,17) in a data-driven manner in our cohort of 10 subjects. Figure 3A demonstrates the results of the spatio-temporal cluster permutation. Six clusters (less than the critical alpha level of 0.05) were found. Among these six clusters, two (#1 and #2) contained temporal and/or frontal sensors within approximately 100-250 ms; one (#1) contained 20 left temporal sensors at 96-276 ms and the other (#2) contained 24 right front-temporal sensors at 105-266 ms. Thus, we defined the MMR time range The results of the spatio-temporal cluster permutation analysis. Six significant spatial and temporal clusters are shown in white circles within the averaged t-statistics (absolute value) and in orange shading within the averaged waveforms, respectively. Blue lines, standard; green lines, deviant; red lines, difference. (B) Grand-averaged ERFs elicited by standard (x std ) and deviant (x dev ) stimuli and the MMR difference sensor waveform (x sub ) from 10 subjects. The MMR time is indicated by the purple line in the standard and deviant ERFs and by the blue shading in the MMR difference sensor waveform. The red and pink lines in the MMR difference sensor waveform represent the MMR sensors from the left and right clusters, respectively. The topographical map represents the peak activity in the bilateral temporal and right frontal sensors (white circles). The onset of the M100 (arrows) is outside of the MMR time range, whereas the offset of the M100 (textured arrows) is included in the MMR time range seen in standard and deviant responses.
as 96-276 ms (n 1 = 96, n 2 = 276) and the MMR sensors as these 44 sensors. The reference standard was defined individually (Figure 2A); where x Ref ∈ R L×(n 2 −n 1 ) and F ∈ R L×N is the matrix that select L = 44 rows corresponding to the MMR sensors out of x sub containing all N = 204 sensors. In other words, the reference standard was the 44 MMR sensors selected from the 204 gradiometers within the MMR time.
We confirmed that the different setting of the cluster-defining threshold (p = 0.005) gave the similar spatio-temporal clusters (Supplementary Figure 2). This means that the clusters obtained were robust.

Component Evaluation: Cosine Similarity
To investigate the resemblance of each component to the reference standard individually, or goodness of fit, we measured cosine similarity (C) as spatial similarity and morphology similarity (M) as temporal similarity.
Cosine similarity refers to the similarity between two column vectors (42,55); whereâ(n) = a(n)/ |a| is the normalized column vector containing the spatial distribution of the reference standard (x Ref ), andb = b/ b is a normalized column vector of A in Eqs. (3, 5, 10, and 11). The symbol T is the transpose ofâ(n). Because of its definition, 0 ≤ C a (n) , b ≤ 1. In the following, we used the maximum of C (C max ) across the MMR time range for the four methods (i.e., subtraction-BSS T/k , subtractioninfomax, weighted-BSS T/k , and weighted-infomax), denoted by C DC max , C ICA max , C DC_ϕ max , and C ICA_ϕ max . C max represents how maximally similar each component is to the reference standard in regard to spatial information.

Component Evaluation: Back-Projection and Morphology Similarity
Temporal similarity should include information about the temporal correlation between each component and the reference standard as well as the amplitude difference between each component and the reference standard. Because the components derived from BSS T/k and infomax (M DC (n) , M ICA (n) , M ϕ DC (n) , and M ϕ ICA (n)) are differently normalized, their ERFs cannot be directly compared according to their amplitudes. Thus, each component was projected back into the sensor space (back-projection) (56). Here, we assumed a general situation for the sake of the following subsection The Cumulative Back-Projection of Salient Components, the cumulative back-projection. When a group of q components, where Q = {q} is selected from 204 components, provides back-projected data in the sensor space (inverted black arrow between Figures 1A,B), where A q # ∈ R 204×q and s q # (n) represents source vectors corresponding to Q. Here, the suffix symbol # indicates DC or ICA. The same formula was applied to the weighted data (inverted black arrow between Figures 1D,E). The ERF was then computed using the subtraction or weighted approach. For the subtraction approach, we applied where x #_std (n, q) and x #_dev (n, q) are ERFs in the sensor space elicited by standard and deviant stimuli, respectively, obtained from each decomposition method (DC or ICA). For the weighted approach, x ϕ # n, q is the ERF obtained from each decomposition method (DC or ICA). Then, corresponding to Eq. (12), we applied x #_sub (q) : = F x #_sub n, q , where x #_sub q and x ϕ # ∈ R L×(n 2 − n 1 ) .
We investigated the correlation between one sensor and the reference standard; where (X, Y) is the inner product. Here, X is one row vector (l) of the reference standard (x Ref ), which corresponds to one sensor, and Y is one row vector (l) of the same sensor of Z, where Z(q) is defined as Eq. (16) or Eq. (17). Notably, x #_sub n, q and x ϕ # n, q ∈ R L×SOA and Z(q) ∈ R L×(n 2 −n 1 ) . Equation (18)   With successful decomposition, it is expected that only a few components will fall within the RU quadrant, and the rest of the components will fall within the LU, RL, and LL quadrants near the borderlines of coordinate origin. In contrast, unsuccessful decomposition will provide a component distribution pattern where no components fall within the RU quadrant, and all components will fall near the LL quadrant. To investigate the distribution pattern of the salient components, principal component analysis (PCA) was applied to z-scored M max and C max . Two individual PCA components were obtained, with most of the variance being captured by the subspace of the first PCA component (more than 84%; Table 1). The center of the distribution of salient components, taken as the cross-point of the first and second PCA components, and the slope of the first PCA component were obtained.
If a z-score > 1.96 (95%) was set, the number of salient components was small ( Table 1), especially in the weighted-BSS T/k . PCA seemed unreliable when the input data were <5; thus, a z-score > 1.65 (90%) was applied.

The Cumulative Back-Projection of Salient Components
To investigate the contribution of each component to the MMR, components were cumulatively projected back into sensor space (subsection Component Evaluation: Back-Projection and Morphology Similarity), and the spatio-temporal resemblance was compared with the reference standard (subsections Component Evaluation: Cosine Similarity and Component Evaluation: Back-Projection and Morphology Similarity). It is expected that the more components that contribute to the MMR are cumulatively back-projected, the more the back-projected sensors resemble the reference standard. The order of cumulation was determined after sorting by the first PCA component axis (Supplementary Figure 3). Salient components were selected for cumulative back-projection because components below thresholds (inconsequential components in the LL quadrant) are expected to contribute little to the MMR. Z q all = F x #_sub n, q all , n 1 ≤ n ≤ n 2 , Thus, the denominator of Eq. (21) represents M ave of the reference standard. Corresponding to each method, RC becomes RC DC , RC ICA , RC DC_ϕ , and RC ICA_ϕ .
As it is expected that RC will decrease as c increases, we applied the exponential function approximation to the plotted data in the c-RC plane; where x implies c, and y represents RC with coefficients α and β.

Statistics
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.

RESULTS
The analysis comprised four parts (Figure 2, double squares): (i) defining the reference standard based on the spatiotemporal cluster permutation from the sensor-space analysis; (ii) qualitative evaluation of each component based on its similarity to the reference standard; (iii) statistical assessment of component distribution patterns with the z-scored scatter plot; and (iv) the relative contribution of each component.

Spatio-Temporal Cluster Permutation and Reference Standard
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 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).

Statistical Assessment of the Component Distribution Pattern
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  , and (iv) for subtraction-BSS T/k , weighted-BSS T/k , subtraction-infomax, and weighted-infomax, respectively. Upper panels: source waveforms (s sub or s ϕ ), lower panels: scatter plots of M max and C max for each component. The color map in the scatter plot indicates the value of C max (from 0 to 1). In each decomposition result, three components are depicted in different colors (red, green, and yellow) with their corresponding topographical maps. In the scatter plots, the arrows with the same color correspond to the components. These components were selected based on the order of cumulation of the first three components (see subsection The Cumulative Back-Projection of Salient Components). Hereafter, the topographical map takes an arbitrary unit due to matrix A in Eq. (2). 0.5, 1.6 ± 0.6, and −0.60 ± 0.4, respectively, for weightedinfomax [ 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 spatiotemporal 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 backprojection 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 subtractioninfomax [ Figure 9(iii), yellow arrow, corresponding to pseudocomponent 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).

The Cumulative Back-Projection and Relative Contribution
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.

DISCUSSION
In the current multi-channel MEG study, we demonstrated that our novel weighted-BSS T/k method using only deviant epochs (deviant concatenation) could extract an MMR confined to one or a few dominant components (Figures 4, 6 Table 2). In the subtraction-BSS T/k /weighted-BSS T/k , the salient components showed positive spatio-temporal correlations with the MMR (Figures 7,  8, and Supplementary Figure 3). However, ICA decomposed the MMR into an assembly of minor or pseudo components with negative spatio-temporal correlations. Specifically, our method avoids having to use the conventional subtraction approach to reveal the MMR. Our method may help with the use of the MMR in basic and clinical research.

The Conventional Subtraction Approach to Reveal the MMR
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.

Periodical Arrangements and Weight Assignments
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  Table 1 for each value. ** < 0.005. 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)(36)(37)(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).

Significance of our Approach
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 spatiotemporal 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.

Comparison With ICA
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)(61)(62)(63)(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.

Future Perspectives
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.

LIMITATIONS
There are several methodological concerns to our study: (i) The spatio-temporal cluster permutation provided several clusters ( Figure 3A); however, we did not select all of these. We selected the most reliable clusters that covered 100-250 ms and the bitemporal sensors (7,17) since the vast majority of EEG studies of MMR generators confirmed these; however, the parietal generator in the later latency (e.g., Figure 3A Clusters #5 and #6) was suggested in several studies (9,10) and should be investigated in a future study. (ii) The SOA of the current study was relatively short so that the brain response could return to the baseline. This short time range may have concatenation artifacts when deviant concatenation. However, in the weighted-BSS T/k method, the amplitude outside of the window was 0.2 (Eq. 9). Therefore, concatenation artifacts, if any, should be limited. (iii) The window function was set as a rectangular window, which may cause a tingling effect. The selection of a window function should be based on a hypothesis; in the current study, we assumed that the crucial time range of MMR is equally distributed at 96-276 ms based on our data-driven approach, even though this time range is not assumed to have a unique significance. However, the non-rectangular window can be used according to the hypothesis. Therefore, we uploaded the source code of weighted-BSS T/k to GitHub (https://github. com/fractionalTypeBSS/BSSTk.git) to enable users to apply it according to their hypothesis and select so that users can use it based on their hypothesis to choose the window function and time range. (iv) The sample size was relatively small for fully describing the performance of our new approach. However, generalization, as well as the validity of our approach, is supported by our additional analysis in a separate cohort [Supplementary Data (2-2) and Supplementary Figure 6].

CONCLUSIONS
We proposed a novel weighted method for extracting the MMR from multi-channel MEG data. Compared with ICA, our weighted-BSS T/k method was more sensitive in highlighting the MMR in one or a few dominant components with positive spatio-temporal correlations. This new approach which used only deviant epochs could replace or complement the conventional subtraction approach. Our method may facilitate the use of the MMR in basic and clinical research and provide a novel approach to analyze complex event-related MEG and EEG data.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of Kyushu University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
TMat, SS, JA, and KK: study conception and design. TMat, TMae, and ST: data collection. TMat and KK: analysis and interpretation of results. TMat: draft manuscript preparation. YG, SK, and MH: revision of manuscript. All authors approved the final version of the manuscript.