Clinical Neurophysiology 110 (1999) 1585±1588

Neural source estimation from a time±frequency component of somatic
evoked high-frequency magnetic oscillations to posterior tibial nerve
stimulation
Kenji Sakuma a,*, Kensuke Sekihara b, Isao Hashimoto a, ca

Department of Integrative Physiology, National Institute for Physiological Sciences, Okazaki 444-8585, Japanb
JST Mind Articulation Project, Yushima 4-9-2, Bunkyo, Tokyo 113-0034, Japanc
Department of Psychophysiology, Tokyo Institute of Psychiatry, 2-1-8 Kasmikitazawa, Setagaya-ku, Tokyo 156-0057, JapanAccepted 15 April 1999

Abstract
Objective: High frequency oscillations (HFOs) evoked by posterior tibial nerve stimulation were recorded using magnetoencephalography (MEG). Time±frequency domain multiple signal classi®cation (TF-MUSIC) algorithm was applied, and the usefulness of this method
was demonstrated.
Methods: Ten normal subjects were studied. To localize sources for the HFOs of those somatosensory evoked ®elds, we applied two kindsof methods: the single moving dipole (SMD) method and the TF-MUSIC method. The SMD method was applied after digitally band-pass®ltering the somatosensory response with a bandwidth of 500±800 Hz. To estimate the locations of sources with the TF-MUSIC algorithm,
we ®rst set the target region on the spectrogram of the somatosensory responses. Then, the procedure described in Section 2.2 was appliedwith this target region.
Results: A clear, isolated region was detected in 6 out of 10 subjects using a time±frequency spectrogram. The averaged distance of thedipole sources between the HFOs and the underlying P37m using the TF-MUSIC algorithm was smaller than using the SMD method.
Conclusions: The TF-MUSIC algorithm is suitable for extracting a target response whose spectrum changes signi®cantly during theobservation. q 1999 Elsevier Science Ireland Ltd. All rights reserved.
Keywords: Somatosensory cortex; High-frequency magnetic oscillations; Posterior tibial nerve; Magnetoencephalography; Time±frequency components;
Multiple signal classi®cation algorithm

1. Introduction
High-frequency primary cortical oscillations characterized by a frequency range of 500±800 Hz superimposedon the ascending slope of the N20 primary response following stimulation of the median nerve have been reported
previously (Yamada et al., 1988; Curio et al., 1994; Hashimoto et al., 1996). We recorded high frequency oscillations(HFOs) evoked by posterior tibial nerve stimulation using
magnetoencephalography (MEG), and demonstrated theintracranial locations of the HFOs sources with respect tobrain anatomy (Sakuma and Hashimoto, 1999). Since theHFOs are of extremely low amplitude with low signal-tonoise ratio (S/N) compared with the underlying P37m, it
was dif®cult to estimate the location of HFOs source* Corresponding author. Tel.: 1 81-564-55-7768/7766; fax: 1 81-56452-7913.
E-mail address: sakuma@nips.ac.jp (K. Sakuma)

