<figure><img src="image_1.png" /></figure>

HAL
open science

# Maximum-likelihood estimation of channel-dependent trial-to-trial variability of auditory evoked brain responses in MEG

Cezary Sieluzycki, Pawel Kordowski

► To cite this version:

Cezary Sieluzycki, Pawel Kordowski. Maximum-likelihood estimation of channel-dependent trial-to-trial variability of auditory evoked brain responses in MEG. BioMedical Engineering OnLine, 2014, 13, pp.75. 10.1186/1475-925X-13-75. hal-01308207

HAL Id: hal-01308207

https://hal.sorbonne-universite.fr/hal-01308207

Submitted on 27 Apr 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Maximum-likelihood estimation of
channel-dependent trial-to-trial variability of
auditory evoked brain responses in MEG

Cezary Sielużycki<fnref n="1" />,<fnref n="2" />* and Paweł Kordowski³

Abstract

**Background:** We propose a mathematical model for multichannel assessment of the trial-to-trial variability of auditory evoked brain responses in magnetoencephalography (MEG).

**Methods:** Following the work of de Munck *et al.*, our approach is based on the maximum likelihood estimation and involves an approximation of the spatio-temporal covariance of the contaminating background noise by means of the Kronecker product of its spatial and temporal covariance matrices. Extending the work of de Munck *et al.*, where the trial-to-trial variability of the responses was considered identical to all channels, we evaluate it for each individual channel.

**Results:** Simulations with two equivalent current dipoles (ECDs) with different trial-to-trial variability, one seeded in each of the auditory cortices, were used to study the applicability of the proposed methodology on the sensor level and revealed spatial selectivity of the trial-to-trial estimates. In addition, we simulated a scenario with neighboring ECDs, to show limitations of the method. We also present an illustrative example of the application of this methodology to real MEG data taken from an auditory experimental paradigm, where we found hemispheric lateralization of the habituation effect to multiple stimulus presentation.

**Conclusions:** The proposed algorithm is capable of reconstructing lateralization effects of the trial-to-trial variability of evoked responses, i.e. when an ECD of only one hemisphere habituates, whereas the activity of the other hemisphere is not subject to habituation. Hence, it may be a useful tool in paradigms that assume lateralization effects, like, e.g., those involving language processing.

**Keywords:** Evoked responses, Habituation, Lateralization, Kronecker product, Maximum-likelihood estimation, MEG, Noise covariance, Trial-to-trial variability

Background

Measurements of evoked brain responses in humans are an established standard in magneto- (MEG) and electroencephalography (EEG). The measured quantity is either the magnetic field above the surface of the head in MEG [1] or the electric field on the surface in EEG [2]. The sources of activity that underlie the measured signal are electric currents generated in a certain area of the cortex. Unfortunately, the background

<footnote>*Correspondence:</footnote>

<footnote>cezary.sieluzycki@icm-institute.org</footnote>

<footnote n="1">¹Special Lab Non-invasive Brain
Imaging, Leibniz Institute for
Neurobiology, Brenneckestr. 6,
39118 Magdeburg, Germany</footnote>

<footnote n="2">²Team Normal and Abnormal
Motor Control, ICM Brain and Spine
Institute, Sorbonne University,
Pierre-and-Marie-Curie University
(Paris 6), INSERM UMR1127, CNRS
UMR7225, Hôpital Pitié Salpetrière,
47 bd de l'Hôpital, 75013 Paris,
France
Full list of author information is
available at the end of the article</footnote>
noise present in MEG and EEG measurements significantly exceeds the amplitude of
the evoked brain activity of interest. The model commonly used to describe these data
is the so-called signal-plus-noise (SPN) model. It assumes that: 1) the brain response
to identical stimuli does not change from trial to trial, i.e. from one stimulation to
another, and 2) the noise contaminating the response has a Gaussian distribution with
zero mean. Hence, simple averaging of the sequentially recorded responses with respect
to stimulus onset is supposed to yield a good estimate of the “real” brain response to the
stimulus.

In the framework of the SPN model, the signal $R_{ij}^{(k)}$ measured on the *i*-th channel, at the *j*-th time moment, in the *k*-th trial, can be described as

$$
R_{ij}^{(k)} = R_{ij} + \epsilon_{ij}^{(k)}, \qquad (1)
$$

where $i = 1, \dots, I$, $j = 1, \dots, J$, $k = 1, \dots, K$. $R_{ij}$ is the brain response to the stimulus,
which is assumed to be constant from trial to trial, and $\varepsilon$ is a Gaussian noise with zero
expected value, whose amplitude is typically larger than the amplitude of the evoked brain
response. Due to the assumption of the Gaussian, zero-mean nature of the noise, $\varepsilon$ is
supposed to vanish on averaging across a large number of $R^{(k)}$s.

Averaging across trials is useful when one is interested in the estimation of the evoked
response, but in this approach information on the response variability across single tri-
als is lost. In fact, the trial-to-trial variability is well documented by physiological studies.
The analysis of single evoked responses has been addressed in various ways. For example,
Bartnik *et al.* extracted and parametrized single evoked potentials in EEG by means of the
wavelet transform [3]. Effern *et al.* [4], Quian Quiroga & van Luijtelaar [5], and Bagshaw &
Warbrick [6] tried to increase the signal-to-noise ratio of single trials by wavelet denois-
ing. Various statistical approaches—e.g., [7-9]—including new dedicated tests [10] were
also proposed. De Munck *et al.* used a maximum-likelihood estimation [11], in which
the spatial and temporal covariance of the contaminating noise was employed. Several
coherence aspects were studied, for example, in [12-14]. Westerkamp & Aunon proposed
a linear multielectrode filter to extract the single-trial responses [15], while Georgiadis
*et al.* used the Kalman recursive filtering to obtain an estimate of a single-trial response by
employing the information from the preceding trials [16]. Matching pursuit in the eval-
uation of the trial-to-trial variability of the auditory responses in MEG was used, e.g.,
in [17-20].

This paper aims at improving the quality of the estimators of auditory evoked brain
responses by taking into account their trial-to-trial variability on the multichannel level.
The current work is an extension of the work of de Munck et al. [11], where the spatio-
temporal response pattern of each trial was weighted by a scalar coefficient assumed
identical for all MEG channels. We propose to enhance that approach by estimating a
weighting coefficient for each individual channel. This enables distinguishing between,
e.g., habituation trends which differ for distinct groups of channels thus reflecting differ-
ent trial-to-trial variabilities of the underlying sources—something which is not possible
in the approach of de Munck et al. [11]. We apply the proposed method to simulated data
as well as an illustrative real MEG data set from an auditory experiment, and discuss its
limitations in the Discussion.
Methods

