- Email: [email protected]
Application of multivariate autoregressive modeling for analyzing the interaction between EEG and EMG in humans
International Congress Series 1270 (2004) 249 – 253
www.ics-elsevier.com
Application of multivariate autoregressive modeling for analyzing the interaction between EEG and EMG in humans Tomohiro Shibata a,b,c,*, Yuichi Suhara a, Tatsuhide Oga d, Yoshino Ueki e, Tatsuya Mima e, Shin Ishii a,c a
Graduate School of Information Science, Nara Institute of Science and Technology, Nara 630-0192, Japan b ATR Computational Neuroscience Laboratories, Japan c CREST, Japan Science and Technology Agency, Japan d Department of Neurology, Graduate School of Medicine, Kyoto University, Japan e Human Brain Research Center, Graduate School of Medicine, Kyoto University, Japan
Abstract. Understanding the network of the human motor control system in noninvasive ways is beneficial not only for designing human interfaces, but also to clinical applications. This article presents applications of multivariate autoregression (MVAR) modeling for analyzing the interaction between EEG and EMG, which is challenging because of their different modalities. In contrast to previous research employing the MVAR modeling by means of frequency-domain analysis, our approach emphasizes time-domain analysis. We examined one normal subject and one mirrormovement (MM) patient. The task was a weak isotonic contraction of the right abductor pollicis brevis muscle in the normal subject, and the left extensor carpi radialis brevis in the MM patient. For each subject, three channels consisting of two EEG signals and one EMG signal were analyzed. The EMG signals were from the bilateral primary sensorimotor cortices. By using the Bayesian Information Criterion (BIC), and by choosing the appropriate data length, the model order was determined in a stable fashion. Our results provided plausible information on EEG – EMG networks: (1) Information-transmission-delay time that seems physiologically appropriate, and (2) relative contribution from the ipsi- and contralateral corticospinal pathway, which is opposite in the normal subject in comparison to MM patients. D 2004 Elsevier B.V. All rights reserved. Keywords: Multivariate autoregression; Time-domain analysis; EEG; EMG; BIC
1. Introduction To understand the mechanism of the brain and other systems in humans, we need to understand not only localized functions, but also effective connectivities in their networks. It has been reported that multivariate autoregressive (MVAR) modeling holds some * Corresponding author. Graduate School of Information Science, Nara Institute of Science and Technology, 8916-5 Takayama-cho, Nara 630-0192, Japan. Tel.: +81-743-72-5981; fax: +81-743-72-5989. E-mail address: [email protected] (T. Shibata). 0531-5131/ D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ics.2004.05.048
250
T. Shibata et al. / International Congress Series 1270 (2004) 249–253
potential for such a purpose [1– 3]. We have been investigating applications of the modeling to analyze the interaction between EEG and EMG, since understanding the network of the human motor control system in noninvasive ways is beneficial not only for designing human interfaces, but also to clinical applications. Moreover, it is challenging from the viewpoint of engineering because of their different modalities. Mima et al. [3] detected the information flow from the left primary sensorimotor cortices (C3) to the right abductor pollicis brevis muscle in normal humans by means of frequency-domain analysis with MVAR modeling. In this article, we also employ the MVAR modelling while emphasizing time-domain analysis. We used data from not only normal subjects, but also a mirror-movement (MM) patient [4]. The comparison of these two cases is expected to validate our approach, because it is known that the contralateral cortex controls limbs in normal humans, while it is opposite in MM patients. 2. Methods We studied one normal subject and one MM patient, both of whom were right-handed. The protocol was approved by the Institutional Review Board and subjects gave written informed consent for the experiment. The task was a weak isotonic contraction of the right abductor pollicis brevis muscle in the normal subject, and the left extensor carpi radialis brevis in the MM patient. Subjects were instructed not to touch their thumb to the other fingers to avoid unnecessary tactile afferent feedback. Their target hand was kept covered. EEG and EMG signals were recorded during the tonic contraction for 2 min. EEG signals were recorded with a linked-ear reference, and surface EMG was recorded using pairs of electrodes and rectified. EEG and EMG signals were amplified in a bandpass of 1 –200 Hz and digitized at 1000 Hz in the normal subject, and at 2000 Hz in the MM patient. Preliminary analysis using a Fast Fourier Transform (FFT) algorithm revealed that the EEG signal was significantly coherent with the EMG and was localized over the contralateral primary sensorimotor area: C3 and C4 in both subjects. Therefore, C3 and C4 were used for analysis with the MVAR model throughout this article. The EEG signals were Hjorth-transformed to increase the spatial resolution [5]. More detailed information on the data acquisition method is described elsewhere [3]. To determine the order of the MVAR model, we employed the Bayesian Information Criterion (BIC) [6] instead of the Akaike Information Criterion (AIC) [7], since we found that AIC did not work as described in the next section. The BIC is an efficient asymptotic approximation to a Bayesian model selection criterion, taking parameters of uncertainty into account. For the time-domain analysis, we employed a method for analyzing feedback systems [1], and investigated impulse responses of identified systems. The determination of the model order was significantly affected by data length, presumably due to the tradeoff between asymptoticity in estimation and stationarity in data. To cope with this tradeoff, we sought an appropriate block length which enabled a model order to be determined consistently over all blocks segmented from the data. 3. Results Fig. 1 presents typical AIC (left) and BIC (right) values as a function of the model order determined with the same data. As shown in this figure, the AIC was unable to find the
T. Shibata et al. / International Congress Series 1270 (2004) 249–253
251
Fig. 1. Typical samples of AIC (left) and BIC (right) as functions of the model order.
minimum value within the order range that we sought, whereas the BIC robustly found it. Therefore, the BIC was employed to acquire the following results. Table 1 shows the means and the standard deviations of the determined model order, in which the block length was varied. In the case of the normal subject, the model order was consistently determined over blocks when the block lengths were 2 and 5 s. Regarding the data for the MM subject, the model order was consistently determined over blocks when the block lengths were 40 and 50 s. Fig. 2 illustrates that, for both subjects, the peak of the impulse response from the left EEG was significantly higher than that from the right EEG (Mann –Whitney U-test; p < 0.01). The impulse responses consisted of a leading noisy part followed by a certain response which was seemingly generated by a lower order linear system. The leading noisy part revealed white (Ljung – Box; p < 0.01), and its time length was almost the same with its corresponding model order. The model order, such as 20.3 (15.5) in the normal (MM) subject, consistently determined over blocks, can be interpreted as the informationtransmission-delay time from the brain to the target muscle, since it is physiologically plausible. Thus, we chose a block length of 2 s for the normal subject, and of 40 s for the MM patient. Fig. 3 shows the mean coherence spectra between the EEGs and the EMG of both the normal (A) and the MM patient (B). Thin lines were computed by the FFT method and thick lines by the 3-channel MVAR model; each thin line and the corresponding thick line
Table 1 Mean and standard deviation (S.D.) of the determined model order (in ms) Normal
MM
Block length [s]
Model order [ms]
Block length [s]
Model order [ms]
1 2 3 4 5
18.7 F 4.44 20.3 F 0.468 30.0 F 7.60 34.9 F 4.00 36.0 F 0.175
10 20 30 40 50
5.27 F 1.50 8.09 F 0.465 14.6 F 1.99 15.5 F 0.00 15.5 F 0.00
252
T. Shibata et al. / International Congress Series 1270 (2004) 249–253
Fig. 2. The lower quartile, median and upper quartile values of the impulse responses; from the left EEG (A) and the right EEG (B) to the right EMG in the normal subject; from the left EEG (C) and the right EEG (D) to the left EMG in the MM patient.
were very similar. In the case of the normal subject, the left EEG had larger coherence than the right EEG, while this phenomenon was reversed in the MM patient. 4. Discussion By emphasizing the time-domain analysis, we have demonstrated that the MVAR modeling approach provided plausible information about EEG – EMG networks: (1) Information-transmission-delay time, and (2) relative contribution from the ipsi- and contralateral corticospinal pathways. It is noteworthy that the known ipsilateral corticospinal contribution in the MM patient could not be detected by the coherence spectra (Fig. 3C,D), but could be found using the impulse responses (Fig. 2C,D). The fact that the AIC did not work very well for our data indicates that the EEG –EMG network in our data cannot be completely modeled by the true MVAR process. Nevertheless, the following results support our approach: (1) the model order consistently determined by BIC over blocks with a specific block length seems physiologically appropriate, and (2) coherence spectra computed by the 3-channel MVAR model were
Fig. 3. Mean coherence spectra computed by the FFT method (thin line) and by the 3-channel MVAR model (thick line) between the EEGs and the EMG; from the left EEG (A) and the right EEG (B) to the right EMG in the normal subject; from the left EEG (C) and the right EEG (D) to the left EMG in the MM patient.
T. Shibata et al. / International Congress Series 1270 (2004) 249–253
253
very similar to the one by the FFT method. The 10-fold difference in the appropriate block lengths for the normal and the MM subjects (Table 1) seems to be improved in the data acquisition step, since very low coherence spectra were observed between the EEGs and the EMG in the MM patient (Fig. 3C,D). References [1] H. Akaike, On the use of a linear model for the identification of feedback systems, Ann. Inst. Stat. Math. 20 (1968) 425 – 439. [2] M.J. Kaminski, K.J. Blinowska, A new method of the description of the information flow in the structures, Biol. Cybern. 65 (1991) 203 – 210. [3] T. Mima, T. Matsuoka, M. Hallet, Information flow from the sensorimotor cortex to muscle in humans, Clin. Neurophysiol. 112 (2001) 122 – 126. [4] S.F. Farmer, et al., Abnormal cortex – muscle interactions in subjects with X-linked Kallman’s syndrome and mirror movement, Brain 127 (2004) 385 – 397. [5] B. Hjorth, Source derivation simplifies topographical EEG interpretation, Am. J. E.E.G. Technol. 20 (1980) 121 – 132. [6] G. Schwartz, Estimating the dimension of a model, Ann. Stat. 8 (1978) 147 – 164. [7] H. Akaike, A new look at the statistical model identification, IEEE Transaction on Automatic Control AC 19 (6) (1974) 716 – 723.
www.ics-elsevier.com
Application of multivariate autoregressive modeling for analyzing the interaction between EEG and EMG in humans Tomohiro Shibata a,b,c,*, Yuichi Suhara a, Tatsuhide Oga d, Yoshino Ueki e, Tatsuya Mima e, Shin Ishii a,c a
Graduate School of Information Science, Nara Institute of Science and Technology, Nara 630-0192, Japan b ATR Computational Neuroscience Laboratories, Japan c CREST, Japan Science and Technology Agency, Japan d Department of Neurology, Graduate School of Medicine, Kyoto University, Japan e Human Brain Research Center, Graduate School of Medicine, Kyoto University, Japan
Abstract. Understanding the network of the human motor control system in noninvasive ways is beneficial not only for designing human interfaces, but also to clinical applications. This article presents applications of multivariate autoregression (MVAR) modeling for analyzing the interaction between EEG and EMG, which is challenging because of their different modalities. In contrast to previous research employing the MVAR modeling by means of frequency-domain analysis, our approach emphasizes time-domain analysis. We examined one normal subject and one mirrormovement (MM) patient. The task was a weak isotonic contraction of the right abductor pollicis brevis muscle in the normal subject, and the left extensor carpi radialis brevis in the MM patient. For each subject, three channels consisting of two EEG signals and one EMG signal were analyzed. The EMG signals were from the bilateral primary sensorimotor cortices. By using the Bayesian Information Criterion (BIC), and by choosing the appropriate data length, the model order was determined in a stable fashion. Our results provided plausible information on EEG – EMG networks: (1) Information-transmission-delay time that seems physiologically appropriate, and (2) relative contribution from the ipsi- and contralateral corticospinal pathway, which is opposite in the normal subject in comparison to MM patients. D 2004 Elsevier B.V. All rights reserved. Keywords: Multivariate autoregression; Time-domain analysis; EEG; EMG; BIC
1. Introduction To understand the mechanism of the brain and other systems in humans, we need to understand not only localized functions, but also effective connectivities in their networks. It has been reported that multivariate autoregressive (MVAR) modeling holds some * Corresponding author. Graduate School of Information Science, Nara Institute of Science and Technology, 8916-5 Takayama-cho, Nara 630-0192, Japan. Tel.: +81-743-72-5981; fax: +81-743-72-5989. E-mail address: [email protected] (T. Shibata). 0531-5131/ D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ics.2004.05.048
250
T. Shibata et al. / International Congress Series 1270 (2004) 249–253
potential for such a purpose [1– 3]. We have been investigating applications of the modeling to analyze the interaction between EEG and EMG, since understanding the network of the human motor control system in noninvasive ways is beneficial not only for designing human interfaces, but also to clinical applications. Moreover, it is challenging from the viewpoint of engineering because of their different modalities. Mima et al. [3] detected the information flow from the left primary sensorimotor cortices (C3) to the right abductor pollicis brevis muscle in normal humans by means of frequency-domain analysis with MVAR modeling. In this article, we also employ the MVAR modelling while emphasizing time-domain analysis. We used data from not only normal subjects, but also a mirror-movement (MM) patient [4]. The comparison of these two cases is expected to validate our approach, because it is known that the contralateral cortex controls limbs in normal humans, while it is opposite in MM patients. 2. Methods We studied one normal subject and one MM patient, both of whom were right-handed. The protocol was approved by the Institutional Review Board and subjects gave written informed consent for the experiment. The task was a weak isotonic contraction of the right abductor pollicis brevis muscle in the normal subject, and the left extensor carpi radialis brevis in the MM patient. Subjects were instructed not to touch their thumb to the other fingers to avoid unnecessary tactile afferent feedback. Their target hand was kept covered. EEG and EMG signals were recorded during the tonic contraction for 2 min. EEG signals were recorded with a linked-ear reference, and surface EMG was recorded using pairs of electrodes and rectified. EEG and EMG signals were amplified in a bandpass of 1 –200 Hz and digitized at 1000 Hz in the normal subject, and at 2000 Hz in the MM patient. Preliminary analysis using a Fast Fourier Transform (FFT) algorithm revealed that the EEG signal was significantly coherent with the EMG and was localized over the contralateral primary sensorimotor area: C3 and C4 in both subjects. Therefore, C3 and C4 were used for analysis with the MVAR model throughout this article. The EEG signals were Hjorth-transformed to increase the spatial resolution [5]. More detailed information on the data acquisition method is described elsewhere [3]. To determine the order of the MVAR model, we employed the Bayesian Information Criterion (BIC) [6] instead of the Akaike Information Criterion (AIC) [7], since we found that AIC did not work as described in the next section. The BIC is an efficient asymptotic approximation to a Bayesian model selection criterion, taking parameters of uncertainty into account. For the time-domain analysis, we employed a method for analyzing feedback systems [1], and investigated impulse responses of identified systems. The determination of the model order was significantly affected by data length, presumably due to the tradeoff between asymptoticity in estimation and stationarity in data. To cope with this tradeoff, we sought an appropriate block length which enabled a model order to be determined consistently over all blocks segmented from the data. 3. Results Fig. 1 presents typical AIC (left) and BIC (right) values as a function of the model order determined with the same data. As shown in this figure, the AIC was unable to find the
T. Shibata et al. / International Congress Series 1270 (2004) 249–253
251
Fig. 1. Typical samples of AIC (left) and BIC (right) as functions of the model order.
minimum value within the order range that we sought, whereas the BIC robustly found it. Therefore, the BIC was employed to acquire the following results. Table 1 shows the means and the standard deviations of the determined model order, in which the block length was varied. In the case of the normal subject, the model order was consistently determined over blocks when the block lengths were 2 and 5 s. Regarding the data for the MM subject, the model order was consistently determined over blocks when the block lengths were 40 and 50 s. Fig. 2 illustrates that, for both subjects, the peak of the impulse response from the left EEG was significantly higher than that from the right EEG (Mann –Whitney U-test; p < 0.01). The impulse responses consisted of a leading noisy part followed by a certain response which was seemingly generated by a lower order linear system. The leading noisy part revealed white (Ljung – Box; p < 0.01), and its time length was almost the same with its corresponding model order. The model order, such as 20.3 (15.5) in the normal (MM) subject, consistently determined over blocks, can be interpreted as the informationtransmission-delay time from the brain to the target muscle, since it is physiologically plausible. Thus, we chose a block length of 2 s for the normal subject, and of 40 s for the MM patient. Fig. 3 shows the mean coherence spectra between the EEGs and the EMG of both the normal (A) and the MM patient (B). Thin lines were computed by the FFT method and thick lines by the 3-channel MVAR model; each thin line and the corresponding thick line
Table 1 Mean and standard deviation (S.D.) of the determined model order (in ms) Normal
MM
Block length [s]
Model order [ms]
Block length [s]
Model order [ms]
1 2 3 4 5
18.7 F 4.44 20.3 F 0.468 30.0 F 7.60 34.9 F 4.00 36.0 F 0.175
10 20 30 40 50
5.27 F 1.50 8.09 F 0.465 14.6 F 1.99 15.5 F 0.00 15.5 F 0.00
252
T. Shibata et al. / International Congress Series 1270 (2004) 249–253
Fig. 2. The lower quartile, median and upper quartile values of the impulse responses; from the left EEG (A) and the right EEG (B) to the right EMG in the normal subject; from the left EEG (C) and the right EEG (D) to the left EMG in the MM patient.
were very similar. In the case of the normal subject, the left EEG had larger coherence than the right EEG, while this phenomenon was reversed in the MM patient. 4. Discussion By emphasizing the time-domain analysis, we have demonstrated that the MVAR modeling approach provided plausible information about EEG – EMG networks: (1) Information-transmission-delay time, and (2) relative contribution from the ipsi- and contralateral corticospinal pathways. It is noteworthy that the known ipsilateral corticospinal contribution in the MM patient could not be detected by the coherence spectra (Fig. 3C,D), but could be found using the impulse responses (Fig. 2C,D). The fact that the AIC did not work very well for our data indicates that the EEG –EMG network in our data cannot be completely modeled by the true MVAR process. Nevertheless, the following results support our approach: (1) the model order consistently determined by BIC over blocks with a specific block length seems physiologically appropriate, and (2) coherence spectra computed by the 3-channel MVAR model were
Fig. 3. Mean coherence spectra computed by the FFT method (thin line) and by the 3-channel MVAR model (thick line) between the EEGs and the EMG; from the left EEG (A) and the right EEG (B) to the right EMG in the normal subject; from the left EEG (C) and the right EEG (D) to the left EMG in the MM patient.
T. Shibata et al. / International Congress Series 1270 (2004) 249–253
253
very similar to the one by the FFT method. The 10-fold difference in the appropriate block lengths for the normal and the MM subjects (Table 1) seems to be improved in the data acquisition step, since very low coherence spectra were observed between the EEGs and the EMG in the MM patient (Fig. 3C,D). References [1] H. Akaike, On the use of a linear model for the identification of feedback systems, Ann. Inst. Stat. Math. 20 (1968) 425 – 439. [2] M.J. Kaminski, K.J. Blinowska, A new method of the description of the information flow in the structures, Biol. Cybern. 65 (1991) 203 – 210. [3] T. Mima, T. Matsuoka, M. Hallet, Information flow from the sensorimotor cortex to muscle in humans, Clin. Neurophysiol. 112 (2001) 122 – 126. [4] S.F. Farmer, et al., Abnormal cortex – muscle interactions in subjects with X-linked Kallman’s syndrome and mirror movement, Brain 127 (2004) 385 – 397. [5] B. Hjorth, Source derivation simplifies topographical EEG interpretation, Am. J. E.E.G. Technol. 20 (1980) 121 – 132. [6] G. Schwartz, Estimating the dimension of a model, Ann. Stat. 8 (1978) 147 – 164. [7] H. Akaike, A new look at the statistical model identification, IEEE Transaction on Automatic Control AC 19 (6) (1974) 716 – 723.