precisely. We applied a time±frequency multiple signal
classi®cation (TF-MUSIC) algorithm to HFOs, and could
demonstrate the superiority of this method compared to asingle moving dipole solution method (SMD).
2. Materials and methods
Ten normal subjects (5 females and 5 males; mean age
26.6 years, range 21±35 years) were studied. They gave
their consent to the experimental procedures. Electricalstimuli with 0.2 ms duration were delivered to the rightposterior tibial nerve at the ankle (cathode proximal). The
stimulus intensity was about 3 times the sensory thresholdand elicited a mild twitch of toe. The stimuli were deliveredat regular intervals with a repetition rate of 4 Hz. Magneticrecordings (bandpass 0.1±1200 Hz) were taken from thevertex centering around the Cz of the international 10-20
system with a 37-channel biomagnetometer (Magnes,1388-2457/99/$ - see front matter q 1999 Elsevier Science Ireland Ltd. All rights reserved.
PII: S 1388-245 7(99)00120-0

CLINPH 98179

1586

K. Sakuma et al. / Clinical Neurophysiology 110 (1999) 1585±1588

Biomagnetic Technologies Inc., San Diego, CA) in amagnetically shielded room.
An epoch of 60 ms duration (10 ms pre- and 50 ms poststimulus) was digitized at a 4167 Hz/channel sampling rate
and 9999 responses were averaged off-line. DC offset wasbased on the pre-stimulus period. A local spherical model
was ®tted to the digitized head shape over the recording areafor each subject.
2.1. SMD method
For separation of the HFOs from the underlying P37m,
the wide-band (0.1±1200 Hz) recorded responses were digitally ®ltered with a 500±800 Hz band-pass ®lter. The location, orientation, and amplitude parameters of a dipole
source model in the spherical head model were iteratively
adjusted to obtain the least-squares error ®t between anobserved ®eld and a theoretical ®eld produced by the dipole
source model (Sarvas, 1987). The correlation between thetheoretical ®eld generated by the model and the observed®eld was calculated and only equivalent current dipoles(ECDs) with a correlation above 0.86 were analyzed.
2.2. TF-MUSIC algorithm
The time±frequency MUSIC algorithm allows us to estimate the locations of neural sources from any time±
frequency region of interest (Sekihara et al., 1998a, b). In
this algorithm, we ®rst calculate Cb(t, f), the time±frequency
representation matrix for the measured data. Its diagonal
elements are the auto-time±frequency distributions of thechannel recordings and its off-diagonal elements are thecross-time-frequency distributions between different channel recordings. The matrix Cb(t, f) can be calculated using
both linear and quadratic time±frequency representations.
Let us de®ne the time±frequency representation matrix forthe source activities as Cs(t, f). Its diagonal elements are theauto-time±frequency distributions of the source activities
and its off-diagonal elements are the cross-time±frequency
distributions between different source activities.
Let us also de®ne the matrices Yb and Ys as those obtainedby averaging Cb(t, f) and Cs(t, f) over the target time±
are locafrequency region
RRV for which the neural sources
RRlized; i.e.
YC(t,f)
dtdf
andY
Ct;
f  dt df,
V bsV sRR bindicatestheintegral
overthetarget
time±
whereVfrequency region. Assuming that the noise and signal areuncorrelated in the target region and the noise is white Gaussian, we ®nally derive the relationshipY b  Lc Y s LTc  1 r2 I1

where Lc is the lead-®eld matrix for the true source locations,r 2 is the total noise power in the target region and I is the unitmatrix.
Let us denote the noise-level eigenvectors of Yb as uj (j PV 1 1,¼, M). Here PV is the number of sources whoseactivities have time±frequency components in the target

region. Therefore, when the matrix Ys is a full-rank matrix,
Eq. (1) leads toLTc uj  0forj  PV 1 1; ¼; M2