Assumptions

Based upon the realistic assumption that the trial-to-trial variation of evoked brain
response should be taken into account when estimating the response amplitude, the
following extension of the SPN model was proposed by de Munck *et al*. [11]:

$$
R_{ij}^{(k)} = \alpha^{(k)} R_{ij} + \varepsilon_{ij}^{(k)}, \qquad (2)
$$

where, for each trial k, α⁽ᵏ⁾ is a coefficient scaling the spatio-temporal amplitude pattern
of the response, R, in this trial. In that work, the spatial and temporal covariance of the
contaminating noise was estimated in order to minimize its influence on the estimates of
α and R.

In this work, we propose to replace the scalar α⁽ᵏ⁾ with a diagonal matrix ψ⁽ᵏ⁾ in order
to trace the trial-to-trial amplitude variation for individual channels i = 1, . . . , I :

$$
\psi^{(k)} = \begin{bmatrix} \psi_1^{(k)} & 0 & \dots & 0 \\ 0 & \psi_2^{(k)} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \psi_I^{(k)} \end{bmatrix}. \tag{3}
$$

Therefore, the order of the matrix $\psi^{(k)}$ equals the number of channels $I$. In this way, in every single trial $k$, each of the channels $i$ is assigned to the scalar coefficient corresponding to the evoked response amplitude variations on this channel. Thus, the fundamental equation of the considered process, (2), is replaced by

$$
R_{ij}^{(k)} = \psi_{i}^{(k)} R_{ij} + \varepsilon_{ij}^{(k)}. \qquad (4)
$$

Additionally, the following assumption is made

$$
\sum_k (\psi^{(k)})^2 = \begin{bmatrix} K & 0 & \dots & 0 \\ 0 & K & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & K \end{bmatrix}, \qquad (5)
$$

which, since *K* is the number of trials, says that in the case of non-variability of the
response, each of the components on the leading diagonal of ψ⁽ᵏ⁾ is equal to one.

The model (4) poses further difficulties, because the dimensionality of the estimation
process increases, yet we hope it offers a more accurate reflection of reality. Since it is
known that the brain response evoked by a characteristic sensory stimulus arises in a par-
ticular area of the cortex (see, e.g., [21-24]), the magnetic fields measured by the channels
corresponding to this areaᵃ are stronger than those acquired from less active parts of the
brain. Thus, it is reasonable to distinguish between individual channels when estimating
the variability of the evoked response. For example, in their study of the habituation of the
P300 wave for visual stimuli, Ravden & Polich observed that habituation was strongest at
the Fz and Cz electrode sites [25].

Following the approach of de Munck *et al*. [26], the following additional assumptions concerning the noise space-time correlation are made:

* the spatio-temporal covariance matrix, C, of the noise is the Kronecker product (see, e.g., [27]), ‘⊗’ of the spatial, X, and the temporal, T, covariance matrices,

$$
C = X \otimes T, \tag{6}
$$
i.e., for $X_{ii'}$ denoting the scalar being the spatial covariance between channel $i$ and
channel $i'$, and $T_{jj'}$ being the temporal covariance between sample $j$ and sample $j'$,
$C$ is of the following form

$$
\begin{bmatrix}
X_{11} T_{11} \dots X_{11} T_{1J} & \dots & X_{1I} T_{11} \dots X_{1I} T_{1J} \\
\vdots & \ddots & \vdots \\
X_{11} T_{J1} \dots X_{11} T_{JJ} & \dots & X_{1I} T_{J1} \dots X_{1I} T_{JJ}
\end{bmatrix}
,
\qquad (7)
$$

• the noise of different trials is uncorrelated.

Although the dimensions of $C = X \otimes T$ are not smaller than the dimensions of the
spatio-temporal covariance of the noise built traditionally, assumption (6) permits split-
ting $C$ into $X$ and $T$. Since the noise covariance matrix has to be estimated and inverted,
this splitting is a clear advantage thanks to the fact that $(X \otimes T)^{-1} = X^{-1} \otimes T^{-1}$.

It is clear from (6) that a common pattern of both *X* and *T* across trials *k* is assumed,
as well as that *X* is time invariant and that *T* is space invariant, which of course is a
simplification (see [28]), yet, a computationally efficient one.

Estimation of evoked response varying across trials

