Kuniharu Kishida

*Department of Information Science, Faculty of Engineering, Gifu University Yanagido, Gifu, Japan*

## **1. Introduction**

24 Will-be-set-by-IN-TECH

194 Advances in Brain Imaging

[25] Terol-Villalobos, I.R.: Morphological Image Enhancement and Segmentation, in Advances in Imaging and Electron Physics, P. W. Hawkes Editor, Academic Press, (2001)

[26] Mendiola-Santibañez, J.D and Terol-Villalobos, I. R. : Morphological contrast mappings on partition based on the flat zone notion. Computación y Sistemas, 6 (2002b) 25–37. [27] M. R. Herbert, D. A. Ziegler, C. K. Deutsch,L. M. O'Brien, N. Lange, A. Bakardjiev, J. Hodgson, K. T. Adrien, S. Steele, N. Makris, D. Kennedy1, G. J. Harris and V. S. Caviness Jr ,"Dissociations of cerebral cortex, subcortical and cerebral white matter volumes in

[28] Calderón-González P.L, Parra-Rodriguez M.A., Libre-Rodríguez J.J., Gutiérrez J.V. Análisis espectral de la coherencia cerebral en la enfermedad de Alzheimer. Rev Neurol

[29] Dunkin JJ, Leuchter AF, Newton TF, Cook IA. Reduced EEG coherence in dementia: state

[30] Papovik E, Haynes, L.W Noradrenergic but not cholinergic innervation of the embryonic cortical neuroepithelium rescues germinal and postmitotic cells in heterochronic

[31] Chiaki Itami, Fumitaka Kimura, Tomoko Kohono, Masato Mtsuoka, Masumi Ichikawa, Tdaharu Tsumoto et Al. Brain-derived neurotrophic factor-dependen unmasking of "silent" synapses in the developing mouse barrel cortex. Neurosciencie 2003 (22) :

[32] Goshima Y, Ito T, Sasaki Y, Nakamura F. Semaphorins as signals for cell repulsion and

autistic boys", *Brain* , vol. 126, pp. 1182–1192, 2003.

or trait marker? Biol Psychiatry 1994; 35: 870-9

invasion. J Clin investigation 2002 (109): 993-998.

cocultures. Brain. Res 2000 (853): 227-235

207–273.

2004; 38 (5): 422-427

13069-13074.

Magnetoencephalography (MEG) can monitor the activation of a neuronal population with millisecond temporal resolution, and offer new noninvasive information about basic activities of the human brain. MEG is usable for the description of spontaneous brain activity and the detection of timely events such as stimulus evoked fields. The most recent advances in the field of MEG concerning cortical responses to stimulation are issued from development of multichannel recordings. We can study the temporal order of several cortical areas from averaged waveforms of MEG, e.g., auditory, visual, and somatosensory evoked responses. For three decades, dynamical information contained in MEG has been analyzed mainly in terms of averaged waveforms.

However, lots of brain activities are included in original MEG data. For spontaneous fields there are huge activities in cortical areas. Generally, it is difficult to estimate source locations of current dipole data from MEG data correctly at the present stage, if the number of current dipoles becomes large. The estimation of their locations is an underdetermined problem. In the case of spontaneous fields it is hard to obtain intracerebral communication between active cortical areas. Hence, we consider an overdetermined problem in the case of evoked fields, since their locations of active cortical areas are limited. Furthermore, we should find a possibility to obtain dynamical information by signal processing of fluctuations around concatenate average waveforms, which have been discarded without signal processing in the conventional MEG analysis. The fluctuations may be induced magnetic fields to communicate between cortical areas.

In this chapter, we will show this new direction of MEG analysis by way of example of somatosensory evoked field (SEF), and intracerebral communication between active somatosensory cortices can be expressed by system identification of fluctuations. That is, intracerebral communication between the primary somatosensory cortex in the contralateral hemisphere (cSI) and the contralateral secondary somatosensory cortex (cSII) in 2Hz median nerves stimuli will be discussed from correlation functions of SEF fluctuations by the identification method with feedback model (see Section 3.6). In this chapter, intracerebral communication between cSI and cSII will be studied not from averaged waveforms but from fluctuations around concatenate averaged waveforms by following steps:

to extend the first approach for nonlinear cases owing to brain complexity, and by using BSS method the latter approach will be served to solve the nonlinear difficulty by elimination of concatenate averaged waveforms. In complex systems the latter approach with BSS is superior to the former approach, when the total number of data is limited. When a current dipole is usually generated by thousands of cells, we can observe its magnetic field by a superconducting quantum interference device (SQUID). Hence, SQUID data are macroscopic from the viewpoint of neuronal population activity, and the statistical nature of system size expansion can be taken as macroscopic variables of SQUID data. As pointed out in the system size expansion method developed in statistical physics (Kishida et al., 1976; Kubo et al., 1973; Tomita & Tomita, 1974; Van Kampen, 1961), the macroscopic variables have extensive property and the mean value of them obeys the nonlinear evolution equation, but their fluctuations follow the linear Gaussian equation in the normal case. This linearity of normal case is suitable

In this chapter we will introduce how to separate fluctuations from evoked magnetic field data by using the decorrelation method of blind source separation (BSS). Fluctuations of MEG rather than averaged waveforms should be analyzed for the better understanding of brain

In the present chapter we will study MEG data for 2Hz periodical median nerve stimuli, since activities of the secondary somatosensory cortices may be found as the somatosensory evoked field (Hamada et al., 2002; Wikström et al., 1996). If there is a current dipole in a brain region, a dipole pattern generated from it can be shown by the Biot-Savart law in SQUID channels (see Fig. 2 or 5). In the special case where there is the one-to-one correspondence between SEF current dipoles and BSS components, fluctuations of SEF can be separated easily. In general, there are no one-to-one correspondences between a current dipole and a component of BSS. Let the dishonest BSS of Section 3.3 be a BSS processing for original MEG data with deterministic evoked waveforms. The honest BSS of Section 3.4 is defined as a BSS processing after elimination of concatenate evoked waveforms from MEG. It is difficult to select honest BSS components related to SEF from all BSS components with various dipole patterns in topographic field maps without repeated line spectrum in the power spectral density (PSD), since there are no line spectrum in PSD in honest BSS components. There remains an open problem to be solved to locate SEF fluctuations in the cases without the one-to-one correspondence. However, fluctuations of SEF can be determined from the honest BSS components, when components with nonzero waveforms of the dishonest BSS are transformed into the honest BSS components by using a transformation matrix. That is, the transformation matrix connecting two types of BSS components will be introduced for selecting the fluctuations of evoked magnetic fields (see Eq. (11)). Hence, SEF fluctuations of the somatosensory evoked field (see Eq. (13)) will be separated from MEG data by using the T/k type decorrelation method of honest BSS, after elimination of deterministic concatenate evoked waveforms. Finally, we can obtain correlation functions of SEF fluctuations from the honest BSS. Our identification method will work efficiently for SEF fluctuations in Section 4.5 and can show intracerebral communication in terms of impulse responses between cSI and

for the latter approach with BSS.

functional activities as in Section 4.

cSII.

**2.3 Difficulties in blind source separation**

1) Extraction of SEF fluctuations from MEG will be discussed in Section 3.3.

2) Inverse problem will be shown in Section 3.5 or 4.4.

3) Intracerebral communication will be discussed in Section 4.5 by feedback model identification.

#### **2. Difficulties in system identification of MEG and countermeasures**

#### **2.1 Deterministic waveforms**

Somatosensory activities are known as cSI, bilateral secondary somatosensory cortices (SIIs) and posterior parietal cortices in somatosensory stimuli with inter-stimulus intervals of more than one second (Kakigi et al., 2000; Wikström et al., 1996), and under more than 1 Hz periodic median nerve stimulus, the somatosensory activity is mainly observed at cSI (Wikström et al., 1996). In the previous paper (Kishida, 2009a), responses of evoked magnetic fields were studied for 5Hz periodical median nerve stimuli. In 5Hz periodic median nerve stimulus, there are the 5Hz repetitive line spectra of the somatosensory evoked field (SEF). This suggests that averaged waveforms are deterministic. Statistical properties of SEF are summarized in Kishida (2009a) as SEF includes the deterministic part of concatenate averaged waveforms and the random part of fluctuations around them.

From the Wold decomposition theorem (Mathematical Society of Japan, 1987), a stationary process is uniquely decomposed into a non-deterministic part and a deterministic part. We can define correlation functions in a stationary process, and they are described by time lag. They can be specified as two types: One type can be obtained from fluctuations and the other pseudo type is calculated from periodic functions. Under the condition of small averaged waveforms, dynamical activities of current dipoles of primary somatosensory cortexes were studied from MEG by the decorrelation method of the BSS method in Kishida (2009b). The usage of the decorrelation method can make an improvement on the signal-to-noise ratio of components of interest.

#### **2.2 Nonlinearity**

Assuming that the brain dynamics is linear, the brain functional activities can be identified by using multivariate autoregressive model. To obtain dynamical information from MEG and/or Electroencephalography data, multivariate autoregressive model has mainly been used for identification of networks or connections between cortices (Cantero et al., 2009; Florin et al., 2010; Hui et al., 2010). Furthermore, we can estimate dynamical interaction between active regions in the brain, providing their pathway is known. It means that dynamical relationship between active regions can be analyzed approximately by the multivariate autoregressive model analysis (Akaike & Nakagawa, 1988). However, the autoregressive model is all-pole type, and has no zeros. As a kind of truncated power series (Kishida, 1997b), autoregressive models are not suitable to be used for evaluating exact transfer functions between regions of cortex. To avoid this difficulty, a feedback model which has rational transfer functions with zeros and poles (Kishida, 1994; 1996) has been used for determination of impulse responses between regions of cortex.

Indeed the brain system is nonlinear. One way is to extend the identification method for a nonlinear case, while the other way is to find a linear part separated from MEG. It is difficult 2 Will-be-set-by-IN-TECH

3) Intracerebral communication will be discussed in Section 4.5 by feedback model

Somatosensory activities are known as cSI, bilateral secondary somatosensory cortices (SIIs) and posterior parietal cortices in somatosensory stimuli with inter-stimulus intervals of more than one second (Kakigi et al., 2000; Wikström et al., 1996), and under more than 1 Hz periodic median nerve stimulus, the somatosensory activity is mainly observed at cSI (Wikström et al., 1996). In the previous paper (Kishida, 2009a), responses of evoked magnetic fields were studied for 5Hz periodical median nerve stimuli. In 5Hz periodic median nerve stimulus, there are the 5Hz repetitive line spectra of the somatosensory evoked field (SEF). This suggests that averaged waveforms are deterministic. Statistical properties of SEF are summarized in Kishida (2009a) as SEF includes the deterministic part of concatenate averaged waveforms

From the Wold decomposition theorem (Mathematical Society of Japan, 1987), a stationary process is uniquely decomposed into a non-deterministic part and a deterministic part. We can define correlation functions in a stationary process, and they are described by time lag. They can be specified as two types: One type can be obtained from fluctuations and the other pseudo type is calculated from periodic functions. Under the condition of small averaged waveforms, dynamical activities of current dipoles of primary somatosensory cortexes were studied from MEG by the decorrelation method of the BSS method in Kishida (2009b). The usage of the decorrelation method can make an improvement on the signal-to-noise ratio of

Assuming that the brain dynamics is linear, the brain functional activities can be identified by using multivariate autoregressive model. To obtain dynamical information from MEG and/or Electroencephalography data, multivariate autoregressive model has mainly been used for identification of networks or connections between cortices (Cantero et al., 2009; Florin et al., 2010; Hui et al., 2010). Furthermore, we can estimate dynamical interaction between active regions in the brain, providing their pathway is known. It means that dynamical relationship between active regions can be analyzed approximately by the multivariate autoregressive model analysis (Akaike & Nakagawa, 1988). However, the autoregressive model is all-pole type, and has no zeros. As a kind of truncated power series (Kishida, 1997b), autoregressive models are not suitable to be used for evaluating exact transfer functions between regions of cortex. To avoid this difficulty, a feedback model which has rational transfer functions with zeros and poles (Kishida, 1994; 1996) has been used for determination of impulse responses

Indeed the brain system is nonlinear. One way is to extend the identification method for a nonlinear case, while the other way is to find a linear part separated from MEG. It is difficult

1) Extraction of SEF fluctuations from MEG will be discussed in Section 3.3.

**2. Difficulties in system identification of MEG and countermeasures**

2) Inverse problem will be shown in Section 3.5 or 4.4.