This equation indicates that source locations can beobtained by checking the orthogonality between the sensor
lead ®eld L(x) and the noise-level eigenvectors of Yb. Thisorthogonality can be evaluated by calculating the following
localizing function,hi3J x  1=lmin L xT ZN ZNT L x; L xT L x
where the matrix ZN is de®ned as ZN  [uPV11,¼, uM] andl min[.,.] indicates the generalized minimum eigenvalue ofthe matrix pair given in parenthesis. Eq. (3) is calculated in a
volume where sources can exist, and each location whereJ(x) reaches a peak is chosen as the location of one dipole
source.
2.3. Evaluation of localization accuracy
The distances between ECD location of P37m and HFOs
calculated with the SMD and between ECD location ofP37m and HFOs calculated with TF-MUSIC werecompared and the accuracy of ECD location between thetwo methods were assessed, since it is well established thatthe sources for HFOs and the underlying N20m for median
nerve stimulation are very close to each other (Hashimoto etal., 1996). These distances were de®ned as ((x 2 x 0 ) 2 1(y 2 y 0 ) 2 1 (z 2 z 0 ) 2) 1/2, where x, y and z indicate P37m
dipole locations and x 0 , y 0 and z 0 , HFOs dipole location.
3. Results
In wide-band recordings, the main P37m de¯ection wasclearly identi®ed in all 10 subjects. The ECDs of P37m peak
(SMD P37m) were estimated at the somatosensory foot areain all subjects. After ®ltering with a 500±800 Hz band-pass®lter, low amplitude high-frequency magnetic oscillations
superimposed on P37m were detected in 5 out of 10subjects. However, since the values of the correlation for
the ECD sources of the HFOs were smaller than 0.86 in one
subject, we calculated the ECDs of HFOs in only 4 subjects
using the SMD method. With the SMD method, the ECD
localizations for the high-frequency peaks and the underlying P37m were close with each other (Table 1). There
were no signi®cant statistical differences between the localizations for the x, y and z axes (paired t test) (Table 1). The
spectrogram was calculated for the frequency range above
400 Hz from the averaged wide-band somatosensory evokedmagnetic ®elds (SEFs). The region containing the HFOs
was selected as the target regions of interest. A clear andisolated time±frequency component was observable in 6 outof 10 subjects. The time±frequency region of interest started
at 29±33 ms (30:8 ^ 1:5ms) with 500±630 Hz (585 ^ 44:6Hz) in frequency, and ended at 35±40 ms (38:3 ^ 1:8 ms)

K. Sakuma et al. / Clinical Neurophysiology 110 (1999) 1585±1588

1587

Table 1
Individual source locations of P37m and HFOs and the distance between the locations of P37m and HFOs using the SMD and TF-MUSIC alogorithm aSubjectKTMKYI
KSAG
TO
Average
SDSMD P37m

SMD HFOs

TF-MUSIC HFOsDistance (cm)xyzxyzxyz

SMD

TF-MUSIC

1.52
1.090.90
1.92
20.22
0.17

1.41
1.89
1.35
1.19
1.27
1.36

10.25
10.29
10.73
10.66
10.23
10.95

N.E.
1.03
0.68
1.25
N.E.
20.27

N.E.
2.41
1.36
2.54
N.E.
0.54

N.E.
9.44
10.33
10.31
N.E.
9.54

0.80
2 0.46
1.65
3.11
0.49
0.18

1.57
1.52
2.26
1.95
2.22
1.38

10.14
10.40
10.50
9.57
10.02
10.65

N.E.
1.00
0.46
1.55
N.E.
1.70
1.17
0.57

0.74
1.59
1.20
1.78
1.21
0.30
1.14
0.55

a

Comparison of values between SMD P37m and SMD HFOs, not signi®cant. Comparison of values between SMD P37m and TF-MUSIC HFOs, notsigni®cant. Comparison of values between Distance SMD and Distance TF-MUSIC, not signi®cant. N.E., not estimated.

with 650±800 Hz (725 ^ 50:1 Hz) in frequency. These time
ranges were coincident with the ascending and descending
slopes of the P37m. Fig. 1A shows a typical record from one
subject where the ECDs of HFOs could be estimated byusing both the SMD and the TF-MUSIC. Fig. 1B shows atypical result from another subject showing clear isolated
region in the time±frequency spectrogram although theSMD algorithm failed to estimate the ECD localization
due to a low signal-to-noise ratio with a low correlation
value. We then applied the TF-MUSIC algorithm to thesetarget regions. Typical TF-MUSIC localizing function maps
are shown in Fig. 2. The contours in this ®gure show therelative values of the TF-MUSIC localization function (Eq.
(3)) projected onto the transverse, coronal and sagittal
planes. Although the average distance of the dipole sources
between the HFOs and the underlying P37m was not significantly different between SMD method and TF-MUSIC