To derive the maximum-likelihood estimators of X, T, ψ⁽ᵏ⁾ and R, one needs to find the distribution of εᵢⱼ⁽ᵏ⁾. Due to the assumption of the zero-mean, Gaussian nature of the noise and the aforementioned assumptions concerning the noise correlation, the covariance between two noise realizations (one on channel i, at sample j, in trial k, and the other on channel i', at sample j', in trial k') equals [26]

$$
E \left( \varepsilon_{ij}^{(k)} \varepsilon_{i'j'}^{(k')} \right) = X_{ii'} T_{jj'} \delta_{kk'}, \quad (8)
$$

where $\delta_{kk'}$ is the Kronecker delta defined as

$$
\delta_{kk'} = \begin{cases} 1 & \text{for } k = k' \\ 0 & \text{for } k \neq k'. \end{cases} \tag{9}
$$

Thus, the determinant of the overall covariance equals [27]

$$
\det \left( E \left( \varepsilon_{ij}^{(k)} \varepsilon_{i'j'}^{(k')} \right) \right) = (\det X)^{JK} (\det T)^{IK} . \quad (10)
$$

Moreover, since

$$
\sum_{i,j,k,i',j',k'} (R_{ij}^{(k)} - \psi_i^{(k)} R_{ij}) (X^{-1})_{ii'} (T^{-1})_{j'j} \delta_{kk'} (R_{i'j'}^{(k')} - \psi_{i'}^{(k')} R_{i'j'}) = \\
= \sum_k \operatorname{tr} \left( (R^{(k)} - \psi^{(k)} R) T^{-1} (R^{(k)} - \psi^{(k)} R)^\mathrm{T} X^{-1} \right), \quad (11)
$$

the probability density distribution $f_{\epsilon}$ of $\epsilon$ is expressed by

$$
f_{\epsilon}(R, \psi^{(k)}, X, T) = \exp \left( -\frac{1}{2} \sum_{k} \operatorname{tr} \left( (R^{(k)} - \psi^{(k)} R) T^{-1} (R^{(k)} - \psi^{(k)} R)^{\mathrm{T}} X^{-1} \right) \right). \quad (12) \\
\cdot (2\pi)^{-JK/2} (\det T)^{-IK/2} (\det X)^{-JK/2}.
$$
In the context of the maximum-likelihood estimation of X, T, ψ(k) and R, (12) represents the likelihood function. For the sake of keeping the formulas simple, all estimators are denoted as non-dashed (i.e. by X instead of X̂, for example).

### Calculation of the estimates of X, T, ψ and R

Differentiating (12) with respect to X, T, ψ(k) and R, and setting the obtained derivatives to zero, one obtains the following estimators (see Appendix for detailed derivations).

$$X = \frac{1}{IK} \sum_k (R^{(k)} - \psi^{(k)} R) T^{-1} (R^{(k)} - \psi^{(k)} R)^T, \quad (13)$$

$$T = \frac{1}{IK} \sum_k (R^{(k)} - \psi^{(k)} R)^T X^{-1} (R^{(k)} - \psi^{(k)} R), \quad (14)$$

$$\psi^{(k)} = \mathrm{dg} \left( (X^{-1} \odot (RT^{-1}R^T))^{-1} \mathrm{diag} \left( X^{-1}R^{(k)}T^{-1}R^T \right) \right), \quad (15)$$

$$R = \left( \sum_k \psi^{(k)} X^{-1} \psi^{(k)} \right)^{-1} \sum_k \psi^{(k)} X^{-1} R^{(k)}, \quad (16)$$

where ⊙ denotes the Hadamard product, i.e. $(Y \odot Z)_{mn} = Y_{mn}Z_{mn}$, diag denotes the operator resulting in a vector that contains elements of the leading diagonal of the matrix that diag acts on, and dg denotes the operator resulting in a diagonal matrix whose leading diagonal contains elements of the leading diagonal of the matrix (or elements of the vector, as in this case) that dg acts on and zeros elsewhere.

The above system of matrix equations (13)–(16) has no closed form solution. Hence, in order to find the approximated estimates of X, T, ψ and R, we propose to solve it iteratively, starting with X, T and ψ(k) being identity matrices, and with $R = K^{-1} \sum_k R^{(k)}$. In such an approach, consecutive updates of X, T, ψ and R are computed based on their values from the previous iterations in which formulae (13)–(16) are computed consecutively. At the end of each iteration, each ψᵢ is scaled according to (5).

### Material

Real MEG data were acquired from a healthy male volunteer (age 42, normal hearing), who is a member of the Special Lab Non-invasive Brain Imaging. The subject was exposed binaurally to 500-ms long 1-kHz sine tones at 90 dB SPL. The difference between the onsets of the consecutive stimuli was 2 s. The data were sampled at the frequency of 1017.25 Hz and filtered between 0.5 and 400 Hz. 160 trials were extracted out of 224 as artifact free. They were offset corrected for each channel with respect to the mean magnetic field over the 200-ms long prestimulus baseline time window.

### Results

#### Simulations

In order to test to which extent the proposed methodology correctly assesses the channel-dependent trial-to-trial variability of the spatio-temporal response pattern R, we simulated two equivalent current dipoles (ECDs), one in each auditory cortex, whose time courses were habituating differently over trials. On the sensor level, the simulated signal can then be expressed as

$$R_{ij}^{(k)} = L_1(i) h_1(k) s_1(j) + L_2(i) h_2(k) s_2(j) + \varepsilon_{ij}^{(k)}, \quad (17)$$
where $L_1(i)$ denotes the leadfield of source 1 to sensor $i$, $h_1(k)$ represents the trial-to-trial variation of the signal $s_1(j)$ of source 1, and likewise $L_2(i)$, $h_2(k)$, and $s_2(j)$. Of course, $L_1(i) h_1(k) s_1(j) + L_2(i) h_2(k) s_2(j)$ is not identical to $\psi_i^{(k)} R_{ij}$ in (4), thus (4) is an approximation rather than the exact representation. Therefore, the method *formally* applies, i.e. it is able to perfectly disentangle $h_1$ from $h_2$, only when the two leadfields are mutually orthogonal, i.e. when $L_1(i)L_2(i) = 0$ for all $i$; then, the Hadamard product of the two related forward fields $L_1(i) h_1(k) s_1(j)$ and $L_2(i) h_2(k) s_2(j)$ is a zero field. However, since the two simulated sources were clearly separated in space, we assumed that this was practically the case.

Based on a magnetic-resonance-imaging (MRI) scan of our subject (see Section *Material*), we selected one location per auditory cortex where the M100/N100 response (see, e.g., [29-31], and references therein) to simple 1-kHz tones is supposed to be generated. Assuming the spherical head model, at each of these locations we seeded a tangential ECD of unit moment (figure 1), whose amplitude time course was simulated as the $[x, y, z]$-moment multiplied by a normalized time series generated as $s_1(j) = \sin(3j/J)$ for the left hemisphere and $s_2(j) = \sin(3j/J + \pi/25)$ for the right hemisphere, to simulate the M100 wave at its peak proximity and to account for a possible phase (or peak

<figure><img src="image_3.png" /></figure>

**Figure 1** Simulated ECD locations. Two tangential ECDs of unit length were seeded (at the back ends of the red arrows) in the auditory cortices. The corresponding forward field is shown in the bottom-right panel. The spherical head model was applied, and the green sphere was fitted using FieldTrip [32], based on 192 surface points selected from head digitization. Only the points close to the two auditory cortices were used for the fit. The blue lines intersect at the anterior commissure.
latency) difference between the time courses of the two sources; see figure 2. In the sim-
ulation, *J* = 51, which covered 50 ms at the sampling frequency of 1017.25 Hz. These
trial-invariant time courses were next multiplied across trials by scalar coefficients orig-
inating from two linear trends (one for each of the two ECDs) with negative slopes, to
model the habituation of the two sources across trials. Namely, each of the two linear
trends was a vector containing *K* weighting coefficients. Those vectors were simulated as
*h*₁(*k*) = −*k* + 2*K* for *s*₁ and *h*₂(*k*) = −*k* + 3*K*/2 for *s*₂, and scaled according to (5) each,
resulting, for *K* = 160, in the slopes: −0.0041 for the left and −0.0060 for the right source.

The forward field resulting from the two sources was calculated for each trial *k* and contaminated by noise taken from the −50–0-ms prestimulus baseline period of individual trials of a real MEG measurement of the same subject, whose MRI was used to seed the simulated ECDs. Both the forward-field signal and the noise were normalized with respect to their norms across all channels *i*, samples *j* and trials *k*, and then the noise was scaled such that the overall signal-to-noise (SNR) ratio equaled 0.1, a very demanding but realistic value, which we based on an estimate of SNR from a real data set.

Figure 3 shows the estimates obtained from the iterative solution of equations (13)–(16) after 20 iterations. The top panel shows the estimated trial-invariant pattern of the response, *R*, at *j* = 25, i.e. between the two M100 peaks originating from the two ECDs. A strong correspondence to the forward field of the simulated ECDs (see figure 1) can be seen. Since the estimated *ψ*, which scales *R* across trials (see (4)), is typically noisy, 1st-degree polynomials were fitted (see [11,17,18]) to the *ψ* estimates on each channel in order to reveal the main tendency in the overall trial-to-trial variability. Figure 3 (middle) shows the slopes of the fitted polynomials for each channel. We only show the values for those sensors at which the value of the slopes was significantly different from zero at the Bonferroni-corrected significance level of 0.05/245, where 245 is the number of artifact-free channels (we excluded 3 sensors from the 248-channel set of our 4D MEG system due to strong artifacts on those channels; let us remind here that we used a real-measurement noise in these simulations). Sensors marked with empty circles denote the slopes which were not significantly different from zero. Note that those are predominantly channels with poor SNR, as revealed by the top panel. The significant slopes reveal the difference in habituation between the two auditory cortices, with the mean equal to −0.0041 for

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

**Figure 3 Results of simulations.** Top: the trial-invariant spatial response pattern *R* at *j* = 25, i.e. for the time point at the middle of the analyzed time window (see also figure 2) in a 2D projection of the MEG sensors layout (open black circles), for the simulated data generated as the superposition of the habituating activity (see figure 2) of the ECDs depicted in figure 1 and the contaminating spatially and temporally correlated background noise, at the overall SNR = 0.1. Middle: Spatial distribution of the slope coefficients of 1-st degree polynomials fitted to the multichannel estimates of ψ reflecting the trial-to-trial variability of *R*. Sensors marked with blue-filled circles denote the negative slopes significantly different from zero at the 0.05/245 significance level, where 245 is the number of analyzed channels, revealing a clear habituation. The exact color of a circle reflects the value of the slope. Bottom: The estimate of the temporal covariance *T* of the noise.
the left-hemisphere sensors and −0.0056 for the right hemisphere, compared to −0.0041
and −0.0060 at the source level. In both cases, the corresponding standard deviation σ
was equal to 0.0007. Hence, the habituation lateralization on the source level was well
reconstructed on the sensor level despite the demanding overall SNR of 0.1.

The bottom panel of figure 3 shows the estimate of the temporal covariance, *T*, of the contaminating noise, which reveals some nonstationarity [33,34]. We do not show the somewhat peculiar pattern of the estimate of *X*, since it stems from the particular spiral arrangement of the sensors labels in the 4D MEG system, which makes visual inspection of *X* impaired.

For illustrative purposes, figure 4 shows, for a selected channel revealing a strong signal
and a clear habituation, the estimated ψ(k) (black dots) and the consecutively fitted 1st-
degree polynomial.

In addition, we simulated a scenario in which only one of the two sources, namely
the one in the right hemisphere, habituated. For that source, habituation strength was
the same as before, whereas for the left-hemisphere EDC, $h_1(k) = 1$ independent
of $k$. Figure 5 shows results of the estimation, in full analogy to those from figure 3.
Note that the algorithm managed to reconstruct the habituation on the sensors that
were located above the right hemisphere, while those for the left hemisphere remained
essentially unaffected. This time the mean of the significantly non-zero slopes of the 1st-
degree polynomials fitted to the estimates of $\psi$ for the right-hemisphere sensors was
equal to -0.0066 (with $\sigma = 0.0012$), still close to the slope of the simulated $h_2$, i.e.
to -0.0060.

In order to test the performance of the method in situations in which the leadfields of
the ECDs are far from mutual orthogonality, we simulated an additional scenario in which
a third ECD was placed in the right hemisphere; see figures 6 and 7. The time course of
this third ECD was simulated as $s_3(j) = \sin(3j/5 + \pi/15)$ normalized with respect to the
peak value. Due to two ECDs in the right hemisphere in this scenario, the overall ECD
signal strength in that hemisphere was larger, compared to that in the left hemisphere.
Moreover, due to the aforementioned scaling of the signal and the noise resulting in an
overall SNR of 0.1, the SNR for the left-hemisphere sensors was now weaker than in the

<figure><img src="image_6.png" /><figcaption>**Figure 4 Habituation.** 1st-degree polynomial (blue line) fitted to $\psi^{(k)}$ for a right-hemisphere channel revealing a strong signal and a clear habituation (see also figure 3). For this channel, the fitted line $p(k) = -0.0058k + 1.2681$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5 Results of simulations.** Like in figure 3 but for another simulation in which the left-hemisphere ECD did not habituate. Note the spatial selectivity of the habituation estimate on the sensor level (the middle figure).</figcaption></figure>

previous simulations. To examine the influence of the trial-to-trial characteristics of the third source on the overall habituation pattern for the right hemisphere, we assumed that this additional source did not habituate, i.e. that its $h_3(k) = 1$ for all $k$ ($h_1$ and $h_2$ were like in the first simulation; see figure 2). Since, as aforementioned, ECDs whose leadfields are not distinct enough, i.e. not orthogonal or almost orthogonal, can not be efficiently separated, the estimated overall habituation pattern for the right hemisphere is weakened in this scenario; see figure 8. This is due to the fact that the third source did not habituate. The observed pattern is a mixture of the trial-to-trial characteristics of the two sources
<figure><img src="image_8.png" /><figcaption>**Figure 6 Simulated ECD locations.** Like in figure 1 but with the additional (third) ECD in the right auditory cortex.</figcaption></figure>

in the right hemisphere; the mean of the significantly non-zero slopes of the 1st-degree
polynomials fitted to the estimates of $\psi$ for sensors above the right-hemisphere was equal
to -0.0036 (with $\sigma = 0.0006$). This will of course be the case for any real-data sce-
nario with multiple sources located close to each other. However, since the exact number
of ECDs is not known in practice, we do not consider such an “integrative”, and hence
impaired, resolution of the method a serious flaw. The fact that we do not know, and hence
do not want to assume explicitly, the exact number of ECDs for real data is essentially the
reason for estimating the trial-to-trial variability on the sensor rather than source level.

Since, as mentioned above, the SNR for the left-hemisphere sensors in this scenario
was smaller than in the two previous simulations, the habituation pattern for the ECD
of the left auditory cortex was now revealed less clearly. Compared to the result shown
in figure 3, the habituation pattern of the left hemisphere is now weaker (although
some sensors still reveal significant decay, with the mean significant slope -0.0049,
$\sigma = 0.0003$).

Real MEG data

Figure 9 shows the findings for real MEG data (see Section *Material*) after 20 iterations
(see Section *Calculation of the estimates ...*). The analyzed time window was 70–120 ms,
<figure><img src="image_9.png" /><figcaption>**Figure 7 Simulated ECD time courses.** As in figure 2 but for the three sources shown in figure 6, with the green color corresponding to the additional ECD in the right hemisphere.</figcaption></figure>

which covered the positive/negative part of the M100 waveform for the left/right channel
with the best SNR.

The top panel of figure 9 shows the estimated trial invariant response *R* at *j* = 25 cor-
responding to *t* = 95 ms. This pattern looks very close to the average signal at this time
point. In the middle panel, a clear habituation pattern is visible for several sensors located
above the right auditory cortex, as they reveal significantly negative slopes at the 0.05/245
significance level. The mean of those significant slopes is −0.0041 (σ = 0.0005). In this
subject, the trial-averaged signal reveals an overall worse SNR above the left, compared
to the right hemisphere, which might be the reason for this particular lateralized result.

The bottom panel shows the estimated temporal covariance matrix *T* of the noise,
which reveals some tiny traces of nonstationarity.

Figure 10 for the real-data analysis is fully analogous to figure 4 for the simulation. Alter-
natively to the 1st-degree polynomial, one might fit an exponential function. However,
given the observed level of the inter-trial variance of ψ and the fact that in the real-data
scenario the exact characteristics of habituation is not known, we decided to follow the
Ockham’s razor principle and hence opt for linear fits. Besides, exponential fits resulted
in exactly the same channels revealing significant habituation.

**Discussion**

We have introduced an extension of the maximum-likelihood approach for the estima-
tion of the trial-to-trial variations of evoked brain responses in MEG based on the work
of de Munck *et al*. [11] to obtain sensor-specific estimates. In [11], the trial-invariant
spatio-temporal response pattern was weighted in each trial by a scalar coefficient which
was identical for all sensors. Here, we further developed this approach by estimating a
coefficient for each individual sensor. We checked its applicability using simulated data
and found that the lateralization of the habituation simulated on the source level could
be reconstructed on the sensor level, even for a very demanding SNR. The analysis of
illustrative real MEG data revealed an auditory-cortex related habituation pattern, too,
showing that the method is capable of finding significant habituation in data from real
experiments.
<figure><img src="image_10.png" /><figcaption>**Figure 8 Results of simulations.** Like in figures 3 and 5 but for the scenario depicted in figures 6 and 7.</figcaption></figure>

The estimation of the channel-dependent $\psi$ is computationally much more demanding compared to the scenario with $\alpha$ assumed identical to all channels [11]. This is so because the dimensionality of the problem increases, while the probability space becomes reduced by one dimension—the number of channels. This makes the overall estimation difficult and possibly less robust compared to the case with the scalar $\alpha^{(k)}$. Assuming that all the cortical processes were characterized by the same trial-to-trial variability pattern, the use of the simpler approach of de Munck *et al*. [11] would be more appropriate. However, it is not realistic to assume a priori that all processes involved will habituate identically, especially not in scenarios with complex stimulation paradigms. In fact, some may even
<figure><img src="image_11.png" /></figure>

**Figure 9 Results for the real MEG data.** Top: the estimate of *R* [fT] at *j* = 25 corresponding to *t* = 95 ms. Middle: the slopes of 1st-degree polynomials fitted to the estimates of *ψ*. Sensors marked with blue-filled circles reveal a clear habituation at the restrictive 0.05/245 significance level. Bottom: the estimate of the temporal covariance *T* of the noise.

reveal the opposite behavior, i.e. sensitization (see, e.g., [35]). Even for simple auditory stimuli, the habituation pattern may reveal differences between the two hemispheres—something that could not be addressed with the model of de Munck *et al*. [11]. Hence, in such cases, topographical representations of the trial-to-trial variability will be more appropriate.

An important issue to be taken into consideration when applying the advanced approach presented here is the transition from the source-level model to the sensor-level model. If all sources were characterized by exactly the same trial-to-trial variability, the
<figure><img src="image_12.png" /><figcaption>**Figure 10 Habituation in real-data analysis.** The blue line depicts the 1st-degree polynomial $p(k) = -0.0042k + 1.2079$ fitted to $\psi^{(k)}$ for a right-hemisphere channel revealing a strong signal and a clear habituation (see figure 9).</figcaption></figure>

model with scalar $\alpha^{(k)}$ [11] would be sufficient, as it corresponds to linear scaling of the forward equation in each trial. Such simple scaling is, however, not possible if each source has its own trial-to-trial variation, as assumed in the simulations (see Section *Simulations*). Hence, our model, given by (4), is not exact. We are aware of this inaccuracy, yet, for sources that are sufficiently separated in space, as it is the case when a single ECD per auditory cortex is assumed, this inaccuracy is rather small. In this case, one can treat the problem almost as separate forward problems, one for each source. Of course, the closer in space the sources, the more problematic this issue becomes, which is why the applicability of the proposed approach in scenarios with differently habituating sources in close proximity to each other would be questionable. Alternatively, one could propose a linear expansion into several $\alpha^{(k)}R$ products with scalar $\alpha^{(k)}$'s, where each term would correspond to the forward field of a single underlying source. However, this approach would require knowledge of the exact number of sources a priori, something which is usually not available. Therefore, we consider our model to be a reasonable approximation in paradigms where one can assume spatially distinct sources.

We critically examined the convergence of the algorithm and found that the estimates of $\psi$ and $R$ show some tiny oscillations with consecutive iterations, i.e. that the ratio of the norms of the estimates in two consecutive iterations does not decay monotonically. Nonetheless, in numerous simulations as well as for real data, we found that 20 iterations is a number which, however somewhat arbitrary, results in the spatial distribution of the estimate of $R$ resembling that of the average signal very well and in $\psi$ nicely revealing the habituation patterns of the underlying sources. A good reason for not proceeding further is—apart from the related time cost—the fact that for real data the algorithm might result in patterns of $R$ which would differ from the trial-averaged signal a bit too strongly. This may be caused by some information flow (leakage) between the estimates of $\psi$ and $R$. Hence, we infer that about 20 iterations offer a reasonable trade-off between the quality of the estimate of $\psi$ and that of the $R$ estimate.

The moderate nonstationarity in the estimated temporal covariance of the noise, depicted in figures 3, 5, 8, and 9, may stem from the intrinsic properties of the noise. A recent work by Kipiński *et al.* [34], which examined nonstationarity of MEG data using
modern statistical tests adopted from econometrics, revealed variance-nonstationarity of the analyzed signals. However, the nonstationarity of *T* that we observed may also stem from a non-perfect separation of the stochastic noise from the quasi-deterministic evoked response. These considerations correspond to the ongoing debate on the mechanism of the generation of evoked responses and the associated relation of the phase of spontaneous rhythms with respect to the stimulus onset (see, e.g., [36,37], and references therein).

## Endnote

ªWe do not discuss here the difference between magnetometers or axial and planar gradiometers in this respect. An interested reader is referred to the famous paper of Hämäläinen *et al.* [1].

## Appendix: Detailed derivations

### Estimation of *X* and *T*

Since the likelihood function $f_{\epsilon}$ attains maximum for the same arguments as its natural logarithm, we can differentiate $\ln(f_{\epsilon})$. We will calculate the derivative of the natural logarithm of (12) with respect to *X*, and set it to zero. Thus, we will differentiate

$$ \frac{\partial \ln f_{\epsilon}}{\partial X} = -\frac{\partial}{\partial X} \left( \frac{1}{2} \sum_{k} \operatorname{tr} \left( (R^{(k)} - \psi^{(k)} R) T^{-1} (R^{(k)} - \psi^{(k)} R)^{\mathrm{T}} X^{-1} \right) + \ln \left( (2\pi)^{J/K/2} (\det T)^{K/2} (\det X)^{J/K/2} \right) \right). \quad (18) $$

Let us notice that for any symmetric, non-singular matrix *X*, it is true (see [27], p. 180) that

$$ \frac{\partial}{\partial X} \det X = (\det X) X^{-1}. \quad (19) $$

Furthermore, since here *X* is square, it is true (see [27], p. 178) that

$$ \frac{\partial}{\partial X} \operatorname{tr}(AX^{-1}) = -(X^{-1}AX^{-1})^{\mathrm{T}}, \quad (20) $$

where

$$ A = (R^{(k)} - \psi^{(k)} R) T^{-1} (R^{(k)} - \psi^{(k)} R)^{\mathrm{T}}. \quad (21) $$

Therefore, using basic properties of logarithm and the fact that both *X* and *T* are symmetric, we can write

$$ \frac{\partial \ln f_{\epsilon}}{\partial X} = \frac{1}{2} \sum_{k} \left( X^{-1} (R^{(k)} - \psi^{(k)} R) T^{-1} (R^{(k)} - \psi^{(k)} R)^{\mathrm{T}} X^{-1} \right) - \frac{JK}{2} X^{-1}. \quad (22) $$

After setting the expression above to zero, we obtain (13).

Performing very similar derivation, one obtains (14).

### Estimation of $\psi^{(k)}$ and *R*

Using the facts (see Table 4 on p. 178 in [27]) that

$$ \frac{\partial}{\partial M} \operatorname{tr}(AM) = A^{\mathrm{T}} \quad (23) $$

and

$$ \frac{\partial}{\partial M} \operatorname{tr}(MAM^{\mathrm{T}}B) = B^{\mathrm{T}}MA^{\mathrm{T}} + BMA, \quad (24) $$
it is straightforward to obtain

$$
\begin{align}
-\frac{1}{2} \frac{\partial}{\partial \psi^{(k')}} \sum_k \operatorname{tr} \left( (R^{(k)} - \psi^{(k)} R) T^{-1} (R^{(k)} - \psi^{(k)} R)^\mathrm{T} X^{-1} \right) &= \nonumber \\
&= -\operatorname{dg} (X^{-1} R^{(k')} T^{-1} R^\mathrm{T}) + \operatorname{dg} (X^{-1} \psi^{(k')} R T^{-1} R^\mathrm{T}), \tag{25}
\end{align}
$$

where dg is the operator which replaces every square matrix *M* with a new one, consisting
of the leading diagonal of *M* and having zeros elsewhere. Then, the following algebraic
manipulations lead to equation (15). For the first term of the right-hand side (RHS) of
equation (25) we can write

$$
\forall i \left( \mathrm{dg} \left( X^{-1} R^{(k')} T^{-1} R^{\mathrm{T}} \right) \right)_{ii} = \left( \mathrm{diag} \left( X^{-1} R^{(k')} T^{-1} R^{\mathrm{T}} \right) \right)_{i}, \quad (26)
$$

where $i = 1, \dots, I$, $()_{ii}$ denotes the $i$th element of the leading diagonal of an $I \times I$ matrix,
and $()_i$ denotes the $i$th element of an $I \times 1$ column vector.

For the second term of the RHS of equation (25) we can write

$$
\forall i \left( \mathrm{dg} \left( X^{-1} \psi^{(k')} RT^{-1} R^{\mathrm{T}} \right) \right)_{ii} = \sum_{t,l,m} X_{it}^{-1} \psi_t^{(k')} R_{tl} T_{lm}^{-1} (R^{\mathrm{T}})_{mi}. \quad (27)
$$

where $t = 1, \dots, I$, $l = 1, \dots, J$, $m = 1, \dots, J$.

Since the RHS of equation (25) must be set to zero and since it contains diagonal terms
only, we have *I* equations, one for each *i*, in which the RHS of equation (26) is equal to the
RHS of equation (27).

Next, we can rewrite the RHS of equation (27) as

$$
\sum_t X_{it}^{-1} \left( \sum_{l,m} R_{tl} T_{lm}^{-1} (R^T)_{mi} \right) \psi_t^{(k')} = \sum_t \beta_{it} \psi_t^{(k')}, \quad (28)
$$

where $\beta$ is an $I \times I$ matrix whose elements are given by

$$
\beta_{it} = X_{it}^{-1} \sum_{l,m} R_{tl} T_{lm}^{-1} (R^T)_{mi} = X_{it}^{-1} (R T^{-1} R^T)_{ti} . \quad (29)
$$

Now, since $T^{-1}$ is symmetric, so is $RT^{-1}R^T$, and hence we can write

$$
\beta_{it} = X_{it}^{-1} (RT^{-1}R^{\mathrm{T}})_{it} = (X^{-1} \odot (RT^{-1}R^{\mathrm{T}}))_{it} . \quad (30)
$$

Let us then denote the RHS of equation (26) by $\gamma_i$, i.e. the $i$th element of a column vector
$\gamma$. Since the RHS of (28) is assumed to be equal to $\gamma_i$, we can write

$$
\sum_t \beta_{it} \psi_t^{(k')} = \gamma_i \tag{31}
$$

and therefore, in the matrix notation,

$$
\beta \mathrm{diag}(\psi^{(k')}) = \gamma. \tag{32}
$$

Then, of course,

$$
\mathrm{diag}(\psi^{(k')}) = \beta^{-1}\gamma, \tag{33}
$$

and hence

$$
\psi^{(k')} = \mathrm{dg}(\beta^{-1}\gamma), \quad (34)
$$
which, given the RHS of (30) for the elements of $\beta$ and the RHS of (26) for the elements of $\gamma$, yields equation (15).

Very similarly, i.e. using (23) and (24), one easily obtains (16).

**Abbreviations**

ECD: Equivalent current dipole; EEG: Electroencephalography; MRI: Magnetic resonance imaging; MEG: Magnetoencephalography; RHS: Right-hand side [of an equation]; SNR: Signal-to-noise ratio; SPN: Signal-plus-noise [model].

**Competing interests**

The authors declare that they have no competing interests.

**Authors' contributions**

CS proposed the model and the estimation algorithm, derived the estimators of $X$ and $T$, performed most of the analyses, and wrote the paper. PK derived the estimators of $\psi$ and $R$, performed some analyses, and contributed substantially to interpretation of data; he also critically revised the model and the manuscript. Both authors read and approved the final manuscript.

**Author's information**

¹Special Lab Non-invasive Brain Imaging, Leibniz Institute for Neurobiology, Brenneckestr. 6, 39118 Magdeburg, Germany. ²Team Normal and Abnormal Motor Control, ICM Brain and Spine Institute, Sorbonne University, Pierre-and-Marie-Curie University (Paris 6), INSERM UMR1127, CNRS UMR7225, Hôpital Pitié Salpêtrière, 47 bd de l'Hôpital, 75013 Paris, France. C. Sielużycki started working on this project while he was with the Special Lab Non-invasive Brain Imaging, Leibniz Institute for Neurobiology, Magdeburg. He is now with the Team Normal and Abnormal Motor Control, ICM Brain and Spine Institute, Sorbonne University, Pierre-and-Marie-Curie University (Paris 6), INSERM UMR1127, CNRS UMR7225, Paris. ³Department of Biomedical Physics, Institute of Experimental Physics, Faculty of Physics, University of Warsaw, ul. Hoża 69, 00-681 Warszawa, Poland. P. Kordowski is currently a trainee in the Special Lab Non-invasive Brain Imaging, Leibniz Institute for Neurobiology, Magdeburg.

**Acknowledgments**

The authors would like to thank the anonymous reviewers as well as Dr Reinhard König and Prof. Wojciech Kordecki for their critical comments which helped to significantly improve the manuscript. The authors acknowledge support from the German Academic Exchange Service (DAAD) within its project based personnel exchange programme (PPP) between Germany and Poland, 2013–2014. C. Sielużycki was supported by the German Research Foundation and is now supported by the Fondation pour la Recherche Médicale (FRM). P. Kordowski is supported by the Foundation for Polish Science: International PhD Projects Programme co-financed by the EU European Regional Development Fund.

**Author details**

¹Special Lab Non-invasive Brain Imaging, Leibniz Institute for Neurobiology, Brenneckestr. 6, 39118 Magdeburg, Germany. ²Team Normal and Abnormal Motor Control, ICM Brain and Spine Institute, Sorbonne University, Pierre-and-Marie-Curie University (Paris 6), INSERM UMR1127, CNRS UMR7225, Hôpital Pitié Salpêtrière, 47 bd de l'Hôpital, 75013 Paris, France. ³Department of Biomedical Physics, Institute of Experimental Physics, Faculty of Physics, University of Warsaw, ul. Hoża 69, 00-681 Warszawa, Poland.

Received: 29 December 2013 Accepted: 10 April 2014

Published: 17 June 2014

**References**

1. Hämäläinen M, Hari R, Ilmoiemi RJ, Knuutila J, Lounasmaa OV: **Magnetoencephalography—theory, instrumentation, and applications to non-invasive studies of the working human brain.** *Rev Modern Phys* 1993, **65**:413–497.

2. Schomer DL, Lopes da Silva FH: *Niedermeyer's Electroencephalography: Basic Principles, Clinical Applications, and Related Fields*. Philadelphia, PA: Lippincott Williams & Wilkins; sixth edition; 2010.

3. Bartnik EA, Blinowska KJ, Durka PJ: **Single evoked potential reconstruction by means of wavelet transform.** *Biol Cybernet* 1992, **67**:175–181.

4. Effern A, Lehnertz K, Fernández G, Grunwald T, David P, Elger CE: **Single trial analysis of event related potentials: non-linear de-noising with wavelets.** *Clin Neurophys* 2000, **111**:2255–2263.

5. Quiroga Quian R, van Luijtelaar ELJM: **Habituation and sensitization in rat auditory evoked potentials: a single-trial analysis with wavelet denoising.** *Int J Psychophys* 2002, **43**:141–153.

6. Bagshaw AP, Warbrick T: **Single trial variability of EEG and fMRI responses to visual stimuli.** *Neuro Image* 2007, **38**:280–292.

7. Gasser T, Möcks J, Verleger R: **SELAVCO: a method to deal with trial-to-trial variability of evoked potentials.** *Electroencephalograph Clinic Neurophys* 1983, **55**:717–723. [Technical Section].

8. Bechtereva NP, Gogolitsin YL, Pakhomov SV: **Single-trial decomposition: a new approach to the analysis of neuronal reactions during psychological tasks.** *Int J Psychophys* 1986, **3**:275–286.

9. Riera-Díaz JJ, Carballo-González JA, Biscay-Lirio R, Valdés-Sosa PA: **The estimation of event related potentials affected by random shifts and scalings.** *Int J Bio-Med Comput* 1995, **38**(2):109–120.

10. Möcks J, Gasser T, Tuan PD: **Variability of single visual evoked potentials evaluated by two new statistical tests.** *Electroencephalograph Clinic Neurophys* 1984, **57**:571–580.
11. de Munck JC, Bijma F, Gaura P, Sielużycki CA, Branco MI, Heethaar RM: **A maximum-likelihood estimator for trial-to-trial variations in noisy MEG/EEG data sets.** *IEEE Trans Biomed Eng* 2004, **51**(12):2123-2128.

12. Lachaux J, Lutz A, Rudrauf D, Cosmelli D, Le Van Quyen M, Martinerie J, Varela F: **Estimating the time-course of coherence between single-trial brain signals: an introduction to wavelet coherence.** *Neurophysiologie Clinique/Clinic Neurophys* 2002, **32**:157-174.

13. Truccolo W, Ding M, Knuth K, Nakamura R, Bressler S: **Trial-to-trial variability of cortical evoked responses: implications for the analysis of functional connectivity.** *Clinic Neurophys* 2002, **113**:206-226.

14. Mäkinen V, Tiitinen H, May P: **Auditory event-related responses are generated independently of ongoing brain activity.** *Neuro Image* 2005, **24**:961-968.

15. Westerkamp JJ, Aunon JL: **Optimum multielectrode a posteriori estimates of single-response evoked potentials.** *IEEE Trans Biomed Eng* 1987, **34**:13-22.

16. Georgiadis SD, Ranta-aho PO, Tarvainen MP, Karjalainen PA: **Single-trial dynamical estimation of event-related potentials: A Kalman filter-based approach.** *IEEE Trans Biomed Eng* 2005, **52**(8):1397-1406.

17. Sielużycki C, König R, Matysiak A, Kuś R, Ircha D, Durka PJ: **Single-trial evoked brain responses modeled by multivariate matching pursuit.** *IEEE Trans Biomed Eng* 2009, **56**:74-82.

18. Sielużycki C, Kuś R, Matysiak A, Durka PJ, König R: **Multivariate matching pursuit in the analysis of single-trial latency of the auditory M100 acquired with MEG.** *Int J Bioelectromagnetism* 2009, **11**(4):155-160.

19. Bénar CG, Papadopoulo T, Torrésani B, Clerc M: **Consensus Matching Pursuit for multi-trial EEG signals.** *J Neurosci Methods* 2009, **180**:161-170.

20. Jörn M, Sielużycki C, Matysiak MA, Żygierewicz J, Scheich H, Durka PJ, König R: **Single-trial reconstruction of auditory evoked magnetic fields by means of Template Matching Pursuit.** *J Neurosci Methods* 2011, **199**:119-128.

21. Ilmoniemi RJ: **Magnetic source imaging.** In *Biological Effects of Electric and Magnetic Fields*, Volume 2. Edited by Carpenter DO, Ayrapetyan S. San Diego, CA: Academic Press; 1994:49-78.

22. Lütkenhöner B, Steinsträter O: **High-precision neuromagnetic study of the functional organization of the human auditory cortex.** *Audiol Neuro-Otol* 1998, **3**:191-213.

23. Ossenblok P, Wilts G, Numminen J, Peters MJ, Lopes da Silva FH: **Locating the cortical sources of somatosensory evoked responses by integration of EEG and MEG.** *Funct Neurosci EEG Suppl* 1996, **46**:183-191.

24. Pantev C, Hoke M, Lehnertz K, Lütkenhöner B, Anogianakis G, Wittkowski W: **Tonotopic organization of the human auditory cortex revealed by transient auditory evoked magnetic field.** *Electroencephalograph Clinic Neurophysiol* 1988, **69**:160-170.

25. Ravden D, Polich J: **Habituation of P300 from visual stimuli.** *Int J Psychophysiol* 1998, **30**:359-365.

26. de Munck JC, Huizenga HM, Waldorp LJ, Heethaar RA: **Estimating stationary dipoles from MEG/EEG data contaminated with spatially and temporally correlated background noise.** *IEEE Trans Signal Process* 2002, **50**(7):1565-1572.

27. Magnus JR, Neudecker H: **Matrix differential calculus with applications in statistics and econometrics.** New York: Wiley; 1988.

28. Bijma F, de Munck JC, Heethaar RM: **The spatiotemporal MEG covariance matrix modeled as a sum of Kronecker products.** *NeuroImage* 2005, **27**(2):402-415.

29. Näätänen R, Picton T: **The N1 wave of the human electric and magnetic response to sound: a review and analysis of component structure.** *Psychophysiology* 1987, **24**:375-425.

30. Hari R: **The neuromagnetic method in the study of the human auditory cortex.** In *Auditory Evoked Magnetic Fields and Potentials, Volume 6 of Advances in Audiology*. Edited by Grandori F, Hoke M, Romani GL. Basel: Karger; 1990:222-282.

31. Zacharias N, Sielużycki C, König R, Kordecki W, Heil P: **The M100 component of evoked magnetic fields differs by scaling factors: Implications for signal averaging.** *Psychophysiology* 2011, **48**(8):1069-1082.

32. Oostenveld R, Fries P, Maris E, Schoffelen JM: **FieldTrip: Open Source Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data.** *Comput Intell Neurosci* 2011, **2011**. [doi:10.1155/2011/156869].

33. Bijma F, de Munck JC, Huizenga HM, Heethaar RM: **A mathematical approach to the temporal stationarity of background noise in MEG/EEG measurements.** *Neuro Image* 2003, **20**:233-243.

34. Kipiński L, König R, Sielużycki C, Kordecki W: **Application of modern tests for stationarity to single-trial MEG data: Transferring powerful statistical tools from econometrics to neuroscience.** *Biol Cybernet* 2011, **105**(3-4):183-195.

35. Callaway III E: **Habituation of averaged evoked potentials in man.** In *Habituation, Volume 2*. Edited by Peeke HVS, Herz MJ. New York: Academic Press; 1973:153-174.

36. Yeung N, Bogacz R, Holroyd CB, Cohen JD: **Detection of synchronized oscillations in the electroencephalogram: An evaluation of methods.** *Psychophysiology* 2004, **41**:822-832.

37. de Munck JC, Bijma F: **How are evoked responses generated? The need for a unified mathematical framework.** *Clinic Neurophysiol* 2010, **121**:127-129.

doi:10.1186/1475-925X-13-75

Cite this article as: Sielużycki and Kordowski: Maximum-likelihood estimation of channel-dependent trial-to-trial variability of auditory evoked brain responses in MEG. *BioMedical Engineering OnLine* 2014 **13**:75.