and the random part of fluctuations around them.

identification.

**2.1 Deterministic waveforms**

components of interest.

between regions of cortex.

**2.2 Nonlinearity**

to extend the first approach for nonlinear cases owing to brain complexity, and by using BSS method the latter approach will be served to solve the nonlinear difficulty by elimination of concatenate averaged waveforms. In complex systems the latter approach with BSS is superior to the former approach, when the total number of data is limited. When a current dipole is usually generated by thousands of cells, we can observe its magnetic field by a superconducting quantum interference device (SQUID). Hence, SQUID data are macroscopic from the viewpoint of neuronal population activity, and the statistical nature of system size expansion can be taken as macroscopic variables of SQUID data. As pointed out in the system size expansion method developed in statistical physics (Kishida et al., 1976; Kubo et al., 1973; Tomita & Tomita, 1974; Van Kampen, 1961), the macroscopic variables have extensive property and the mean value of them obeys the nonlinear evolution equation, but their fluctuations follow the linear Gaussian equation in the normal case. This linearity of normal case is suitable for the latter approach with BSS.

In this chapter we will introduce how to separate fluctuations from evoked magnetic field data by using the decorrelation method of blind source separation (BSS). Fluctuations of MEG rather than averaged waveforms should be analyzed for the better understanding of brain functional activities as in Section 4.

#### **2.3 Difficulties in blind source separation**

In the present chapter we will study MEG data for 2Hz periodical median nerve stimuli, since activities of the secondary somatosensory cortices may be found as the somatosensory evoked field (Hamada et al., 2002; Wikström et al., 1996). If there is a current dipole in a brain region, a dipole pattern generated from it can be shown by the Biot-Savart law in SQUID channels (see Fig. 2 or 5). In the special case where there is the one-to-one correspondence between SEF current dipoles and BSS components, fluctuations of SEF can be separated easily. In general, there are no one-to-one correspondences between a current dipole and a component of BSS.

Let the dishonest BSS of Section 3.3 be a BSS processing for original MEG data with deterministic evoked waveforms. The honest BSS of Section 3.4 is defined as a BSS processing after elimination of concatenate evoked waveforms from MEG. It is difficult to select honest BSS components related to SEF from all BSS components with various dipole patterns in topographic field maps without repeated line spectrum in the power spectral density (PSD), since there are no line spectrum in PSD in honest BSS components. There remains an open problem to be solved to locate SEF fluctuations in the cases without the one-to-one correspondence. However, fluctuations of SEF can be determined from the honest BSS components, when components with nonzero waveforms of the dishonest BSS are transformed into the honest BSS components by using a transformation matrix. That is, the transformation matrix connecting two types of BSS components will be introduced for selecting the fluctuations of evoked magnetic fields (see Eq. (11)). Hence, SEF fluctuations of the somatosensory evoked field (see Eq. (13)) will be separated from MEG data by using the T/k type decorrelation method of honest BSS, after elimination of deterministic concatenate evoked waveforms. Finally, we can obtain correlation functions of SEF fluctuations from the honest BSS. Our identification method will work efficiently for SEF fluctuations in Section 4.5 and can show intracerebral communication in terms of impulse responses between cSI and cSII.

In the T/k type decorrelation method, the absolute sum of off-diagonal elements of normalized correlation matrices is minimized at times corresponding to 2Hz and its higher harmonic frequencies. The main processing in the decorrelation method is removing the

where the superscript *T* denotes the transposition of matrix and *z*(*n*) with zero mean is

The Jacobi-like algorithm has been used to solve the simultaneous diagonalization problem approximately on k normalized correlation matrices (Cardoso & Souloumiac, 1996). This process is to determine a square matrix *U* in the problem of minimization of the cost function

where (*UCzz*(*τm*)*UT*)*ij* denotes the *ij*-element of matrix *UCzz*(*τm*)*UT*. Finally, we have *A*−1

**x**(*n*) = *A*(30)**s**(30)

(30)**s**

From the Wold decomposition theorem (Mathematical Society of Japan, 1987), a stationary process is uniquely decomposed into **x**(*n*) = **u**(*n*) + *δ***x**(*n*), where *δ***x**(n) is non-deterministic and **u**(*n*) is deterministic. We can define correlation functions in a stationary process, and they are specified into two types. When concatenate averaged waveforms are very small, the decorrelation method and the identification method used in Kishida (2009b) work efficiently, since correlation functions can be described by time lag. Inversely the decorrelation method and the identification method do not work well, when concatenate averaged waveforms are not small. In the remaining part of this paper we will first eliminate concatenation of averaged waveforms from MEG data, and then apply the decorrelation method and the identification method to eliminated MEG data. Such approach was called "honest" in the (9) of discussion

Supposing there is a subject who has no or negligible small fluctuations of SEF, and that he has

the subject in BSS calculation of Eq. (7). Therefore BSS of Eq. (7) may be called as "dishonest",

averaged waveforms of SEF. However, certain fluctuations of SEF happen to be in **s**

) and **s**(30)(*n*)=(**s**

(30) *e* (*n*) + *Aa*

(30)**s**

(30) *s* (*n*)*T* **s**

*z*(*n*)*z*(*n* + *τm*)*T m* = 1, ..., *k*, (5)

2, (6)

(30) *a* (*n*), (8)

(30) *e* (*n*)*T* **s**

(*n*). (7)

(*k*) :=

(30) *a* (*n*)*T*)*T*

(30) *e* (*n*) of

*V*−1**x**(*n*). Here V is the covariance matrix given by *E*{**x**(*n*)**x**(*n*)*T*}.

|(*UCzz*(*τm*)*UT*)*ij*|

off-diagonal elements of the correlation matrices at fractional type time lag *τm*,

*N*−1 ∑ *n*=0

> *k* ∑ *m*=1 ∑ *i*�=*j*

For the classification of MEG we have used the T/k type of BSS with e.g. *k* = 30:

(30) *s* (*n*) + *Ae*

*N*

*J*(*U*) =

Then, MEG is decomposed by the fractional BSS with *k* = 30:

(30)**s**

(30) *Aa* (30)

**x**(*n*) = *As*

**3.4 Honest BSS and transformation matrix**

(30) *Ae*

*Czz*(*τm*) = 1

orthonormalized as √