algorithm, it should be noted that HFO dipole locations
were estimated in 6 out of 10 subjects with the TFMUSIC algorithm as against in 4 out of 10 subjects withthe SMD method (Table 1).
4. Discussion
In this study, the TF-MUSIC algorithm was applied toestimating the neural sources of the high-frequency
magnetic oscillations evoked by posterior tibial nerve
stimulation. The dipole locations of the HFOs were successfully estimated using the TF-MUSIC algorithm even whenthe estimations failed with the SMD method. This is mainlybecause only the high-frequency oscillatory activities visualized in a time±frequency spectrogram were used for calculation. The dipole strength of the HFOs is as small as oneFig. 1. (A) A typical recording from one subject where the ECD of HFOs could be estimated using both the SMD algorithm and the TF-MUSIC algorithm(subject TO). Wide-band recorded somatosensory evoked magnetic ®elds following right posterior tibial nerve stimulation. Thirty-seven waveforms aresuperimposed to illustrate the polarity-reversed P37m de¯ections (a). High-frequency oscillations after digital ®ltering (500±800 Hz) the wide-band somatosensory magnetic ®elds illustrated in (a). The high-frequency oscillations are roughly coincident with the P37m (b). Averaged spectrogram for the frequency
range above 400 Hz obtained from the MEG data is shown in (c). A rectangle with a dotted line is the target region for the HFOs. (B) A typical result fromanother subject showing a clear isolated region in time±frequency spectrogram although the ECD of HFOs could not be estimated with the SMD algorithm dueto a low signal-to-noise ratio with a low correlation value (correlation  0.74) (subject AG).

1588

K. Sakuma et al. / Clinical Neurophysiology 110 (1999) 1585±1588

Fig. 2. (A,B) TF-MUSIC localization function projected onto transverse, coronal and sagittal planes with the target region chosen in Fig. 1A and B,respectively.

fourteenth of the N20m for median nerve stimulation(Hashimoto et al., 1996). In addition, it is widely known
that somatosensory evoked potentials evoked by posterior
tibial nerve stimulation are more dif®cult to record thanthose evoked by median nerve stimulation. Although weused `massive' averaging over nearly 10 000 stimulus
presentations to obtain the desirable S/N, it is impossible
to eliminate all magnetic noise such as external ambient
noise, system noise and brain noise. Since such a magnetic
noise has no phase-locking to the stimulus, the clear isolated
activity detected in the averaged time±frequency spectrogram should be the signal free from the noise. In fact, theTF-MUSIC has theoretical as well as empirical advantages
for eliminating noise from the evoked magnetic ®elds (Sekihara et al., 1998a,b). With this algorithm, the time±
frequency spectrogram of evoked magnetic response canbe checked by visual inspection for inclusion in the analysis.
Therefore, the TF-MUSIC algorithm is suitable for extracting a target response whose spectrum changes signi®cantly
during the observation.

References
Curio G, Mackert BM, Burghoff M, Koetitz R, Abraham-Fuchs K, Harer
W. Localization of evoked neuromagnetic 600 Hz activity in the cerebral somatosensory system. Electroenceph clin Neurophysiol
1994;91:483±487.
Hashimoto I, Mashiko T, Imada T. Somatic evoked high-frequency
magnetic oscillations re¯ect activity of inhibitory interneurons in thehuman somatosensory cortex. Electroenceph clin Neurophysiol
1996;100:189±203.
Sakuma K, Hashimoto I. High-frequency magnetic oscillations evoked byposterior tibial nerve stimulation. NeuroReport 1999;:2.
Sarvas J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys Med Biol 1987;32:11±22.
Sekihara K, Nagarajan S, Poeppel D, Miyauchi S, Takino R, Fujimaki N,
Miyashita Y. Estimating neural sources from each time±frequency
component of MEG data. Neuroimage 1998a;7:S673.
Sekihara K, Nagarajan S, Poeppel D, Miyashita Y. Time±frequency MEGMUSIC algorithm. IEEE Trans Med Imaging 1998b; in press.
Yamada T, Kameyama S, Fuchigami Y, Nakazumi Y, Dickins QS, Kimura
J. Changes of short latency somatosensory evoked potential in sleep.
Electroenceph clin Neurophysiol 1988;70:126±136.