where *A*(30) = (*As*

corresponding to Eq. (1).

of Kishida (2009b).

*J*,

*U* √ *V*−1.

#### **3. Methods**

#### **3.1 Magnetoencephalography (subjects and stimuli)**

Five healthy subjects participated in this study. After explaining the nature of the study, informed consent was obtained. The experimental procedures were in accordance with the Declaration of Helsinki. Their median nerves were stimulated electrically with a constant voltage, square-wave pulse of 0.5 ms duration delivered at the right wrist. Stimulus frequency was periodical 2Hz or the inter-stimulus interval (ISI) was 500 ms. MEG data were recorded with a 64-channel whole-head MEG system (NeuroSQUID Model 100; CTF Systems Inc.). SQUID was the axial gradiometer type. Details were reported in Kishida (2009a).

#### **3.2 Data analysis**

MEG signals were digitized at 1250 Hz and filtered with a 300 Hz on-line low-pass filter. MEG data during 175 s per one subject were used as a single sweep for data analysis. The number of median nerve stimuli was 350. The sampling time (Δ*t*) is 0.8 ms and let magnetic field of SQUIDs be

$$\mathbf{x}(t) = \mathbf{x}(n\Delta t) =: \mathbf{x}(n)$$

at time *t* or discrete time *n*. In periodical stimului MEG data are classified into 3 groups as

$$\mathbf{x}(n) = \mathbf{x}\_{\mathbf{s}}(n) + \mathbf{x}\_{\mathbf{e}}(n) + \mathbf{x}\_{\mathbf{a}}(n),\tag{1}$$

where **x***s*, **x***e* and **x***a* are spontaneous and evoked magnetic fields and artifact noises. The SEF evoked magnetic field consists of a deterministic part and fluctuations around it:

$$\mathbf{x}\_{\varepsilon}(n) = L\_{\varepsilon} \mathbf{Q}\_{\varepsilon}(n) = \mathbf{u}(n) + L\_{\varepsilon} \mathbf{Q}\_{\varepsilon}^{\natural}(n),\tag{2}$$

where *Le* is the lead field from the evoked current dipoles, **Qe**(*n*), **u**(*n*) := *r*350**w** is defined by concatenating 350 copies of the average waveforms **w**, and **Q***h e* (*n*) are SEF fluctuations in the evoked magnetic field (see Eq. (16)).

#### **3.3 Fractional type of BSS**

The decorrelation method was developed by Molgedey & Schuster (1994), Ziehe et al. (2000), Murata et al. (2001) and briefly summarized in Kishida et al. (2003). However, BSS performance is strongly dependent on the choice of time lag. The temporal decorrelation method of BSS has an open problem in choice of the time lag (Hironaga & Ioannides, 2007; Kishida, 2008; Tang et al., 2005; Ziehe et al., 2000). Here it should be noted that the decorrelation method is effective for MEG data in a stationary process.

Since PSD of MEG has a line spectrum of 2Hz and those of higher harmonic modes in 2Hz periodical median nerve stimuli, improvements with the fractional type time lag *τm* defined by

$$
\tau\_m = \left(\frac{1250}{2}\right) / m = \frac{625}{m}, \ m = 1, 2, \dots, k \tag{3}
$$

could be made on BSS (Kishida, 2008; 2009a). Hence, we can determine a mixing matrix *A*(*k*) by the blind source separation with the T/k type (fractional) time lag:

$$\mathbf{x}(n) = A\_{(k)} \mathbf{s}^{(k)}(n). \tag{4}$$

4 Will-be-set-by-IN-TECH

Five healthy subjects participated in this study. After explaining the nature of the study, informed consent was obtained. The experimental procedures were in accordance with the Declaration of Helsinki. Their median nerves were stimulated electrically with a constant voltage, square-wave pulse of 0.5 ms duration delivered at the right wrist. Stimulus frequency was periodical 2Hz or the inter-stimulus interval (ISI) was 500 ms. MEG data were recorded with a 64-channel whole-head MEG system (NeuroSQUID Model 100; CTF Systems Inc.).

MEG signals were digitized at 1250 Hz and filtered with a 300 Hz on-line low-pass filter. MEG data during 175 s per one subject were used as a single sweep for data analysis. The number of median nerve stimuli was 350. The sampling time (Δ*t*) is 0.8 ms and let magnetic field of

**x**(*t*) = **x**(*n*Δ*t*) =: **x**(*n*) at time *t* or discrete time *n*. In periodical stimului MEG data are classified into 3 groups as

where **x***s*, **x***e* and **x***a* are spontaneous and evoked magnetic fields and artifact noises. The SEF

where *Le* is the lead field from the evoked current dipoles, **Qe**(*n*), **u**(*n*) := *r*350**w** is defined by

The decorrelation method was developed by Molgedey & Schuster (1994), Ziehe et al. (2000), Murata et al. (2001) and briefly summarized in Kishida et al. (2003). However, BSS performance is strongly dependent on the choice of time lag. The temporal decorrelation method of BSS has an open problem in choice of the time lag (Hironaga & Ioannides, 2007; Kishida, 2008; Tang et al., 2005; Ziehe et al., 2000). Here it should be noted that the

Since PSD of MEG has a line spectrum of 2Hz and those of higher harmonic modes in 2Hz periodical median nerve stimuli, improvements with the fractional type time lag *τm* defined

/*m* = 625

could be made on BSS (Kishida, 2008; 2009a). Hence, we can determine a mixing matrix *A*(*k*)

**x**(*n*) = *A*(*k*)**s**(*k*)

**x***e*(*n*) = *Le***Q***e*(*n*) = **u**(*n*) + *Le***Q***h*

evoked magnetic field consists of a deterministic part and fluctuations around it:

concatenating 350 copies of the average waveforms **w**, and **Q***h*

decorrelation method is effective for MEG data in a stationary process.

1250 2

by the blind source separation with the T/k type (fractional) time lag:

*τm* =

evoked magnetic field (see Eq. (16)).

**3.3 Fractional type of BSS**

by

**x**(*n*) = **x***s*(*n*) + **x***e*(*n*) + **x***a*(*n*), (1)

*e* (*n*), (2)

*e* (*n*) are SEF fluctuations in the

*m* , *m* = 1, 2, . . . , *k* (3)

(*n*). (4)

SQUID was the axial gradiometer type. Details were reported in Kishida (2009a).

**3. Methods**

**3.2 Data analysis**

SQUIDs be

**3.1 Magnetoencephalography (subjects and stimuli)**

In the T/k type decorrelation method, the absolute sum of off-diagonal elements of normalized correlation matrices is minimized at times corresponding to 2Hz and its higher harmonic frequencies. The main processing in the decorrelation method is removing the off-diagonal elements of the correlation matrices at fractional type time lag *τm*,

$$\mathbb{C}\_{zz}(\tau\_m) = \frac{1}{N} \sum\_{n=0}^{N-1} \mathbf{z}(n)\mathbf{z}(n + \tau\_m)^T \quad m = 1, \dots, k,\tag{5}$$

where the superscript *T* denotes the transposition of matrix and *z*(*n*) with zero mean is orthonormalized as √ *V*−1**x**(*n*). Here V is the covariance matrix given by *E*{**x**(*n*)**x**(*n*)*T*}. The Jacobi-like algorithm has been used to solve the simultaneous diagonalization problem approximately on k normalized correlation matrices (Cardoso & Souloumiac, 1996). This process is to determine a square matrix *U* in the problem of minimization of the cost function *J*,

$$f(\mathcal{U}) = \sum\_{m=1}^{k} \sum\_{i \neq j} |(\mathcal{U}\mathcal{C}\_{zz}(\tau\_m)\mathcal{U}^T)\_{ij}|^2 \,\tag{6}$$

where (*UCzz*(*τm*)*UT*)*ij* denotes the *ij*-element of matrix *UCzz*(*τm*)*UT*. Finally, we have *A*−1 (*k*) := *U* √ *V*−1.

For the classification of MEG we have used the T/k type of BSS with e.g. *k* = 30:

$$\mathbf{x}(n) = A\_{(30)} \mathbf{s}^{(30)}(n). \tag{7}$$

Then, MEG is decomposed by the fractional BSS with *k* = 30:

$$\mathbf{x}(n) = A\_{(30)}^{s} \mathbf{s}\_{s}^{(30)}(n) + A\_{(30)}^{\varepsilon} \mathbf{s}\_{\varepsilon}^{(30)}(n) + A\_{(30)}^{a} \mathbf{s}\_{a}^{(30)}(n),\tag{8}$$

where *A*(30) = (*As* (30) *Ae* (30) *Aa* (30) ) and **s**(30)(*n*)=(**s** (30) *s* (*n*)*T* **s** (30) *e* (*n*)*T* **s** (30) *a* (*n*)*T*)*T* corresponding to Eq. (1).

#### **3.4 Honest BSS and transformation matrix**

From the Wold decomposition theorem (Mathematical Society of Japan, 1987), a stationary process is uniquely decomposed into **x**(*n*) = **u**(*n*) + *δ***x**(*n*), where *δ***x**(n) is non-deterministic and **u**(*n*) is deterministic. We can define correlation functions in a stationary process, and they are specified into two types. When concatenate averaged waveforms are very small, the decorrelation method and the identification method used in Kishida (2009b) work efficiently, since correlation functions can be described by time lag. Inversely the decorrelation method and the identification method do not work well, when concatenate averaged waveforms are not small. In the remaining part of this paper we will first eliminate concatenation of averaged waveforms from MEG data, and then apply the decorrelation method and the identification method to eliminated MEG data. Such approach was called "honest" in the (9) of discussion of Kishida (2009b).

Supposing there is a subject who has no or negligible small fluctuations of SEF, and that he has averaged waveforms of SEF. However, certain fluctuations of SEF happen to be in **s** (30) *e* (*n*) of the subject in BSS calculation of Eq. (7). Therefore BSS of Eq. (7) may be called as "dishonest",

**3.5 Estimation of SEF fluctuations of current dipoles**

Kishida, 2006). That is, the covariance matrix, **D** = *E*{**x***h*

when lead fields of evoked magnetic field are selected by

Section 4.4 in the case of *q* = 5, where p=4 from Eq. (22).

assumed to be expressed by a feedback model defined by

when estimation of current dipole is true.

**3.6 Identification of feedback model**

have from Eqs. (9) and (13)

or,

eigenvectors:

The problem of multi-dipole estimation can be solved efficiently by the conventional least square method (Hämäläinen et al., 1993) with the combinatory optimization problem for location of current dipoles. Estimated current dipoles were equivalently evaluated by the least mean square method with a cost function of distance between signal spaces (Yokota &

**D** =

the *p* dimensional signal subspace **V***s* and the other subspace is noise subspace **V***n*:

For *q* current dipole estimation the cost function is defined (Yokota & Kishida, 2006) by

64 ∑ *k*=1

where *λk* is an eigenvalue and *φk* is its eigenvector. Let us divide into two subspaces; one is

**V***s* = (*φ*1, *φ*2,..., *φp*), **V***n* = (*φp*+1, *φp*+2,..., *φ*64).

**Lr**ˆ := (**l***r*ˆ1 , **l***r*ˆ2 ,..., **l***r*ˆ*q* ).

The supporting reason for cost function is that the rank of [(**V***s*, **Lr**ˆ)*T*(**V***s*, **Lr**ˆ)] becomes lower,

Hence, we can determine q current dipoles which represent dipole patterns of SEF. That is, we

*e Ae* (25)*δ***s**

where the symbol † means the Moore Penrose type of generalized inverse matrix and *Le* is a lead field evaluated from q current dipoles. Equation (16) will be expressed as Eq. (23) in

The identification method mentioned in Kishida (1996; 1997a; 2009b) could be applied to evaluate dynamical properties between current dipoles in cortices. Let current diploes

*y*1(*n*) = *F*12(*z*−1)*y*2(*n*) + *F*1(*z*−1)*f*1(*n*)

(25)*δ***s**

(25) *e* (*n*),

*y*2(*n*) = *F*21(*z*−1)*y*1(*n*) + *F*2(*z*−1)*f*2(*n*), (17)

*e* (*n*) = *Ae*

*e* (*n*) = *L*†

*Le***Q***h*

**Q***h*

corresponding to cSII and cSI be *y*1(*n*) and *y*2(*n*) selected from **Q***h*

*λkφkφT*

*e* (*n*)**x***h*

*J*(**r**ˆ1,...,**r**ˆ*q*ˆ) = det[(**V***s*, **Lr**ˆ)*T*(**V***s*, **Lr**ˆ)], (15)

*e* (*n*)*T*}, can be decomposed into

*k* , (14)

(25) *e* (*n*), (16)

*e* (*n*). Their activities are

since these fictitious fluctuations cause errors in identification method. As further mentioned in Discussion, subject D corresponds to the one who has no or negligible small fluctuations of SEF. To avoid these fictitious fluctuations, we should use the honest BSS instead of the dishonest BSS of Eq. (7), since there are no deterministic parts in the honest BSS. On the other hand, it is difficult to select honest BSS components of evoked field with dipole patterns in topographic field maps, since there is no line spectrum in PSD in honest BSS components. To overcome these difficulties, a transformation matrix between components of BSS with line spectrum and those of honest BSS will be introduced.

First let us eliminate concatenation of averaged waveforms from MEG data: *δ***x**(*n*) := **x**(*n*) − **u**(*n*) where **u**(*n*) := *r*350**w**. After eliminating concatenation of averaged waveforms the honest BSS and the identification algorithm can be applied. Fluctuations of MEG are decomposed by the honest BSS with e.g. *k* = 25 as the same method as in Section 3.3:

$$\begin{split} \delta \mathbf{x}(n) &= A\_{(25)}^{h} \delta \mathbf{s}^{(25)}(n) \\ &= \mathbf{x}(n) - \mathbf{u}(n) \\ &= \mathbf{x}\_{\mathcal{S}}(n) + L\_{\mathcal{C}} \mathbf{Q}\_{e}^{h}(n) + \mathbf{x}\_{\mathcal{A}}(n) \\ &= A\_{(25)}^{s} \delta \mathbf{s}\_{s}^{(25)}(n) + A\_{(25)}^{e} \delta \mathbf{s}\_{e}^{(25)}(n) + A\_{(25)}^{a} \delta \mathbf{s}\_{a}^{(25)}(n). \end{split} \tag{9}$$

In the honest BSS, there are no deterministic parts, and it would be difficult to find SEF BSS components, since there is no line spectrum of the evoked magnetic field in PSDs of honest BSS components and there were no one-to-one correspondences between dipoles and BSS components in four subjects. We must find a way to select SEF components of honest BSS.

From the two relations between Eq. (9) and Eq. (7) we can obtain approximately the next equation,

$$A\_{\left(30\right)}\mathbf{s}^{\left(30\right)}\left(\mathfrak{n}\right) = A\_{\left(25\right)}^{h}\delta\mathbf{s}^{\left(25\right)}\left(\mathfrak{n}\right),\tag{10}$$

when the amplitude of **u**(*n*) is negligible small in comparison with that of **x**(*n*). A transformation matrix is defined by

$$T := (A^h\_{(25)})^{-1} A\_{(30)}.\tag{11}$$

The *j*th column of *T* shows a distribution of the *j*th component of **s**(30)(*n*), *s* (30) *j* (*n*), on *δ***s**(25)(*n*) is denoted by *T*∗,*j*:

$$
\delta \mathbf{s}^{(25)}(n) = T \mathbf{s}^{(30)}(n). \tag{12}
$$

That is, *T*∗,*j* denotes the transformation from one component of **s**(30)(*n*) to a distribution on *δ***s**(25)(*n*). When components of honest BSS satisfy the condition that absolute value of *T*∗,*j* is close to one, the honest BSS components can be found from the dishonest BSS components. Therefore, fluctuations of SEF are given from the fractional BSS with *k* = 25 in the honest way by

$$\mathbf{x}\_{\varepsilon}^{h}(n) = A\_{(25)}^{\varepsilon} \delta \mathbf{s}\_{\varepsilon}^{(25)}(n). \tag{13}$$

#### **3.5 Estimation of SEF fluctuations of current dipoles**

The problem of multi-dipole estimation can be solved efficiently by the conventional least square method (Hämäläinen et al., 1993) with the combinatory optimization problem for location of current dipoles. Estimated current dipoles were equivalently evaluated by the least mean square method with a cost function of distance between signal spaces (Yokota & Kishida, 2006). That is, the covariance matrix, **D** = *E*{**x***h e* (*n*)**x***h e* (*n*)*T*}, can be decomposed into eigenvectors:

$$\mathbf{D} = \sum\_{k=1}^{64} \lambda\_k \phi\_k \phi\_{k\text{\textquotedbl{}}}^T \tag{14}$$

where *λk* is an eigenvalue and *φk* is its eigenvector. Let us divide into two subspaces; one is the *p* dimensional signal subspace **V***s* and the other subspace is noise subspace **V***n*:

$$\mathbf{V}\_s = (\phi\_1, \phi\_2, \dots, \phi\_p), \mathbf{V}\_n = (\phi\_{p+1}, \phi\_{p+2}, \dots, \phi\_{64}).$$

For *q* current dipole estimation the cost function is defined (Yokota & Kishida, 2006) by

$$J(\mathbf{\hat{r}}\_1, \dots, \mathbf{\hat{r}}\_{\vec{\mathbb{Q}}}) = \det[(\mathbf{V}\_{\text{s}}, \mathbf{L}\_{\mathbf{\hat{r}}})^T (\mathbf{V}\_{\text{s}}, \mathbf{L}\_{\mathbf{\hat{r}}})] \,\,\,\tag{15}$$

when lead fields of evoked magnetic field are selected by

$$\mathbf{L}\_{\mathbf{f}} := (\mathbf{l}\_{\mathbb{P}\_1}, \mathbf{l}\_{\mathbb{P}\_2}, \dots, \mathbf{l}\_{\mathbb{P}\_q})\,.$$

The supporting reason for cost function is that the rank of [(**V***s*, **Lr**ˆ)*T*(**V***s*, **Lr**ˆ)] becomes lower, when estimation of current dipole is true.

Hence, we can determine q current dipoles which represent dipole patterns of SEF. That is, we have from Eqs. (9) and (13)

$$L\_{\varepsilon} \mathbf{Q}\_{\varepsilon}^{h}(n) = A\_{(25)}^{\varepsilon} \delta \mathbf{s}\_{\varepsilon}^{(25)}(n),$$

$$\mathbf{Q}\_{\varepsilon}^{h}(n) = L\_{\varepsilon}^{\dagger} A\_{(25)}^{\varepsilon} \delta \mathbf{s}\_{\varepsilon}^{(25)}(n),\tag{16}$$

or,

6 Will-be-set-by-IN-TECH

since these fictitious fluctuations cause errors in identification method. As further mentioned in Discussion, subject D corresponds to the one who has no or negligible small fluctuations of SEF. To avoid these fictitious fluctuations, we should use the honest BSS instead of the dishonest BSS of Eq. (7), since there are no deterministic parts in the honest BSS. On the other hand, it is difficult to select honest BSS components of evoked field with dipole patterns in topographic field maps, since there is no line spectrum in PSD in honest BSS components. To overcome these difficulties, a transformation matrix between components of BSS with line

First let us eliminate concatenation of averaged waveforms from MEG data: *δ***x**(*n*) := **x**(*n*) − **u**(*n*) where **u**(*n*) := *r*350**w**. After eliminating concatenation of averaged waveforms the honest BSS and the identification algorithm can be applied. Fluctuations of MEG are decomposed by

*e* (*n*) + **x***a*(*n*)

(25) *δ***s**

In the honest BSS, there are no deterministic parts, and it would be difficult to find SEF BSS components, since there is no line spectrum of the evoked magnetic field in PSDs of honest BSS components and there were no one-to-one correspondences between dipoles and BSS components in four subjects. We must find a way to select SEF components of honest BSS. From the two relations between Eq. (9) and Eq. (7) we can obtain approximately the next

(*n*) = *Ah*

when the amplitude of **u**(*n*) is negligible small in comparison with that of **x**(*n*). A

(*n*) = *T***s**(30)

That is, *T*∗,*j* denotes the transformation from one component of **s**(30)(*n*) to a distribution on *δ***s**(25)(*n*). When components of honest BSS satisfy the condition that absolute value of *T*∗,*j* is close to one, the honest BSS components can be found from the dishonest BSS components. Therefore, fluctuations of SEF are given from the fractional BSS with *k* = 25 in the honest way

(25)*δ***s**

*T* := (*Ah*

(25)*δ***s**(25)

(25) *e* (*n*) + *Aa*

(25) *δ***s** (25) *a* (*n*).

(*n*), (10)

(30)

(*n*). (12)

(25) *e* (*n*). (13)

*j* (*n*), on *δ***s**(25)(*n*)

(25))−1*A*(30). (11)

(9)

(25) *s* (*n*) + *Ae*

spectrum and those of honest BSS will be introduced.

*δ***x**(*n*) = *Ah*

transformation matrix is defined by

equation,

is denoted by *T*∗,*j*:

by

= *As* (25) *δ***s**

the honest BSS with e.g. *k* = 25 as the same method as in Section 3.3:

*δ***s**(25)(*n*)

*A*(30)**s**(30)

The *j*th column of *T* shows a distribution of the *j*th component of **s**(30)(*n*), *s*

**x***h*

*e* (*n*) = *Ae*

*δ***s**(25)

(25)

= **x**(*n*) − **u**(*n*) = **x***s*(*n*) + *Le***Q***h*

> where the symbol † means the Moore Penrose type of generalized inverse matrix and *Le* is a lead field evaluated from q current dipoles. Equation (16) will be expressed as Eq. (23) in Section 4.4 in the case of *q* = 5, where p=4 from Eq. (22).

#### **3.6 Identification of feedback model**

The identification method mentioned in Kishida (1996; 1997a; 2009b) could be applied to evaluate dynamical properties between current dipoles in cortices. Let current diploes corresponding to cSII and cSI be *y*1(*n*) and *y*2(*n*) selected from **Q***h e* (*n*). Their activities are assumed to be expressed by a feedback model defined by

$$\begin{aligned} y\_1(n) &= F\_{12}(z^{-1})y\_2(n) + F\_1(z^{-1})f\_1(n) \\ y\_2(n) &= F\_{21}(z^{-1})y\_1(n) + F\_2(z^{-1})f\_2(n) \end{aligned} \tag{17}$$

**4. Results**

were invalid.

4 x 10−13

−2

−4

−2

0

w(T)

2

4 x 10−13

0

w(T)

2

**4.1 MEG pre-processing and SEF**

0 100 200 300 400 500 −4

Time(ms)

0 100 200 300 400 500

Time(ms)

(c)

(a)

After zero preset signal processing mentioned in Kishida (2009a), the average of 350 MEG responses was used. The averaged waveforms **w**(*n*) are obtained from SQUID time series data and shown in Fig. 1 for four subjects, where **w**(*n*) were defined by Eq. (1) of Kishida (2009a). Since the head of the other subject E moved during MEG measurement, the data

−2

−2

0

w(T)

Fig. 1. Averaged waveforms of SEF. (a) subject A, (b) subject B, (c) subject C and (d) subject D.

In Fig. 1(a), electric stimulus time is at time 101.6 ms (*n* = 127) and the first peak of averaged waveforms at time 120.8 ms (*n* = 151) is known as N20 with latency 19.2 ms (Hämäläinen et al., 1993; Kakigi et al., 2000; Wikström et al., 1996). N20 has a dipole pattern of cSI in the field map of Fig. 2(a). In 2Hz median nerve stimulus, averaged waveforms usually have a dipole pattern of cSI, since activities of cSI are found in 2Hz median nerve stimulus (Wikström et al., 1996). However, a dipole pattern of SIIs was not found in averaged waveforms at a glance, even if activities of SIIs exist in 2Hz median nerve stimulus (Wikström et al., 1996). This situation is similar to the case where stars are seen in the nighttime, but not in the daytime.

2

4 x 10−13

0

w(T)

2

4 x 10−13

0 100 200 300 400 500 −4

Time(ms)

(b)

0 100 200 300 400 500 −4

Time(ms)

(d)

where *z*−1 is the time shift operator: *z*−1*yk*(*n*) = *yk*(*n* − 1) (*k* = 1 or 2), *f*1(*n*) and *f*2(*n*) are Gaussian white random current dipoles in two regions of thalamus, *F*12(*z*−1) is an intracerebral transfer function from cSI to cSII, and *F*21(*z*−1) is a reverse intracerebral transfer function from cSII to cSI. *F*1(*z*−1) is a transfer function from one region of thalamus to cSII and *F*2(*z*−1) is that from another region of thalamus to cSI.

As mentioned in Eq. (A3) of Kishida (2009b), an innovation model of minimum phase (Kishida, 1991) can be identified from fluctuations of cSI and cSII:

$$\begin{aligned} \mathbf{x}(n|n) &= A\_H \mathbf{x}(n-1|n-1) + B\gamma(n) \\ \mathbf{y}(n) &= \mathbf{C} \mathbf{x}(n|n). \end{aligned} \tag{18}$$

Three system matrices of Eq. (18) are determined by the identification method from correlation function matrices, since correlation function matrices of Hankel matrix in Step 3 of Appendix A of Kishida (2009b) are obtained from output data of fluctuations of cSI and cSII. From Eq. (18) a closed loop transfer function matrix is

$$\mathbb{C}(z^{-1}) = \mathbb{C}(I - A\_H z^{-1})^{-1} \mathbb{B}. \tag{19}$$

Attention should be paid to the existence of an undetermined matrix U in the representation of the innovation model: *A*′ = U−1*A*U, *B*′ = U−1*B* and *C*′ = *C*U, although the closed loop transfer function matrix is invariant for U, e.g., eigenvalues of *A*′ is the same as those of *A*.

As in Eq. (B3) of Kishida (2009b), transfer functions between cSII and cSI are determined from Eq. (19) as

$$
\begin{aligned}
\hat{\mathbf{f}}\_{12}(z^{-1}) &= \mathbf{G}\_{12}(z^{-1}) / \mathbf{G}\_{22}(z^{-1}) \\
\hat{\mathbf{f}}\_{21}(z^{-1}) &= \mathbf{G}\_{21}(z^{-1}) / \mathbf{G}\_{11}(z^{-1}).
\end{aligned}
\tag{20}
$$

Transfer functions of *F*12(*z*−1) and *F*21(*z*−1) of intracerebral communication between cSI and cSII can be obtained via Eqs. (19) and (20) from the feedback structure under either one of the two sufficient conditions of original feedback model (Kishida, 1994). On the other hand, transfer functions of *F*1(*z*−1) and *F*2(*z*−1) are not identified from output data of fluctuations of cSI and cSII, since input data about fluctuations of thalamus are not available.

Our algorithm to obtain transfer functions between cSI and cSII is summarized in the following flowchart:

#### **4. Results**

8 Will-be-set-by-IN-TECH

where *z*−1 is the time shift operator: *z*−1*yk*(*n*) = *yk*(*n* − 1) (*k* = 1 or 2), *f*1(*n*) and *f*2(*n*) are Gaussian white random current dipoles in two regions of thalamus, *F*12(*z*−1) is an intracerebral transfer function from cSI to cSII, and *F*21(*z*−1) is a reverse intracerebral transfer function from cSII to cSI. *F*1(*z*−1) is a transfer function from one region of thalamus to cSII and

As mentioned in Eq. (A3) of Kishida (2009b), an innovation model of minimum phase

**x**(*n*|*n*) = *AH***x**(*n* − 1|*n* − 1) + *Bγ*(*n*)

Three system matrices of Eq. (18) are determined by the identification method from correlation function matrices, since correlation function matrices of Hankel matrix in Step 3 of Appendix A of Kishida (2009b) are obtained from output data of fluctuations of cSI and

Attention should be paid to the existence of an undetermined matrix U in the representation of the innovation model: *A*′ = U−1*A*U, *B*′ = U−1*B* and *C*′ = *C*U, although the closed loop transfer function matrix is invariant for U, e.g., eigenvalues of *A*′ is the same as those of *A*. As in Eq. (B3) of Kishida (2009b), transfer functions between cSII and cSI are determined from

12(*z*−1) = *G*12(*z*−1)/*G*22(*z*−1)

Transfer functions of *F*12(*z*−1) and *F*21(*z*−1) of intracerebral communication between cSI and cSII can be obtained via Eqs. (19) and (20) from the feedback structure under either one of the two sufficient conditions of original feedback model (Kishida, 1994). On the other hand, transfer functions of *F*1(*z*−1) and *F*2(*z*−1) are not identified from output data of fluctuations

Our algorithm to obtain transfer functions between cSI and cSII is summarized in the

honest BSS, Eq.(9) −−−−−−−−−−→ (dishonest BSS)

inverse problem, Eq. (16) −−−−−−−−−−−−−−→

identification method based −−−−−−−−−−−−−−−→ on feedback model, Eq. (17)

of cSI and cSII, since input data about fluctuations of thalamus are not available.

**y**(*n*) = *C***x**(*n*|*n*). (18)

*G*(*z*−1) = *C*(*I* − *AHz*−1)−1*B*. (19)

21(*z*−1) = *G*21(*z*−1)/*G*11(*z*−1). (20)

Extraction of SEF fluctuations. (see Eq. (13) or (22))

Fluctuations of current dipoles at cSI and cSII. (see Eqs. (24) and (25))

Transfer functions between cSI and cSII. (see Figs. 11 and 12)

*F*2(*z*−1) is that from another region of thalamus to cSI.

(Kishida, 1991) can be identified from fluctuations of cSI and cSII:

cSII. From Eq. (18) a closed loop transfer function matrix is

*F*ˆ

*F*ˆ

Eq. (19) as

following flowchart:

MEG data in 2Hz median nerve stimuli

> SEF fluctuations (see Eq. (13))

Correlation functions of SEF fluctuations at cSI and cSII

#### **4.1 MEG pre-processing and SEF**

After zero preset signal processing mentioned in Kishida (2009a), the average of 350 MEG responses was used. The averaged waveforms **w**(*n*) are obtained from SQUID time series data and shown in Fig. 1 for four subjects, where **w**(*n*) were defined by Eq. (1) of Kishida (2009a). Since the head of the other subject E moved during MEG measurement, the data were invalid.

Fig. 1. Averaged waveforms of SEF. (a) subject A, (b) subject B, (c) subject C and (d) subject D.

In Fig. 1(a), electric stimulus time is at time 101.6 ms (*n* = 127) and the first peak of averaged waveforms at time 120.8 ms (*n* = 151) is known as N20 with latency 19.2 ms (Hämäläinen et al., 1993; Kakigi et al., 2000; Wikström et al., 1996). N20 has a dipole pattern of cSI in the field map of Fig. 2(a). In 2Hz median nerve stimulus, averaged waveforms usually have a dipole pattern of cSI, since activities of cSI are found in 2Hz median nerve stimulus (Wikström et al., 1996). However, a dipole pattern of SIIs was not found in averaged waveforms at a glance, even if activities of SIIs exist in 2Hz median nerve stimulus (Wikström et al., 1996). This situation is similar to the case where stars are seen in the nighttime, but not in the daytime.

We can select BSS components of the evoked magnetic field from their averaged waveforms of **s**(30)(*n*) in Fig. 4. Averaged waveforms of BSS components were defined as *ws*(*n*) by Eq. (15) of Kishida (2009a). Cyan dotted and blue broken bold lines which have N20 peaks in Fig. 4 were the 4th and the 22nd BSS components of cSI, and yellow, red and black bold lines at time 214 ms (latency 112.4 ms) were the 13th, 17th and 35th BSS components of SIIs. The

(30) 17 (*n*)*s*

Their dipole patterns of BSS components corresponding to the somatosensory evoked field are shown in the topographic field maps of Fig. 5. Here, the topographic field map of the *j*th component of BSS is expressed by *Aj* as the *j*th column vector of *A*(30). In subject A there were no one-to-one correspondences between two dipoles and five BSS components. There were also no one-to-one correspondences in the other subjects. It was a rare case to have the

0 100 200 300 400 500

Time(ms)

Fig. 4. BSS waveforms of subject A. Cyan bold broken and blue dotted lines are the 4th and the 22nd BSS components of cSI, and three bold lines are the 13th, 17th and 35th BSS

PSD of MEG with elimination of repeated waveforms has no line spectrum of 2Hz and no higher harmonic modes in Fig. 6. Here the SQUID channel of Fig. 6 is the same as that of Fig. 3. From Fig. 6 electrical power noise is clearly found in line structure at 60Hz and has higher

(30) 22 (*n*)*s*

(30)

35 (*n*))*T*. (21)

somatosensory evoked field was given by five components of BSS:

one-to-one correspondences between dipoles and BSS components.

(30) 13 (*n*)*s*

**s** (30) *e* (*n*) = (*s* (30) 4 (*n*)*s*

components of SIIs.

**4.2 Honest BSS**

harmonics.

ws

0

0.5

1

1.5

2

Fig. 2. Dipole pattern in the topographic field map of waveforms of subject A. (a) time 120.8ms (latency 19.2 ms), (b) time 214.4ms (latency 112.8 ms).

SIIs correspond to stars and cSI is the sun. The nighttime is a time when the sign of amplitude of cSI changes. The time is 214.4 ms corresponding to latency 112.8 ms in Fig. 1(a), and the dipole pattern of SIIs is found from **w**(*n*) at the time as in the field map of Fig. 2(b). From Fig. 2(b) cSII is stronger than the ipsilateral secondary somatosensory cortex (iSII). This result is consistent to that of Wikström et al. (1996).

The black bold line in Fig. 1(a) is an averaged waveform at the 14th SQUID channel, location of which will be shown by a black large circle in Fig. 10. PSD of SQUID data at the 14th channel has a line spectrum of 2Hz and those of higher harmonic modes for 2Hz periodical median nerve stimuli as in Fig. 3. The line spectrum of PSD is due to periodical structure of the somatosensory evoked magnetic field. The 2Hz repetitive line spectra of Fig. 3 show an effect of the deterministic and nonlinear part.

Fig. 3. PSD of the 14th SQUID channel of subject A.

10 Will-be-set-by-IN-TECH

Fig. 2. Dipole pattern in the topographic field map of waveforms of subject A. (a) time

SIIs correspond to stars and cSI is the sun. The nighttime is a time when the sign of amplitude of cSI changes. The time is 214.4 ms corresponding to latency 112.8 ms in Fig. 1(a), and the dipole pattern of SIIs is found from **w**(*n*) at the time as in the field map of Fig. 2(b). From Fig. 2(b) cSII is stronger than the ipsilateral secondary somatosensory cortex (iSII). This result is

The black bold line in Fig. 1(a) is an averaged waveform at the 14th SQUID channel, location of which will be shown by a black large circle in Fig. 10. PSD of SQUID data at the 14th channel has a line spectrum of 2Hz and those of higher harmonic modes for 2Hz periodical median nerve stimuli as in Fig. 3. The line spectrum of PSD is due to periodical structure of the somatosensory evoked magnetic field. The 2Hz repetitive line spectra of Fig. 3 show an

1 10 100 −320

f(Hz)

P

B

(b)

P

L R

FLR

L R

A

A

P

(a)

consistent to that of Wikström et al. (1996).

effect of the deterministic and nonlinear part.

−250

−310

Fig. 3. PSD of the 14th SQUID channel of subject A.

−300

−290

−280

P(dB/Hz)

−270

−260

120.8ms (latency 19.2 ms), (b) time 214.4ms (latency 112.8 ms).

L R PA L R

A

We can select BSS components of the evoked magnetic field from their averaged waveforms of **s**(30)(*n*) in Fig. 4. Averaged waveforms of BSS components were defined as *ws*(*n*) by Eq. (15) of Kishida (2009a). Cyan dotted and blue broken bold lines which have N20 peaks in Fig. 4 were the 4th and the 22nd BSS components of cSI, and yellow, red and black bold lines at time 214 ms (latency 112.4 ms) were the 13th, 17th and 35th BSS components of SIIs. The somatosensory evoked field was given by five components of BSS:

$$\begin{array}{l} \mathbf{s}\_{\varepsilon}^{(30)}(n) \\ = \left( s\_{4}^{(30)}(n) s\_{13}^{(30)}(n) \, s\_{17}^{(30)}(n) \, s\_{22}^{(30)}(n) \, s\_{35}^{(30)}(n) \right)^{T} . \end{array} \tag{21}$$

Their dipole patterns of BSS components corresponding to the somatosensory evoked field are shown in the topographic field maps of Fig. 5. Here, the topographic field map of the *j*th component of BSS is expressed by *Aj* as the *j*th column vector of *A*(30). In subject A there were no one-to-one correspondences between two dipoles and five BSS components. There were also no one-to-one correspondences in the other subjects. It was a rare case to have the one-to-one correspondences between dipoles and BSS components.

Fig. 4. BSS waveforms of subject A. Cyan bold broken and blue dotted lines are the 4th and the 22nd BSS components of cSI, and three bold lines are the 13th, 17th and 35th BSS components of SIIs.

#### **4.2 Honest BSS**

PSD of MEG with elimination of repeated waveforms has no line spectrum of 2Hz and no higher harmonic modes in Fig. 6. Here the SQUID channel of Fig. 6 is the same as that of Fig. 3. From Fig. 6 electrical power noise is clearly found in line structure at 60Hz and has higher harmonics.

1 10 100 −320

Fig. 6. PSD of the 14th SQUID channel of subject A, after elimination of concatenation of

the condition |*T*∗,*j*| > 0.6, the dishonest BSS components are transformed into those of honest BSS with *k* = 25: they are transformed into the 4th, 27th, 36th and 43rd components of honest

Their topographic field maps of honest BSS components are shown in Fig. 8. Here, the

Figs. 5 and 8. If the condition |*T*∗,*j*| > 0.6 is satisfied at the least, the correspondence between the dishonest and honest BSS components is plausible. However, there was the exceptional example (see Section 5.1) that the correspondence did not hold even when the condition was

The honest BSS is a blind source separation which relates to the block diagonal of covariance matrix. The block diagonal property is found in cross-correlation functions among honest BSS components of the evoked magnetic field. Cross-correlation functions between the 4th component of honest BSS and the other component of honest BSS are shown in Fig. 9(a), and especially those between the 4th component and the 27th, 36th and 43rd components of honest BSS are shown by yellow broken, blue 1-dot broken and blue solid bold lines in Fig. 9(b). In Fig. 9 large deviations of cross-correlation away from the time origin are found occasionally. Therefore it is concluded that there remain small correlations among honest BSS components

Hence, fluctuations of SEF are given from the fractional BSS with *k* = 25 in the honest way by

(25) 27 (*n*) *δs*

(25)*δ***s**

(25) 36 (*n*) *δs*

(25) *e* (*n*), (22)

(25) 43 (*n*))*T*.

**x***h*

(25) 4 (*n*) *δs*

*e* (*n*) = *Ae*

(25) in the honest way. Similarities are found in patterns of dipoles between

*j* as the *j*th

topographic field map of the *j*th component of honest BSS is also denoted by *Ah*

**4.3 Fluctuations of somatosensory evoked or induced field**

f (Hz)

−310

averaged waveforms from MEG.

BSS with *k* = 25.

column vector of *Ah*

corresponding to cSI and SIIs.

*δ***s**

(25) *e* (*n*)=(*δs*

satisfied.

where

−300

−290

−280

P (dB/Hz)

−270

−260

−250

When components of honest BSS satisfy the condition that absolute value of *T*∗,*j* is close to one, the dishonest BSS components are transformed into honest BSS components with similar topographic field maps. In Fig. 7 *T*∗,4, *T*∗,13, *T*∗,17, *T*∗,22 and *T*∗,35 of the 4th, 13rd, 17th, 22rd and 35th components of BSS with *k* = 30 are denoted by bold solid, bold 1-dot broken, bold broken, dotted and sold lines respectively. When the components of BSS with *k* = 30 satisfy 12 Will-be-set-by-IN-TECH

P

B

(b) *A*13.

A

A

P

P

(d) *A*22.

A

A

P

P

(c) *A*17.

A

P

(e) *A*35.

Fig. 5. Topographic field maps of SEF BSS components in subject A.

When components of honest BSS satisfy the condition that absolute value of *T*∗,*j* is close to one, the dishonest BSS components are transformed into honest BSS components with similar topographic field maps. In Fig. 7 *T*∗,4, *T*∗,13, *T*∗,17, *T*∗,22 and *T*∗,35 of the 4th, 13rd, 17th, 22rd and 35th components of BSS with *k* = 30 are denoted by bold solid, bold 1-dot broken, bold broken, dotted and sold lines respectively. When the components of BSS with *k* = 30 satisfy

L R

L R

L R

P

L R

FLR

L R

A

A

Fig. 6. PSD of the 14th SQUID channel of subject A, after elimination of concatenation of averaged waveforms from MEG.

the condition |*T*∗,*j*| > 0.6, the dishonest BSS components are transformed into those of honest BSS with *k* = 25: they are transformed into the 4th, 27th, 36th and 43rd components of honest BSS with *k* = 25.

Their topographic field maps of honest BSS components are shown in Fig. 8. Here, the topographic field map of the *j*th component of honest BSS is also denoted by *Ah j* as the *j*th column vector of *Ah* (25) in the honest way. Similarities are found in patterns of dipoles between Figs. 5 and 8. If the condition |*T*∗,*j*| > 0.6 is satisfied at the least, the correspondence between the dishonest and honest BSS components is plausible. However, there was the exceptional example (see Section 5.1) that the correspondence did not hold even when the condition was satisfied.

#### **4.3 Fluctuations of somatosensory evoked or induced field**

The honest BSS is a blind source separation which relates to the block diagonal of covariance matrix. The block diagonal property is found in cross-correlation functions among honest BSS components of the evoked magnetic field. Cross-correlation functions between the 4th component of honest BSS and the other component of honest BSS are shown in Fig. 9(a), and especially those between the 4th component and the 27th, 36th and 43rd components of honest BSS are shown by yellow broken, blue 1-dot broken and blue solid bold lines in Fig. 9(b). In Fig. 9 large deviations of cross-correlation away from the time origin are found occasionally. Therefore it is concluded that there remain small correlations among honest BSS components corresponding to cSI and SIIs.

Hence, fluctuations of SEF are given from the fractional BSS with *k* = 25 in the honest way by

$$\mathbf{x}\_{\varepsilon}^{h}(n) = A\_{(25)}^{\varepsilon} \delta \mathbf{s}\_{\varepsilon}^{(25)}(n),\tag{22}$$

where

$$\delta \mathbf{s}\_{\varepsilon}^{(25)}(n) = (\delta \mathbf{s}\_{4}^{(25)}(n) \, \delta \mathbf{s}\_{27}^{(25)}(n) \, \delta \mathbf{s}\_{36}^{(25)}(n) \, \delta \mathbf{s}\_{43}^{(25)}(n))^{T}.$$

P

A

F

P

Fig. 8. Topographic field maps of honest SEF BSS components in subject A.

*e* (*n*) = *L*†

*e Ae* (25)*δ***s**

*Qh* 62(*n*) *Qh* 363(*n*) *Qh* 1210(*n*) *Qh* 1835(*n*) *Qh* 7706(*n*) ,

**L***e* = (*l*62, *l*363, *l*1210, *l*1835, *l*7706),

.

From location of current dipoles in Fig. 10 cSII corresponds to the current dipole at position 1835: fluctuation *y*1(*n*) of cSII with yellow symbol in Fig.10 is mixed by a matrix *L*†

**Q***h*

**Q***h e* (*n*) :=

Here, notations of the above equation were reported in Kishida (2009b).

*l*1(∗, *x*) cos *θ*∗ + *l*1(∗, *y*) sin *θ*∗ *l*2(∗, *x*) cos *θ*∗ + *l*2(∗, *y*) sin *θ*∗ . . . *l*64(∗, *x*) cos *θ*∗ + *l*64(∗, *y*) sin *θ*∗

B

(c) *Ah* 36.

L R PA L R P

A

F

P

(25) *e* (*n*), (23)

*e Ae* (25)

B

(d) *Ah* 43.

L R PA L R

B

(b) *Ah* 27.

L R PA L R

L

L

A

F

B

(a) *Ah* 4.

L R PA L R

L

L

where

with *l*∗ :=

A

F

Fig. 7. Transition from dishonest SEF BSS components with *k* = 30 to honest SEF BSS components with *k* = 25.

Induced fields generated by the current dipoles at the same locations in Fig. 8 may be included in fluctuations.

#### **4.4 Source analysis**

If a brain can be approximated by a spherically symmetric model, and if the primary currents generated in the brain can be described by an equivalent current dipole, a lead field from a current dipole in the brain to SQUID channels was given by Sarvas formula (Sarvas, 1987).  They are expressed by tiny symbol of square in Fig. 10. The number of them was 9081. Large circles in Fig. 10 are 64 SQUIDs. Here the black large circle in Fig. 10 is the 14th SQUID channel mentioned in Figs. 3 and 6.

Five current dipoles were estimated by applying the method of dipole current estimation (see Section 3.5). 66.8 % of the covariance matrix of **D** = *E*{**x***h e* (*n*)**x***h e* (*n*)*T*} could be expressed by five current dipoles. Five current dipoles at position 7706 with angle 9.06 degrees, position 1835 with angle 100.27 degrees, position 1210 with angle 108.72 degrees, position 363 with angle 130.26 degrees and position 62 with 36.98 degrees are shown by symbols with magenta, yellow, blue, green and red respectively in Fig. 10. That is, we have

(b) *Ah* 27.

Fig. 8. Topographic field maps of honest SEF BSS components in subject A.

$$\mathbf{Q}\_{\varepsilon}^{h}(n) = L\_{\varepsilon}^{\dagger} A\_{(25)}^{\varepsilon} \delta \mathbf{s}\_{\varepsilon}^{(25)}(n),\tag{23}$$

where

14 Will-be-set-by-IN-TECH

0 10 20 30 40 50 60 70 -1

Induced fields generated by the current dipoles at the same locations in Fig. 8 may be included

If a brain can be approximated by a spherically symmetric model, and if the primary currents generated in the brain can be described by an equivalent current dipole, a lead field from a current dipole in the brain to SQUID channels was given by Sarvas formula (Sarvas, 1987).  They are expressed by tiny symbol of square in Fig. 10. The number of them was 9081. Large circles in Fig. 10 are 64 SQUIDs. Here the black large circle in Fig. 10 is the 14th SQUID channel mentioned in

Five current dipoles were estimated by applying the method of dipole current estimation (see

five current dipoles. Five current dipoles at position 7706 with angle 9.06 degrees, position 1835 with angle 100.27 degrees, position 1210 with angle 108.72 degrees, position 363 with angle 130.26 degrees and position 62 with 36.98 degrees are shown by symbols with

magenta, yellow, blue, green and red respectively in Fig. 10. That is, we have

*e* (*n*)**x***h*

*e* (*n*)*T*} could be expressed by

Section 3.5). 66.8 % of the covariance matrix of **D** = *E*{**x***h*

Fig. 7. Transition from dishonest SEF BSS components with *k* = 30 to honest SEF BSS

in fluctuations.

Figs. 3 and 6.

**4.4 Source analysis**

components with *k* = 25.

0

0.2

0.4

0.6

0.8

1

$$\mathbf{Q}\_{\varepsilon}^{h}(n) := \begin{pmatrix} \mathbf{Q}\_{62}^{h}(n) \\ \mathbf{Q}\_{363}^{h}(n) \\ \mathbf{Q}\_{1210}^{h}(n) \\ \mathbf{Q}\_{1835}^{h}(n) \\ \mathbf{Q}\_{7706}^{h}(n) \end{pmatrix},$$

$$\begin{aligned} \mathbf{L}\_{\varepsilon} &= (\underline{l}\_{62'} \underline{l}\_{363'} \underline{l}\_{1210'} \underline{l}\_{1835'} \underline{l}\_{7706}), \\ \text{with } \underline{l}\_{\*} := \begin{pmatrix} l\_1(\*, x) \cos \theta\_\* + l\_1(\*, y) \sin \theta\_\* \\ l\_2(\*, x) \cos \theta\_\* + l\_2(\*, y) \sin \theta\_\* \\ \vdots \\ l\_{64}(\*, x) \cos \theta\_\* + l\_{64}(\*, y) \sin \theta\_\* \end{pmatrix}. \end{aligned}$$

Here, notations of the above equation were reported in Kishida (2009b).

From location of current dipoles in Fig. 10 cSII corresponds to the current dipole at position 1835: fluctuation *y*1(*n*) of cSII with yellow symbol in Fig.10 is mixed by a matrix *L*† *e Ae* (25)

−0.1

−0.05

0

L AP

0.05

z(m)

(top) Axial view (bottom) Coronal view.

position 7706 is also expressed by

0.1

0.15

−0.1 −0.05 0 0.05 0.1

−0.1 −0.05 0 0.05 0.1

Fig. 10. Location of SQUIDs and current dipoles of subject A in the MEG head coordinates.

Fluctuation *y*2(*n*) of cSI with magenta symbol in Fig.10 corresponds to the current dipole at

*y*2(*n*) = *Qh*

x(m)

x(m)

P

A

R

R

7706(*n*). (25)

−0.05

y(m)

0

L

0.05

0.1

Fig. 9. Cross-correlation functions between the 4th honest BSS and the other honest BSS of subject A in the honest way. (a) All cross-correlation functions. (b) Three cross-correlation functions.

from  $\delta\_{\varepsilon}^{(25)}(n)$ , and is evaluated from Eq. (16) or (23) as
$$y\_1(n) = Q\_{1835}^{\natural}(n). \tag{24}$$
.

16 Will-be-set-by-IN-TECH

−120 −80 −40 −0 40 80 120

(a)

−120 −80 −40 0 40 80 120

(b) Fig. 9. Cross-correlation functions between the 4th honest BSS and the other honest BSS of subject A in the honest way. (a) All cross-correlation functions. (b) Three cross-correlation

*y*1(*n*) = *Qh*

(25) *e* (*n*), and is evaluated from Eq. (16) or (23) as

n(ms)

1835(*n*). (24)

n(ms)

−0.1

−0.1

functions. from *δ***s**

−0.05

0

0.05

0.1

−0.05

0

0.05

0.1

Fig. 10. Location of SQUIDs and current dipoles of subject A in the MEG head coordinates. (top) Axial view (bottom) Coronal view.

Fluctuation *y*2(*n*) of cSI with magenta symbol in Fig.10 corresponds to the current dipole at position 7706 is also expressed by

$$y\_2(n) = Q\_{7706}^h(n). \tag{25}$$

1 10 100 300 −60

F12

1 10 100 300 −400

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 −0.1

Since no suitable component of cSI in the honest BSS of subject D could be found, there are no

From averaged waveforms of subject D shown in Fig. 1(d), dipole patterns at time 120 and 180 ms are shown in the field maps of Figs. 14(a) and 14(b). The dipole pattern at time 120 or 180 ms was that of cSI. In 2Hz median nerve stimulus waveforms, it corresponds to what usually have a dipole pattern of cSI (Wikström et al., 1996). Then the fourth BSS component with the dipole pattern similar to that of cSI could be produced, as shown in Fig. 15(a). Transformation matrix from the fourth BSS component to components of honest BSS is given by *T*∗,4 in Fig. 15(c). *T*∗,4 denotes the 32nd component of honest BSS with *k* = 25 corresponding to the fourth component of BSS with *k* = 30. However, there is no similarity between dipole pattern

fluctuations of cSI in subject D does not succeed. This means that fluctuations of SEF in subject D were small and they have not been separated from MEG data by the honest BSS. Hence we

32 in Fig. 15(b). Separating the

Fig. 11. Bode diagram of *F*12 from cSI to cSII and its impulse response of subject A.

impulse responses in subject D. This will be explained in Discussion 5.1.

of *A*4 in the topographic field map in Fig. 15(a) and that of *Ah*

Time (sec)

Frequency(Hz)

Impulse Response

Bode diagram

−40

−200

−0.05 0 0.05 0.1 0.15

Amplitude

**5. Discussion**

**5.1 Separation of fluctuations of SEF**

have no impulse responses of subject D.

Phase(deg)

0

200

−20

Magnitude(dB)

0

Finally, fluctuations of cSI and cSII of feedback model were given by

$$\mathbf{y}(n) = \begin{pmatrix} y\_1(n) \\ y\_2(n) \end{pmatrix}.$$

#### **4.5 Intracerebral communication between cSI and cSII**

Coefficient matrices of the data-oriented innovation model of Eq. (18) must satisfy the minimum phase properties which are the pole stability and the zero invertibility (Kishida, 1991; 1996). Otherwise, the data-oriented innovation model cannot produce correct results theoretically. Pole stability and zero invertibility mean that roots of |*I* − *AHz*−1| = 0 and those of |*I* − *AH*(*I* − *BC*)*z*−1| = 0 are outside the unit circle in the complex plane of *z*−1. Or Eigenvalues of matrices *AH* and *AH*(*I* − *BC*) are inside the unit circle in the complex plane. The feedback model of Eq. (17) has excellent advantages for identification of MEG (Kishida, 2005; 2009b).

The number of singular values of Hankel matrix was 10, when the data-oriented innovation model was determined under the condition of minimum phase properties. Transfer functions of *F*12(*z*−1) and *F*21(*z*−1) of intracerebral communication between cSI and cSII could be obtained via Eqs. (19) and (20) from the feedback structure under the condition that a sufficient condition was satisfied in original feedback model (Kishida, 1994; 1996). The bode diagrams of *F*ˆ 12(*z*−1), *F*ˆ 21(*z*−1) and their impulse responses are shown in Figs. 11 and 12. The top and middle panel of figure is the magnitude and the phase of Bode diagram of each transfer function, and the bottom panel of figure is the impulse response. There are two modes from cSI to cSII found in Fig. 11. One is a rapid mode which is within about 5ms and the other is a slow mode with about 30ms order. On the other hand a rapid mode from cSII to cSI within about 5 ms is found in Fig. 12.

The reproducibility of the modes in impulse responses of intracerebral communication between cSI and cSII are now examined. In subject B, cSI and SIIs were represented by the 2nd, 4th, 14th, 19th and 20th components of dishonest BSS from the fractional BSS with *k* = 30. They were transformed into the 2nd, 22nd, 40th, 42nd and 57th components of honest BSS with *k* = 25. Six current dipoles were estimated from the same method of dipole current estimation. 76.2 % of the covariance matrix of *D* of subject B could be expressed by six current dipoles. The number of singular values of Hankel matrix was 9, when the data-oriented innovation model was determined. In Fig. 13, obtained impulse responses are shown by bold dotted lines, while bold solid lines indicate impulse responses of subject A in Fig. 13.

In subject C, cSI and SIIs were represented by the 4th, 12th, 14th, 33rd and 35th components of dishonest BSS from the fractional BSS with *k* = 30. They were transformed into the 10th, 13th, 14th, 35th and 36th components of honest BSS with *k* = 25. Seven current dipoles were also estimated. 75.7 % of the covariance matrix of *D* of subject C could be expressed by seven current dipoles. The number of singular values of Hankel matrix was 10, when the data-oriented innovation model was determined. The impulse responses are shown by bold broken lines in Fig. 13.

In Fig. 13, impulse responses of three subjects are summarized for clarifying the intracerebral communication between cSI and cSII. Two modes from cSI to cSII were found in Fig. 13(a). One rapid mode from cSII to cSI within 5 ms was mainly found in Fig. 13(b).

Fig. 11. Bode diagram of *F*12 from cSI to cSII and its impulse response of subject A.

Since no suitable component of cSI in the honest BSS of subject D could be found, there are no impulse responses in subject D. This will be explained in Discussion 5.1.

#### **5. Discussion**

18 Will-be-set-by-IN-TECH

Coefficient matrices of the data-oriented innovation model of Eq. (18) must satisfy the minimum phase properties which are the pole stability and the zero invertibility (Kishida, 1991; 1996). Otherwise, the data-oriented innovation model cannot produce correct results theoretically. Pole stability and zero invertibility mean that roots of |*I* − *AHz*−1| = 0 and those of |*I* − *AH*(*I* − *BC*)*z*−1| = 0 are outside the unit circle in the complex plane of *z*−1. Or Eigenvalues of matrices *AH* and *AH*(*I* − *BC*) are inside the unit circle in the complex plane. The feedback model of Eq. (17) has excellent advantages for identification of MEG (Kishida,

The number of singular values of Hankel matrix was 10, when the data-oriented innovation model was determined under the condition of minimum phase properties. Transfer functions of *F*12(*z*−1) and *F*21(*z*−1) of intracerebral communication between cSI and cSII could be obtained via Eqs. (19) and (20) from the feedback structure under the condition that a sufficient condition was satisfied in original feedback model (Kishida, 1994; 1996). The bode diagrams

and middle panel of figure is the magnitude and the phase of Bode diagram of each transfer function, and the bottom panel of figure is the impulse response. There are two modes from cSI to cSII found in Fig. 11. One is a rapid mode which is within about 5ms and the other is a slow mode with about 30ms order. On the other hand a rapid mode from cSII to cSI within

The reproducibility of the modes in impulse responses of intracerebral communication between cSI and cSII are now examined. In subject B, cSI and SIIs were represented by the 2nd, 4th, 14th, 19th and 20th components of dishonest BSS from the fractional BSS with *k* = 30. They were transformed into the 2nd, 22nd, 40th, 42nd and 57th components of honest BSS with *k* = 25. Six current dipoles were estimated from the same method of dipole current estimation. 76.2 % of the covariance matrix of *D* of subject B could be expressed by six current dipoles. The number of singular values of Hankel matrix was 9, when the data-oriented innovation model was determined. In Fig. 13, obtained impulse responses are shown by bold dotted lines, while bold solid lines indicate impulse responses of subject A in Fig. 13. In subject C, cSI and SIIs were represented by the 4th, 12th, 14th, 33rd and 35th components of dishonest BSS from the fractional BSS with *k* = 30. They were transformed into the 10th, 13th, 14th, 35th and 36th components of honest BSS with *k* = 25. Seven current dipoles were also estimated. 75.7 % of the covariance matrix of *D* of subject C could be expressed by seven current dipoles. The number of singular values of Hankel matrix was 10, when the data-oriented innovation model was determined. The impulse responses are shown by bold

In Fig. 13, impulse responses of three subjects are summarized for clarifying the intracerebral communication between cSI and cSII. Two modes from cSI to cSII were found in Fig. 13(a).

One rapid mode from cSII to cSI within 5 ms was mainly found in Fig. 13(b).

21(*z*−1) and their impulse responses are shown in Figs. 11 and 12. The top

 *y*1(*n*) *y*2(*n*)  .

**y**(*n*) =

Finally, fluctuations of cSI and cSII of feedback model were given by

**4.5 Intracerebral communication between cSI and cSII**

2005; 2009b).

12(*z*−1), *F*ˆ

about 5 ms is found in Fig. 12.

broken lines in Fig. 13.

of *F*ˆ

#### **5.1 Separation of fluctuations of SEF**

From averaged waveforms of subject D shown in Fig. 1(d), dipole patterns at time 120 and 180 ms are shown in the field maps of Figs. 14(a) and 14(b). The dipole pattern at time 120 or 180 ms was that of cSI. In 2Hz median nerve stimulus waveforms, it corresponds to what usually have a dipole pattern of cSI (Wikström et al., 1996). Then the fourth BSS component with the dipole pattern similar to that of cSI could be produced, as shown in Fig. 15(a). Transformation matrix from the fourth BSS component to components of honest BSS is given by *T*∗,4 in Fig. 15(c). *T*∗,4 denotes the 32nd component of honest BSS with *k* = 25 corresponding to the fourth component of BSS with *k* = 30. However, there is no similarity between dipole pattern of *A*4 in the topographic field map in Fig. 15(a) and that of *Ah* 32 in Fig. 15(b). Separating the fluctuations of cSI in subject D does not succeed. This means that fluctuations of SEF in subject D were small and they have not been separated from MEG data by the honest BSS. Hence we have no impulse responses of subject D.

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

Impulse Response

Time (sec)

Impulse Response

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

Time (sec)

(b)

Fig. 13. Impulse responses between cSI and cSII by using the honest BSS (Kishida, 2011). (a)

Impulse responses from cSI to cSII (*F*12). (b) Impulse responses from cSII to cSI (*F*21).

(a)

−0.15

−0.15

−0.1

−0.05

0

0.05

Amplitude

0.1

0.15

0.2

−0.1

−0.05

0

0.05

Amplitude

0.1

0.15

0.2

0.25

0.3

Fig. 12. Bode diagram of *F*21 from cSII to cSI and its impulse response of subject A.

#### **5.2 Comparison between dishonest BSS and honest BSS**

In our research, honest BSS was used after elimination of concatenate averaged waveforms, and impulse responses from fluctuations of SEF were identified.

There is another way to determine impulse responses: If the dishonest BSS is used first, SEF data can be selected and then it is possible to eliminate concatenate averaged waveforms from SEF data before identification of impulse responses. Two ways are summarized as

$$\begin{array}{ccccc}\mathbf{x}(n) & \implies & \mathbf{x}\_{\varepsilon}(n) \\ \downarrow\_{E} & & \downarrow\_{E} \\ \delta\mathbf{x}(n) & \implies A\_{h}^{\varepsilon}\mathbf{s}\_{\varepsilon}^{h}(n) \neq A\_{d}^{\varepsilon}\mathbf{s}\_{\varepsilon}^{d}(n) - \mathbf{u}(n). \end{array} \tag{26}$$

Here *Ae d* is defined by *Ae* (30) in Eq. (8) and **s***h e* corresponds to *δ***s** (25) *e* in Eq. (9). In Eq. (26), the symbol "=⇒" means application of BSS and selection of SEF components, and the symbol "↓*E*" is elimination of concatenate averaged waveforms.

In comparison with Fig. 13, Fig. 16 shows impulse responses in the dishonest way. When calculating SEF fluctuations in the dishonest way, the same positions of dipole currents have been used. Robustness of impulse responses of Fig. 16 is lower in comparison with Fig. 13.

20 Will-be-set-by-IN-TECH

Bode diagram

1 10 100 300 −35

F21

1 10 100 300 −60

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 −0.04

In our research, honest BSS was used after elimination of concatenate averaged waveforms,

There is another way to determine impulse responses: If the dishonest BSS is used first, SEF data can be selected and then it is possible to eliminate concatenate averaged waveforms from

↓*E* ↓*E*

symbol "=⇒" means application of BSS and selection of SEF components, and the symbol "↓*E*"

In comparison with Fig. 13, Fig. 16 shows impulse responses in the dishonest way. When calculating SEF fluctuations in the dishonest way, the same positions of dipole currents have been used. Robustness of impulse responses of Fig. 16 is lower in comparison with Fig. 13.

*e* (*n*) �= *Ae*

*d***s***d*

*e* corresponds to *δ***s**

*h***s***h*

**x***e*(*n*) (see Eq.(2))

*e* (*n*) − **u**(*n*).

(26)

(25) *e* in Eq. (9). In Eq. (26), the

Fig. 12. Bode diagram of *F*21 from cSII to cSI and its impulse response of subject A.

SEF data before identification of impulse responses. Two ways are summarized as

**x**(*n*) =⇒

(30) in Eq. (8) and **s***h*

**5.2 Comparison between dishonest BSS and honest BSS**

and impulse responses from fluctuations of SEF were identified.

*δ***x**(*n*) (see Eq. (9)) =⇒ *Ae*

is elimination of concatenate averaged waveforms.

Time (sec)

Frequency(Hz)

Impulse Response

−30

−40

−0.02 0 0.02 0.04 0.06

Amplitude

Here *Ae*

*d* is defined by *Ae*

−20

Phase(deg)

0

−25

Magnitude(dB)

−20

Fig. 13. Impulse responses between cSI and cSII by using the honest BSS (Kishida, 2011). (a) Impulse responses from cSI to cSII (*F*12). (b) Impulse responses from cSII to cSI (*F*21).

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 −0.2

(a)

Time (sec)

Impulse Response

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

Time (sec)

(b)

Although separate fluctuations of SEF in subject D could not be separated, the other three subjects could be separated, and the intracerebral communication between cSI and cSII was

Fig. 16. Impulse responses between cSI and cSII by using the dishonest BSS. (a) Impulse

responses from cSI to cSII (*F*12). (b) Impulse responses from cSII to cSI (*F*21).

Impulse Response

−0.15

−0.15

**5.3 SEF with steady state responses**

−0.1

−0.05

0

0.05

Amplitude

0.1

0.15

0.2

−0.1

−0.05

Amplitude

0

0.05

0.1

0.15

Fig. 14. Dipole patterns in the topographic field map of averaged waveforms of subject D in Fig. 1(d). (a) Time 120ms. (b) Time 180ms.

Fig. 15. Dipole pattern of BSS component corresponding to cSI and its transformation. (a) Topographic field map of dishonest BSS, (b) Topographic field map of honest BSS, (c) Transformation from the fourth BSS to components of honest BSS.

22 Will-be-set-by-IN-TECH

Fig. 14. Dipole patterns in the topographic field map of averaged waveforms of subject D in

0 10 20 30 40 50 60 70 −0.8

(c)

Fig. 15. Dipole pattern of BSS component corresponding to cSI and its transformation. (a) Topographic field map of dishonest BSS, (b) Topographic field map of honest BSS, (c)

P

(b)

A

A

P

P

(b) *Ah* 32.

 (a) Time 120ms. (b) Time 180ms.

A

P

P

(a) *A*4.

−0.6 −0.4 −0.2 0 0.2 0.4 0.6

Transformation from the fourth BSS to components of honest BSS.

L R

BFL

AL R

L R PA L R

A

Fig. 16. Impulse responses between cSI and cSII by using the dishonest BSS. (a) Impulse responses from cSI to cSII (*F*12). (b) Impulse responses from cSII to cSI (*F*21).

#### **5.3 SEF with steady state responses**

Although separate fluctuations of SEF in subject D could not be separated, the other three subjects could be separated, and the intracerebral communication between cSI and cSII was

Florin, E., Gross, J., Pfeifer, J., Fink, G. R. & Timmermann, L. (2010). The effect of filtering on

Hamada Y., Otsuka, S., Okamoto, T. & Suzuki, R. (2002). The profile of the recovery cycle in

Hämäläinen, M., Hari, R., Ilmoniemi, R. J., Knuutila, J. & Lounasmaa, O. V. (1993).

Hironaga N. & Ioannides, A. A. (2007). Localization of individual area neuronal activity,

Hui, H. B., Pantazis, D., Bressler, S. L. & Leahy, R. M. (2010). Identifying true cortical interactions in MEG using the nulling beamformer, NeuroImage, 49, 3161-3174. Kakigi R., Hoshiyama, M., Shimoto, M., Naka, D., Yamasaki, H., Watanabe, S., Xiang, J.,

Kishida, K., Kanemoto, S. & Sekiya, T. (1976). Reactor noise theory based on system size

Kishida, K. (1991). Contraction of information in systems far from equilibrium, J. Math. Phys.,

Kishida, K. (1994). A theory of reactor diagnosis in feedback systems, J. Nucl. Sci. Technol., 31,

Kishida, K. (1996). Contraction of information and its inverse problem in reactor system

Kishida, K. (1997a). Numerical study on identification of transfer functions in a feedback

Kishida, K. (1997b). Identification of transfer function in a feedback system and autoregressive

Kishida, K. (1999). Identification of transfer functions in a feedback system and stochastic

Kishida, K., Fukai, H., Hara, T. & Shinosaki, K. (2003). A new approach to blind system identification in MEG data. IEICE Trans. Fundamentals, E86-A, 611-619. Kishida, K. (2005). Blind identification of brain mechanism in MEG. Proc. IEEE International Symposium on Circuits and Systems (ISCAS), Kobe, May, 5694-5697. Kishida, K. (2008). Classification of activities related to 5Hz periodical median nerve stimuli by

Kishida K. (2009a). Evoked magnetic fields of magnetoencephalogrphy and their statistical

Kishida, K. (2009b). Dynamical activities of primary somatosensory cortices studied by

Kishida, K. (2011). Intracerebral communication in 2 Hz periodical median nerve stimuli. J.

system and model reduction, J. Nucl. Sci. Technol., 34, 1115-1120.

Theory and its Applications, Yokohama, Nov., 25-30.

magnetoencephalography. Phys. Rev. E, 80, 051906(1-13).

Japan Biomag. Bioelectromag. 24 (in Japanese), 212-213.

identification, in Advances in Nuclear Science and Technology, edited by Lewins and

analysis, Proc. 11th IFAC symposium on system identification (SYSID' 97), 3,

inverse problem. Proc. 31st ISCIE International Symposium on Stochastic Systems

using the temporal decorrelation method of BSS, in Biomagnetism - Interdisciplinary Research and Exploration, R. Kakigi, K Yokosawa and S. Kuriki Eds. Sapporo;

studies of the working human brain, Rev. Mod. Phys., 65, 413-497.

study, Clin. Neurophysiol., 113, 1787-1793.

magnetic fields, Prog. Neurobiology, 61, 495-523.

expansion, J. Nucl. Sci. Technol., 13, 19-29.

NeuroImage, 34, 1519-1534.

Becker, Plenum Press, 23, 1-68.

Hokkaido University Press, 124-126.

property, Phys. Rev. E, 79, 011922(1-7).

32, 92-98.

526-538.

1437-1442.

Granger causality based multivariate causality measures, NeuroImage, 50, 577-588.

human primary and secondary somatosensory cortex: a magnetoencephalography

Magnetoencephalography - theory, instrumentation, and applications to noninvasive

Maeda, K., Lam, K., Itomi, K. & Nakamura, A. (2000). The somatosensory evoked

identified. Our method is applicable to the normal case, in which there are line spectrum of fundamental frequency of SEF and higher harmonic ones. In the case of more than 20Hz ISI, steady state responses were observed for the somatosensory evoked potentials (Snyder, 1992; Tobimatsu et al., 1999). Their mechanism is different from transient responses of SEF. Supposing the normal scaling does not hold in the case, a new identification formalism is needed for nonlinear model, instead of linear feedback model in Eq. (17).

## **6. Conclusion**

Intracerebral communication between cSI and cSII of human brain in 2Hz median nerves stimuli was studied by MEG. The honest type of blind source separation was used for identification of feedback model between cSI and cSII, current dipole data of which were in a stationary process.

In the special case where there was the one-to-one correspondence between current dipoles and BSS components (Kishida, 2009b), fluctuations of SEF could be separated after selection of BSS components with nonzero waveforms in the dishonest BSS. In general, there are no one-to-one correspondences between them. However, fluctuations of SEF could be determined from the honest BSS components as in Eq. (22), when more than three components of the dishonest BSS with nonzero waveforms were transformed into the honest BSS components by using the transformation matrix. That is, the transformation matrix between components of two types of BSS, Eq. (11), was introduced for selecting the fluctuations of evoked magnetic fields from MEG.

The identification method based on feedback system theory applied to fluctuations of evoked magnetic fields demonstrated two transit modes between cSI and cSII in impulse responses of transfer functions: One was a slow mode of order 30ms and the other was rapid one within 5 ms. It should be noted that the transfer functions are not necessarily elementary but overall process between active regions. If fluctuations generated by more than three current dipoles, we must use a multi-feedback model (Kishida, 1999) to obtain elemental transfer functions. Furthermore, physiological meanings of two modes in impulse responses will be future problems.

#### **7. Acknowledgment**

The author would like to thank Prof. K. Shinosaki for providing MEG data. This research was supported by Grant-in-Aid for Scientific Research (No.22500254) of Japan society for the Promotion of Science.

#### **8. References**

24 Will-be-set-by-IN-TECH

identified. Our method is applicable to the normal case, in which there are line spectrum of fundamental frequency of SEF and higher harmonic ones. In the case of more than 20Hz ISI, steady state responses were observed for the somatosensory evoked potentials (Snyder, 1992; Tobimatsu et al., 1999). Their mechanism is different from transient responses of SEF. Supposing the normal scaling does not hold in the case, a new identification formalism is

Intracerebral communication between cSI and cSII of human brain in 2Hz median nerves stimuli was studied by MEG. The honest type of blind source separation was used for identification of feedback model between cSI and cSII, current dipole data of which were in a

In the special case where there was the one-to-one correspondence between current dipoles and BSS components (Kishida, 2009b), fluctuations of SEF could be separated after selection of BSS components with nonzero waveforms in the dishonest BSS. In general, there are no one-to-one correspondences between them. However, fluctuations of SEF could be determined from the honest BSS components as in Eq. (22), when more than three components of the dishonest BSS with nonzero waveforms were transformed into the honest BSS components by using the transformation matrix. That is, the transformation matrix between components of two types of BSS, Eq. (11), was introduced for selecting the fluctuations of

The identification method based on feedback system theory applied to fluctuations of evoked magnetic fields demonstrated two transit modes between cSI and cSII in impulse responses of transfer functions: One was a slow mode of order 30ms and the other was rapid one within 5 ms. It should be noted that the transfer functions are not necessarily elementary but overall process between active regions. If fluctuations generated by more than three current dipoles, we must use a multi-feedback model (Kishida, 1999) to obtain elemental transfer functions. Furthermore, physiological meanings of two modes in impulse responses will be

The author would like to thank Prof. K. Shinosaki for providing MEG data. This research was supported by Grant-in-Aid for Scientific Research (No.22500254) of Japan society for the

Akaike, H. & Nakagawa, T. (1988). Statistical Analysis and Control of Dynamic Systems, KTK Scientific Pub., originally published (in Japanese) by Saiensusya, 1972. Cantero, J.L., Atienza, M., Cruz-Vadell, A., Suarez-Gonzalez, A. & Gil-Neciga, E. (2009).

Increased synchronization and decreased neural complexity underlie thalamocortical oscillatory dynamics in mid cognitive impairment, NeuroImage, 46, 938-948. Cardoso J. F. & Souloumiac, A. (1996). Jacobi angles for simultaneous diagonalization, SIAM

needed for nonlinear model, instead of linear feedback model in Eq. (17).

**6. Conclusion**

stationary process.

future problems.

**7. Acknowledgment**

Promotion of Science.

J. Mat. Anal. Aplpl., 17, 161-163.

**8. References**

evoked magnetic fields from MEG.

**11**

*Thailand*

Wichian Sittiprapaporn

**Pre-Attentive Processing of Sound Duration**

The human central auditory system has a remarkable ability to establish memory traces for invariant features in the acoustic environments in order to correct the interpretation of natural acoustic sound heard. Even when no conscious attention is paid to the surrounding sounds, changes in their regularity can cause the listener to redirect his or her attention toward the sound heard (Tervaniemi *et al.,* 2001). When engaged in a conversation, listeners tune in to the relevant stream of speech and filter out irrelevant speech input that may be present in the same environment. Nonetheless, attention might be involuntarily diverted to meaningful items coming from an ignored stream, like in the well-known own-name effect (Moray, 1959). This brings up the question of to what extent speech is processed in the ignored streams. In the past decade, there have witnessed a resurgence in the electrophysiological literature of attempts to understand how the brain processes the speech signal (Kraus *et al.,* 1993, 1996; Molfese, 1985). One of the most used and well known paradigms in electrosphysiological research is the so-called oddball paradigm, in which typically two stimuli are presented, in random order. One of the stimuli occurs less frequently than the other and the subject is required to discriminate the infrequent stimulus (deviant, target or oddball) from the frequent one (standard). Two main types of ERPs have been described in the literature as a response to the detection of the deviant: P300 and the MMN (Aaltonen *et al.,* 1994; Kraus *et al.,* 1993, 1996). If the subject is required to respond overtly --- for example, by pressing a button – each time he/she detects the deviant, a positive wave peaking approximately 300 ms after deviant onset is elicited. This wave is called P300 and it is largest over electrode sites in normal adults. Such positivity is thought to reflect voluntary focused attention (context updating, response selection). However, if the subject is not required to respond overtly, and one subtracts the event-related potentials (ERPs) obtained in response to the standard, from the ERPs obtained for the deviant, so-called mismatch negativity (MMN) may be observed, usually peaking between 100 and 300 ms after stimulus onset depending on the characteristics of the difference between standard and deviant stimuli. This component is thought to reflect a pre-attentional detection of deviation, a mismatch between the deviant

**1. Introduction**

and the memory trace formed by the standard.

**Changes: Low Resolution Brain**

*Department of Educational Psychology and Guidance,*

**Electromagnetic Tomography Study**

*Faculty of Education, Mahasarakham University, Maha Sarakham,*


**18. Embolectomy**

Several interventional techniques including mechanical fragmentation, thrombus aspiration, and direct thrombolytic therapy are currently available for the treatment of acute PE. No clear guidelines can be provided regarding the use of interventional techniques for PE as there are no randomized controlled trials to assess these treatments. ACCP guidelines advise consider‐ ation of these interventions in the setting of massive PE with contraindications to systemic thrombolysis, failed thrombolysis, or likelihood of death before thrombolysis (Kearon 2012). Surgical embolectomy is infrequently required for management of acute PE. There are limited data guiding surgical interventions in this setting. Leacche et al (Leacche 2005) reported surgical outcomes in 47 consecutive patients with acute massive PE treated surgically. Indications for surgery included contraindications to thrombolysis (45%), failed medical treatment (10%) and right ventricular dysfunction (32%). Operative mortality and late mortality were 6% and 12% with the majority (83%) of the late mortality related to metastatic cancer. 86% and 83% of patients were alive 1 and 3 years after surgery, respectively (Leacche 2005). Specific instances when surgical intervention may be indicated are the presence of right atrial thrombus, paradoxical arterial embolism, or closure of a patent foramen ovale (Kearon 2012).

**20. Summary**

**Author details**

Jean M. Elwing1,2\* and Ralph J. Panos1,2

cal Center, Cincinnati, USA

Medicine, Cincinnati, USA

**References**

\*Address all correspondence to<EMAIL_ADDRESS>
pects. Blood. 2009 Mar 26;113(13):2878-87.

Acute pulmonary thromboembolism is the third most common acute cardiovascular disease (Giuntini 1995) and is associated with significant morbidity and mortality. Thromboembolic obstruction of the pulmonary vasculature frequently causes dyspnea, hypoxia, and chest pain. In severe cases, PE can cause an acute rise in pulmonary pressures precipitating RV dysfunc‐ tion and hemodynamic instability. The increase in hemodynamic pressures is due to both mechanical obstruction of the pulmonary vascular bed as well as the release of vasoconstricting mediators and derangement of neural regulation of vasomotor tone. Hypotension, RV dysfunction, elevated biomarkers (BNP, troponin I, troponin T) and ECG abnormalities with T-wave inversions are markers of poor prognosis. Rapid evaluation and risk stratification are necessary for effective treatment. Hemodynamically stable PE, i.e. low risk PE, is routinely managed with heparin therapy and clinical monitoring. Hemodynamically unstable PE, i.e. massive PE, is frequently managed with thrombolysis. Optimal management of hemodynam‐ ically stable PE with high risk features, i.e. submassive PE, is unclear at this time. Further study to determine best practice in this group is ongoing. Unfortunately, acute PE has both shortterm and long-term consequences. Additional evaluation and characterization of patients at

Acute Thromboembolic Pulmonary Hypertension

http://dx.doi.org/10.5772/56719

133

1 Pulmonary, Critical Care, and Sleep Medicine Division, Cincinnati Veterans Affairs Medi‐

2 Pulmonary, Critical Care, and Sleep Medicine Division, University of Cincinnati College of

[1] Adam SS, Key NS, Greenberg CS. D-dimer antigen: current concepts and future pros‐

[2] Anderson DR, Kahn SR, Rodger MA, Kovacs MJ, Morris T, Hirsch A, Lang E, Stiell I, Kovacs G, Dreyer J, Dennie C, Cartier Y, Barnes D, Burton E, Pleasance S, Skedgel C, O'Rouke K, Wells PS.Computed tomographic pulmonary angiography vs ventila‐ tion-perfusion lung scanning in patients with suspected pulmonary embolism: a

randomized controlled trial. JAMA. 2007 Dec 19;298(23):2743-53.

risk to develop long term complications of pulmonary emboli are needed.