-1-

Accepted manuscript

Toda A, Imamizu H, Kawato M, Sato MA. (2011). Reconstruction of two-dimensional
movement trajectories from selected magnetoencephalography cortical currents bycombined sparse Bayesian methods. NeuroImage, 54(2), 892-905.

Request for a reprint PDF to imamizu@gmail.com

-2-



Reconstruction of two-dimensional movement trajectories fromselected magnetoencephalography cortical currents bycombined sparse Bayesian methods
Akihiro Toda 1, 2, †, *, Hiroshi Imamizu 3, 4, *, Mitsuo Kawato 1 and Masa-aki Sato 51

Computational Neuroscience Laboratories, Advanced Telecommunications ResearchInstitute International, 2-2-2 Hikaridai, Keihanna Science City, Kyoto 169-0288,



Japan
2

Nagaoka University of Technology, 1603-1 Kamitomioka, Nagaoka, Niigata
940-2188, Japan

3

National Institute of Information and Communications Technology, 2-2-2 Hikaridai,
Keihanna Science City, Kyoto 169-0288, Japan



4

Cognitive Mechanisms Laboratories, Advanced Telecommunications ResearchInstitute International, 2-2-2 Hikaridai, Keihanna Science City, Kyoto 169-0288,
Japan

5


Neural Information Analysis Laboratories, Advanced Telecommunications ResearchInstitute International, 2-2-2 Hikaridai, Keihanna Science City, Kyoto 169-0288,
Japan† Current affiliation: Fundamental Technology Research Center, Honda Research andDevelopment Corporation
* These authors contributed equally to this work.

Address correspondence to Hiroshi Imamizu
2-2-2 Hikaridai, Keihanna Science City, Kyoto 619-0288, Japan

Telephone: +81 774 95 1220, Fax: +81 774 95 1236
e-mail: imamizu@atr.jp

-3-Reconstruction of movements from non-invasively recorded brain activity is a key
technology for brain-machine interfaces (BMIs). However, electroencephalography(EEG) or magnetoencephalography (MEG) inevitably records a mixture of signals
originating from many cortical regions, and thus it is not only less effective than

invasive methods but also poses more difficulty for incorporating neuroscienceknowledge. We combined two sparse Bayesian methods to overcome this difficulty.
First, thousands of cortical currents were estimated on the order of millimeters andmilliseconds by a hierarchical Bayesian MEG inverse method, and then asparse-regression method automatically selected only relevant cortical currents inaccurate reconstruction of movements by a linear weighted sum of their time-series.
Using the combined methods, we reconstructed two-dimensional trajectories of theindex fingertip during pointing movements to various directions by moving the wrist
joint. A good generalization (reconstruction) performance was observed for test
datasets: mean error between the predicted and actual positions was 15 mm, which was

7% of the path length of the required movement. The reconstruction accuracy of the
proposed method was significantly higher than directly using MEG sensor signals.
Moreover, spatial distribution and temporal characteristics of weight values revealed
that the primary sensorimotor, higher motor, and parietal regions mainly contributed to
the reconstruction with expected time-courses. These results suggest that the combinedsparse Bayesian methods provide effective means to predict movement trajectory fromnoninvasive brain activity directly related to sensorimotor control.

-4-Invasive recording techniques using implanted electrodes in the brain have succeeded inreconstructing hand motions from neuronal activity in the primary motor cortex(Hochberg et al., 2006; Nicolelis, 2001; Taylor et al., 2002) and in predicting intended
goals from activity in parietal regions (Musallam et al., 2004). However, non-invasiverecordings are desirable for broad use in BMIs due to their longevity, cost, and safetywith regard to surgery. Moreover, non-invasive recording methods having fine temporal
resolution, such as EEG and MEG, are needed for reconstruction of rapid and smooth
movements. Previous studies using EEG or MEG sensor signals have succeeded incontrolling a computer cursor. For example, it has been demonstrated that subjects canlearn to modulate sensorimotor (mu and beta) rhythms in EEG signals and control acursor (Wolpaw and McFarland, 2004). A classification of patterns of EEG sensor
signals has been used for controlling a cursor in video games (Krepki et al., 2006).
Hand-movement direction (left, right, up or down) can be decoded during a reach by aclassification of MEG and EEG sensor signals above the motor area (Waldert et al.,2008). Reaching targets can be predicted during the planning period from a combinationof many features, such as the power levels of several frequency bands extracted from
EEG sensor signals (Hammon et al., 2007). The trajectories of a joystick position(Georgopoulos et al., 2005) and hand velocity (Bradberry et al., 2010; Bradberry et al.,2009) were reconstructed from the weighted sum of MEG or EEG sensor signals.

However, reconstruction of movements using EEG or MEG sensor signals is
generally thought to be inferior to reconstruction using invasive recordings for the

-5-following reasons. First, the spatial resolution and quality of data obtained by tens or
hundreds of sensors are contaminated by noise and poor in comparison to the data froma multiple-microelectrode array that records neuronal activity more directly than
non-invasive sensors. Second, sensors inevitably record a mixture of signals originatingfrom many cortical regions, layers and neurons. Therefore, neural signals of interest fordecoding not only make up a small portion of sensor signals but are also contaminatedby those from other regions, layers, and neurons of no interest. Third, a successful BMI,
regardless of whether it’s invasive or non-invasive, utilizes signals from restrictedcortical regions that are known to contribute to the functions to be reconstructed

according to previous studies in neuroscience. That is, information related to visualfunctions has been decoded from visual areas (Kamitani and Tong, 2005; Miyawaki etal., 2008), and signals related to motor control have been measured in the primary motor
cortex (Hochberg et al., 2006; Nicolelis, 2001; Taylor et al., 2002). Therefore, it isreasonable to expect increased reconstruction accuracy by using neuronal signals



precisely extracted from target cortical regions. In functional imaging studies, MEG or
EEG source currents have been estimated from sensor data to extract neuronal signals
originating from specific cortical regions. For example, a study (Jerbi et al., 2007)
estimated MEG source currents related to hand movements and found significant phaselocking between the tangential velocity of the hand during a tracking task and

oscillatory activity of the source currents in a specific region of the primary motor
cortex.

-6-

To make the effectiveness of non-invasive BMI closer to that of invasive BMI, we
propose a method for reconstruction of rapid two-dimensional movements (duration
about 0.4 s) from source currents estimated on cortical surfaces. Thousands of corticalcurrents were estimated on the order of millimeters and milliseconds by a hierarchical

Bayesian method that solves an inverse problem (projection from sensors to current
sources) by incorporating functional magnetic resonance imaging (fMRI) activity as ahierarchical prior (Sato et al., 2004; Yoshioka et al., 2008). The merits of thereconstruction using cortical source currents are as follows. First, improvement of
spatial resolution from the number of MEG sensors (several hundred) to that of currentsources (several thousand) is expected to yield rich information on the patterns ofcortical activity. Furthermore, data quality can be improved because it has been
demonstrated that simultaneous estimation of cortical and artifact currents is effective
for removing the artifacts whose source locations can be easily identified, such as eye
and cardiac movements and muscle activity (Fujiwara et al., 2009; Morishige et al.,2009). Second, isolation of cortical currents in a particular region from a mixture ofsignals originating from many regions would help in the extraction of crucial
information for reconstruction. Because different cortical regions are specialized indifferent functions, especially for primary sensory and motor functions, we assume thatit’s possible to achieve a better BMI for a particular brain function by using signals thatare not contaminated by signals from other cortical regions. Third, we can utilize
previous neuroscience knowledge if currents are mapped onto a cortical surface:
Currents relevant to a particular function can be selected based on knowledge about the

-7-relationships between brain functions and regions.
Selection of input signals relevant to useful features/parameters is important forreconstruction using decoding methods or machine-learning techniques. This is because
having too many parameters (e.g., weights for input signals in a multiple linear

regression model) in relation to the number of training datasets is known to lead to poor
generalization performance (over-fitting problem) (Akaike, 1974; Geman et al., 1992).
Recent advanced MEG systems have several hundreds of sensors, and each sensor has
high temporal resolution. Consequently, if we use signal time courses as input signals toa reconstruction model, the number of parameters becomes huge. By contrast, thenumber of training datasets is limited by relatively short time periods (up to several
hours), during which stable signals can be obtained in a comfortable situation for
subjects. Therefore, selection of useful input signals or parameters is essential for
ensuring high generalization performance of the model. We can obtain little benefit
from increasing the resolution of signals if we apply simple maximum likelihoodalgorithms due to the over-fitting problem, but it has been suggested that Bayesianalgorithms such as a sparse regression can effectively and automatically selectappropriate feature sets from thousands of parameters (Nambu et al., 2009; Sato, 2001;
Sato et al., 2004; Ting et al., 2005; Ting et al., 2008). This automatic selection should
not be applied to already mixed sensor signals but to cortical currents that are notcontaminated by signals from other cortical regions. Furthermore, we can apply
knowledge from neuroscience, which could be regarded as a qualitative prior, to the-8-selection of input signals if they are mapped onto the cortical surface.
In this paper, to reconstruct movements from non-invasively measured brain
activity, we combined two sparse algorithms: 1) a hierarchical Bayesian method that
calculates thousands of cortical currents based on sparse estimation of hyperparametersin an inverse filter from MEG signals to cortical currents, and 2) a sparse linear
regression that automatically selects effective cortical currents for the reconstruction.
Our results indicate not only that this combination improves reconstruction accuracy
compared to direct reconstruction from MEG sensor signals but also that we can gain
great advantages in using cortical currents, i.e., automatic selection of effective input



signals and incorporation of neuroscience knowledge in the selection on demand.
Furthermore, we could confirm which region and which time point were important for
the reconstruction by investigating the automatically selected currents on the corticalsurface.

Materials and methods



Subjects
Five male right-handed subjects (21–45 years of age) participated in this study. A
signed informed consent form, approved by the institutional ethics committee, wasobtained from each participant.

Task for subjects

-9-

Subjects moved the right index fingertip from a start position in eight directions
separated by 45º (gray arrows in Fig. 1A) while lying in the supine position in MEG and
MR scanners (Fig. 1B). The subjects' forearms were fixed to a platform, and their finger
joints were immobilized by a brace (Fig. 1C) so that only wrist-joint movements were

allowed.
During the MEG experiment (Fig. 1D), eight light-gray lines arrayed radiallyaround the center of a screen were always visible as a reference of the movement
directions. A small square was displayed at the center as a fixation point at which
subjects were instructed to fix their eyes. Subjects immobilized the hand while thesquare was black. This fixation point changed to white 1.6 s after the beginning of thetrial, and subjects moved their index fingers along one of the reference lines, which
subjects freely chose, and back to the start position until the color became black 0.4 s
later. Because these intervals were fixed, subjects could predict when they should startand end their movements. A motion-tracking system (see Data acquisition) recorded thefingertip position. The position was projected to the plane orthogonal to the body axis
(x-y plane, see Fig. 1B). No visual feedback was given to subjects during the movement,
but their fingertip paths were displayed for 0.5 s as black lines on the screen after the
movement to inform subjects of their movement accuracy. We analyzed MEG data from
1.0 s before to 1.0 s after the movement onset, during which the movement path was notdisplayed. At the end of each trial, a white peripheral square marked the direction inwhich subjects had moved. It remained in the subsequent trials, and subjects were
instructed to move the finger in an unmarked direction. The squares disappeared after

-10-

the eight directions were marked. This marking equalized the number of movements foreach direction. Trails were repeated in a session of 10 minutes, and subjects performed
six sessions. Eye blinks and movements were allowed only during an inter-trial intervalof 0.4 s. Electro-oculogram (EOG) was recorded to detect blinks and eye movements.

fMRI sessions were comprised of alternating blocks of execution and observation
periods (17.5 s each; Fig. 1E). During an execution period, subjects conducted pointing
movements in the direction indicated by a white peripheral square. They were asked tomove their fingers in the same fashion as in the MEG experiment. The path of thefingertip was displayed for 0.5 s on a screen 1.0 s after the end of each movement.



Subjects began the next movement as soon as the path disappeared from the screen.
During an observation period, screen images recorded in the preceding execution period
were replayed, and subjects observed them while immobilizing their hands. Eachsubject underwent four sessions. The order of movement direction was
pseudo-randomized within each session. There were eight execution and observation



periods for each direction (2 periods in a session x 4 sessions) across the fMRI
experiment.

Data acquisition
A motion-tracking system (QuickMag4 type 2; OKK Inc., Japan) was used to measure
finger movements. A marker of the system was attached to the right index fingertip, and

its position was recorded at 60 Hz. A starting position was registered in computer
memory at the beginning of MEG and fMRI experiments when the subject was asked to-11-keep his index finger parallel to his body axis as shown in Figure 1B.
A whole-head 208-channel system (MEG vision- PQ1400RM; Yokogawa ElectricCo., Japan) was used for MEG recording. The sampling frequency was 2 kHz. EOG andelectrocardiogram (ECG) were simultaneously recorded. Before the MEG experiment,the subject's face and head shape were scanned using a hand held laser scanner and astylus marker (FastSCAN Cobra; Polhemus, U.S.A) for later co-registration of MEG
and MRI results. To measure the head position in the MEG sensor system, four
calibration coils were bilaterally mounted on the subject’s temporal skin (two each forthe superior superciliary and anterior subauricular regions). Electromagnetic calibrationof the coil positions was conducted before and after each MEG recording session bypassing alternating currents to the coils.
A 1.5 Tesla MR scanner (MAGNEX ECLIPSE; Shimadzu-Marconi, Japan) wasused to obtain blood oxygen level-dependent contrast functional images. Images
weighted with the apparent transverse relaxation time were obtained with an echo

planer imaging sequence (repetition time, 3.5 s; echo time, 65 ms; flip angle, 90º). The
entire brain was covered in 44 axial slices (3.4-mm thickness; 1-mm gap), each ofwhich was acquired as a 64 x 64 matrix (field of view, 217.6 mm) with a voxel size of
3.4 x 3.4 x 4.4 mm. In total, 164 volumes were acquired in each session. T1-weighted
structural images were acquired with 1 x 1 x 1 mm resolution with a gradient echosequence.

Behavioral data analysis

-12-

The time series of the fingertip position was low-pass filtered with a cutoff frequency of
15 Hz. Data were recorded at 60 Hz, which was the highest sampling rate of the
tracking system, during the experiment. We increased the rate to 200 Hz using splineinterpolation to equalize the rate to that of the re-sampled MEG data (see below).
Movement initiation was defined as the first time the tangential velocity of the fingertip
crossed 5% of the maximum velocity of each trial, and movement termination was thelast time the velocity fell below 5% of the maximum. Mean movement time (interval
between movement initiation and termination) across trials was 406 ms (SD: 61), 486
ms (69), 540 ms (78), 559 ms (62) and 505 ms (47) for each subject. To investigate theeffect of movement direction on movement time, we applied a one-way (target
direction) analysis of variance (ANOVA) to movement time separately for each subject.
The effect was significant in every subject (F(7, 481) = 15.8, P < 0.0001 for the mostsignificant subject). However, according to Tukey’s HSD post-hoc test for multiple
comparisons, movement time for each target is significantly different from that for only

one or two other targets (1.7 targets averaged across targets and subjects, SD: 0.3) at P
< 0.05 level. This suggests relatively uniform distribution of movement time amongtarget directions. The mean movement time for each target across subjects was 497 ms(SD: 88) for the top, 495 ms (96) for the upper-right, 515 ms (109) for the right, 545 ms(130) for the lower-right, 496 ms (113) for the bottom, 459 ms (100) for the lower-left,



488 ms (131) for the left, and 496 ms (110) for the upper-left target.

MEG data preprocessing
We confirmed that head positions did not move more than 4 mm during any session bychecking the position data of calibration coils. Trials were excluded from analysis if thesignal value was larger than 1000 fT or if the EOG-signal value was larger than 20 !A.

-13-

If the MEG sensor's signal value exceeded 1000 fT in more than 10% of total trials, thesensor was considered impaired, and all of the data obtained by that sensor were also
excluded from analysis in each subject. Consequently, when averaged across subjects,
the total number of trials was 561 (SD: 68.3), the number of trials for each movementdirection was 70.2 (SD 3.31), and the number of effective sensors was 186 (SD: 15.5).
MEG-signals were passed through a low-pass filter with a cut-off frequency of 100 Hz
and sampled at 200 Hz. For each trial, the signal value was adjusted so that the meanvalue from 1.0 to 0.5 s before the movement onset became zero. A linear trend wasremoved by the least-squares fit of a straight line to the signal time course andsubtraction of the resulting function from the time course.
fMRI data analysis
Functional imaging data were analyzed using SPM2 (Wellcome Department of
Cognitive Neurology, London, UK; http://www.fil.ion.ucl.ac.uk/spm). We discarded thefirst four volumes of images in each session to allow for T1 equilibration and then

spatially aligned the data to the first remaining volume. The data were spatiallynormalized to the Montreal Neurological Institute (MNI; Montreal, Quebec, Canada)
reference brain and resliced to a 2-mm isotropic voxel size. Data were smoothed
spatially with a Gaussian kernel of 6 mm full-width at half-maximum (FWHM). Voxeltime series were high-pass filtered with a cutoff frequency of 0.002 Hz and low-passfiltered with a cutoff frequency of 0.25 Hz.
Statistical analyses were performed for each subject. Boxcar functions modeled
execution periods and observation periods. They were convolved with the canonicalhemodynamic response function in SPM2 to yield regressors in a general linear model.

-14-

A parameter was estimated for each regressor by the least-squares method. T-statistics
were used for comparison between the estimated parameters (execution - observation)
to yield a t-value for each voxel. Although subjects moved in a fixed direction in each
execution period (Fig. 1E), we did not distinguish movement directions in the analysis.


We used a threshold of P < 0.001 (uncorrected for multiple comparisons) according to
the previous study (Yoshioka et al., 2008) using hierarchical Bayesian estimation. The
yielded statistical parametric map was used as prior information in the estimation of
MEG source currents.

Cortical current estimation by the hierarchical Bayesian methodA polygon model of the cortical surface was constructed based on MR structural images
using Brain Voyager software (Brain Innovation, Maastricht, the Netherlands). For eachsubject, by using the hierarchical Bayesian method, we estimated about 2,500 (Mean:
2,579, SD: 210) single-current dipoles that were equidistantly distributed on andperpendicular to the cortical surface. The method calculated an inverse filter to estimate

the cortical current for each dipole from MEG sensor signals (see Supplementary
Information 1 for details of estimation). fMRI information was imposed on the priorinformation for the estimation through two types of parameters: a variancemagnification parameter ( m0 ) controlling the relative amplitude of the prior current
variance and a confidence parameter ( ! 0 ) controlling the width of the prior distribution.



The values of m0 and ! 0 were set at 100 and 10, respectively, in the current study.
The inverse filter was estimated by using the data of all trials, and the filter was applied
to sensor signals in each trial to calculate cortical currents.
A spatial smoothness constraint on the current distribution, along with the cortical

-15-

surface, was incorporated in the estimation (6 mm full-width at half-maximum, see
Supplementary Information 1). We also incorporated artifact dipoles in the estimationaccording to previous studies (Fujiwara et al., 2009; Morishige et al., 2009). Artifact
dipoles were located at the center of the heart, the right shoulder and wrist joints, theleft and right eyeballs, and the carotid arteries.
We assumed that the pattern of cortical activity changes according to various
phases of movements (e.g., planning, execution and feedback). Therefore, we divided
the MEG-signal time series into 39 time windows (100 ms length with 50 ms overlap)
from 1 s before to 1 s after the movement initiation and calculated an inverse filterseparately for each time window. Cortical currents were estimated every 5 ms (200 Hz)
from MEG data using the filter. In the overlap periods, they are averaged between twotime windows. Because the time series from 1.0 s to 0.5 s before the movement
initiation was used to estimate a baseline of the current variance (see Supplementary
Information 1), source currents during this period were not estimated. Therefore, thetime series of the estimated source currents ranged from 0.5 s before to 1.0 s after the
initiation.
Examination of time courses of estimated cortical currents
To validate the cortical currents estimated by the hierarchical Bayesian methods, we
compared time courses of the currents with those of MEG source current estimated by

previous methods (Cheyne et al., 2006; Huang et al., 2004) or electrocorticogram(ECoG) measurements (Kunieda et al., 2000) related to finger movements. Thesestudies have examined time courses mainly in the primary motor region (M1), the-16-

primary somatosensory region (S1), the supplementary motor area (SMA), and thepremotor region (PM) of the hemisphere contralateral to the moving hand; therefore, we
set our regions of interest (ROIs) in these regions using the Automated Anatomical
Labeling (AAL) map (Tzourio-Mazoyer et al., 2002) in the WFU PickAtlas (Maldjian


et al., 2003) (http://fmri.wfubmc.edu/cms/software) as shown in Figure 3B. We
investigated the peak latencies of the time courses averaged within the individual ROIs
and compared them to those reported in the previous studies.

Estimation of reconstruction model of fingertip position from cortical source currents
We conducted a multivariate linear regression in which the time course of the aboveestimated source currents were independent variables and the corresponding timecourses of the x (horizontal) and y (vertical) coordinates of the fingertip position weredependent variables:x (t ) =y (t ) =N source 10# # wijh ! Ji (t " j ) + w0h
i=1

j =0

,

(1)

,

(2)

N source 10# # wijv ! Ji (t " j ) + w0v
i=1

j =0

where x(t) and y(t) are fingertip positions at time t (from 0.3 s before to 1.0 s after
the movement initiation at regular intervals of 20 ms). Ji is the value of current
estimated at the i-th current source on the cortical surface. N source is the number ofcurrent sources used for reconstruction. We used the currents of the 1,500 sources withthe highest current amplitudes among the 2,500 sources mentioned above. The time

series of the currents ranging from 0.5 s before to 1.0 s after the movement initiation

-17-

were divided into 20-ms time windows and averaged within each window. They were
sampled at 11 time points within 200 ms preceding the currently predicted fingertip
position at t (see Fig. 3C). The estimation was done separately for the x and y positions.
Each of the regression models (Eqs. 1 and 2) has numerous weights (1,500 currentsources x 11 time points) to be estimated, and thus was expected to encounter anover-fitting problem. We examined the effectiveness of the following three methods incontrolling this problem.
A regularized least-squares regression is a technique often used to control theover-fitting problem. We added a regularization (penalty) term to an error function inorder to prevent the weights from reaching large values. This is done to find a set ofweights that minimizes an evaluation function for the x position, for example:
2
1 U T# N source 102( xu (t) ! x̂u (t)) + " " wij .
""2 u =1 t =1
2 i =1 j = 0

( )

(3)

The first term is the error function, in which xu (t) and x̂u (t) are the actual and thereconstructed positions, respectively, at a sampling time t in the u-th trial. U is thenumber of trials in training datasets (see below). T is the number of sampling points ineach trial. The coefficient ! governs the relative importance of the second
regularization term compared with the error term. We estimated reconstruction models
by varying ! at five values (1, 10, 100, 1000 and 10000).
The second method is a sparse regression that has been proposed and adopted inprevious studies (Nambu et al., 2009; Sato, 2001; Sato et al., 2004; Ting et al., 2005;
Ting et al., 2008). This method estimates the weight and the automatic relevancedetermination (ARD) parameters, which represent how the weight contributes to thereconstruction. Based on values of the ARD parameters, a weight value with a small

-18-

degree of contribution is set to nearly zero, and thus ineffective parameters are pruned(see Supplementary Information 2). In our post-hoc analysis, we investigated thenumber of ineffective sources that were pruned by the sparse regression: We calculated
the summation of weight values over time points for each source and counted thenumber of sources whose summation value was lower than 0.1% of the maximum value.
As a result, the number of ineffective sources was 1292.34 (SD: 5.25) averaged acrossdatasets and subjects, suggesting that a small set of weights ((1,500 – 1,292) / 1,500 =
13.9%) plays a crucial role in the reconstruction.
The above two methods automatically decrease weight values using the presentdata and reduce input dimensions. However, prior knowledge can often be used to selectan important set of dimensions in practical problems. Because current sources were
estimated on the cortical surface in our study, we could anatomically select sources inregions that have been known to contribute to motor control according to knowledge inneuroscience. That is, we determined an ROI consisting of SMA, PM, M1, S1 andparietal regions including Brodmann areas 5, 7, 39 and 40 based on the WFU PickAtlasand then estimated reconstruction models (Eqs. 1 and 2) using only currents fromsources in the ROI. The number of sources ( N source ) in the ROI was 287.4 (SD: 35.0)
averaged across subjects.

Evaluation of generalization ability of estimated models


We estimated reconstruction models combining (or not combining) the anatomical
selection of input currents with the regularized least-squares regression or the sparse
regression and evaluated their generalization ability to test datasets. We randomly
assigned sets of source currents and fingertip position data obtained in individual trials

-19-

to training and test datasets. The ratio of the number of trials in the training dataset tothat in the test dataset was 5:95, 20:80 or 70:30. That is, the percentage of the number oftrials in the training dataset to the total data was 5%, 20% or 70%. We made 50 pairs oftest and training datasets for each ratio and carried out the following procedure for each

pair. First, we estimated weights and constants by applying the above regression
analysis to the training dataset. Then we reconstructed the time course of a fingertip
position from source currents in the test dataset and calculated a positional error (E)between the reconstructed and actual positions in the u-th trial according toT

Eu = "( xu (t) ! x̂u (t))2 + ( yu (t) ! ŷu (t))2
T

t =1



.

(4)

Here, ( xu ( t ) , yu ( t )) is the actual position measured by the motion-tracking system at asampling time t, and ( x̂u ( t ) , ŷu ( t )) is the predicted position of the fingertip by themodel. We compared errors averaged across trials and datasets between different
reconstruction models to evaluate their generalization abilities.
Estimation of reconstruction model of fingertip position from MEG sensor signals

We also made a model that reconstructs the position directly from MEG sensor signals.
After removing cardiac artifacts by principal component analysis, we conducted asparse linear regression analysis in the same fashion as described above (Eqs. 1 and 2),
except that N source corresponded to the number of sensors. We compared theirperformance with that of reconstruction models using cortical source currents in termsof a positional error (Eq. 4) and a coefficient of determination. The coefficient was
separately calculated for the horizontal( R ) and vertical ( R ) directions:
h 2
i

v 2
i

-20-

T

2

( )

" t =1( xi (t) ! x̂i (t))
T2
" t =1( xi (t) ! xi ) ,( )" t =1( yi (t) ! ŷi (t))
T2
" t =1( yi (t) ! yi ) .

2
Rih = 1 !

2
Riv = 1 !

T

(5)

2

(6)

Here, ( xi , yi ) is mean position during the i-th trial. However, when the subjects moved
their fingertips toward the top or the bottom target (Fig. 1A), there was little positional
change in the horizontal direction. Therefore, the denominator in the right side ofEquation 5 becomes nearly zero, and thus we cannot calculate a reliable coefficientvalue. Consequently, we excluded from this analysis the coefficients in the horizontal
direction for a trial in which the subject aimed at the top or the bottom target. Similarly,
we excluded the coefficients in the vertical direction (Eq. 6) for a trial in which thesubject aimed at the rightmost or the leftmost target.

Anatomical localization of source currents contributing to reconstruction
We examined which anatomical region contributed to the reconstruction of the fingertip
position by investigating the weight values in the reconstruction model. Weight values
were averaged across time points and datasets and projected on an inflated model of the

cortical surface (Fig. 7A). Inward currents on the cortical surface corresponded topositive values of the estimated source currents, and signed values of weights for
current sources were represented by a pseudo-color.
For quantitative examination, we calculated the summation of absolute weightvalues separately for the following anatomical regions: 1) the frontal region includingareas anterior to the precentral sulcus (not including SMA), 2) the central region
including areas posterior to the precentral sulcus and anterior to the postcentral sulcus

-21-

and SMA, 3) the parietal region including Brodmann areas 5, 7, 39 and 40, 4) the
occipital region (areas 17, 18 and 19), and 5) the temporal region including the polar
gyrus, the superior temporal gyrus, the middle temporal gyrus and the inferior temporalgyrus (see bottom panels of Fig. 7C).


We also investigated the weight values estimated in a regression analysis of thefingertip position and MEG sensor signal. We averaged the values across time points
and datasets. Using a pseudo-color, we represented a topographic pattern of theaveraged weight values in the sensor space (Fig. 7B).
Investigation of temporal change of weight values in reconstruction model

We investigated temporal change of weight values for individual current sources in
representative ROIs (SMA, PM, M1, S1 and parietal regions; see above) as a functionof advanced time from the time point when a fingertip position is reconstructed. Since
we used source currents sampled at 11 time points from -200 to 0 ms for reconstruction(Fig. 3C), there are 11 weights in the temporal dimension for each current source. Time

courses of weight values were averaged across sources within an ROI, and then we
fitted an exponential curve to the time course:

w = a + b ! exp(t / " ) .

(7)

Here, w and t correspond to weight value and time, respectively. Parameters a, b and! were estimated by a least-squares method: a corresponds to the asymptotic level
(weight value at infinitely advanced time), b corresponds to an increase in value at 0 ms
from the asymptotic level, ! is known as a time constant, and - ! corresponds to time
when an increase in the weight value reaches 37% of b-value. Thus, the larger ! -value
indicates an earlier increase in the weight value as illustrated by the bottom right panel

-22-

in Figure 8A. We examined ! -values to investigate the temporal difference in each
ROI’s contribution to reconstruction.

Examination of contribution of movement-evoked magnetic fields to reconstruction
We examined whether signals contributing to the reconstruction of fingertip position

originated from a movement-evoked magnetic field (MEF). We averaged the MEG
signal separately for each sensor at 50 ms after the movement initiation across trials andsubjects, and then we created a template of the MEF’s topographic pattern by rendering
the averaged signal values in sensor space (color-coded contour plot in Fig. 9A). Basedon this template, we divided MEG sensors into two groups: central group, including



sensors near the MEF (dark-gray region in Fig. 9A), and peripheral group, including theother sensors (light-gray region). We created data by shuffling the order of trials in
MEG signals of test datasets in the central group (central-shuffling condition) or in the
peripheral group (peripheral-shuffling condition). Because we did not shuffle the orderof trials in the behavioral data, this shuffling removes any correspondence between

MEG signals and fingertip positions. We estimated source currents and reconstructed
the fingertip positions from the currents of test datasets. The inverse-filter andregression weights were the same as those in the previous analyses, i.e., they were
estimated by using the intact training datasets. The above procedure was done for the 50
test datasets. We compared positional error averaged across trials, datasets and subjects



between the conditions. To investigate “baseline” positional error when there was nocorrespondence between the MEG signal in any sensor and the finger position, we also
created data by shuffling the order of trials in signals of test datasets in all sensors(all-shuffling condition) and then conducted the same analysis.

-23-

Our hypothesis is as follows. If signals from the MEF contribute to thereconstruction, the accuracy would degrade in the central-shuffling condition. However,
it would not significantly degrade in the peripheral-shuffling condition because thereconstruction does not depend on information obtained by the peripheral sensors. By

contrast, if signals from peripheral sensors rather than the central sensors contribute to
the reconstruction, we would observe the opposite results. The degree of degradation
caused by the partial (central or peripheral) shuffling can be investigated by comparing
the observed error to the baseline error in the all-shuffling condition.

Results

MEG sensor signals and fMRI activity
Figure 2A shows the vertical and horizontal position of the fingertip when a subject
moved the fingertip toward the top target and back to the initial position. Figure 2B
shows corresponding MEG sensor signals as a function of time. They were aligned tothe movement initiation (0 ms) and averaged across trials and sessions. We confirmed ahighly dipolar pattern of MEF over the contralateral hemisphere to the moving hand at
50 ms after the initiation (a red rectangle in Fig. 2C) in topographic field patterns
sampled at intervals of 50 ms.
We subtracted fMRI activity when subjects observed the screen without finger
movements (observation periods) from activity when subjects moved their fingers

(execution periods). Figure 2D shows the subtracted fMRI activity in the above subject
at the same threshold as that applied to prior information for MEG source localization

-24-(P < 0.001 uncorrected). The activations were mainly observed in sensory-motor
regions, and the highest three activation peaks were observed in the left precentral
lobule. Their peak coordinates (x, y, z) in the MNI reference brain (see Materials andmethods) were (-20, -12, 78), (-36, -24, 72) and (-40, -10, 68), respectively.



Cortical current estimated by the hierarchical Bayesian method
For each subject, by using the hierarchical Bayesian method, we estimated about 2,500
single-current dipoles that were equidistantly distributed on and perpendicular to thecortical surface. Figure 3A shows estimated source currents of the above subject in unitsof current density on the cortical surface (pAm/mm2). Time courses of current density ateach current source (dipole) were aligned to the movement initiation and averaged
across trials. Absolute values of the time course were then averaged across sources ineach ROI in the left hemisphere (Fig. 3B).
Current density in the SMA, PM and M1 gradually increased toward the initiation,
which is consistent with knowledge that those regions are related to planning orexecution of movements. The marked increase of current density in the S1 after the
initiation is consistent with knowledge that the region is related to somatosensory
information processing. Peak latencies observed in these time courses are consistent
with findings in the previous studies. The latencies in PM and M1 around 50 ms before
the movement initiation were consistent with those observed in the previous study of

MEG source current (Cheyne et al., 2006; Huang et al., 2004). Peak latencies in M1 andS1 around 100 - 200 ms after the initiation were also found in studies using MEG

-25-source currents (Cheyne et al., 2006) and ECoG (Kunieda et al., 2000). Amplitudes ofthe estimated current density were comparable to those in previous studies investigating
MEG source currents in sensorimotor (Huang et al., 2004) and visual (Yoshioka et al.,2008) regions.



Performance of reconstruction models of fingertip position from cortical source
currents
We estimated reconstruction models based on a multivariate liner regression while
controlling the over-fitting problem by using different methods, and then we compared
the positional errors in test datasets between the methods (see Materials and methods).



Figure 4 shows positional errors averaged across trials, datasets and subjects for the
regularized least-squares regression method and the sparse regression method, with orwithout anatomical constraint. These results indicate that performance improves as thenumber of training datasets increases, regardless of the regression type, since the error
level decreased as the ratio of the number of trials in the training dataset to that in thetotal data increased from 5% (Fig. 4A) and 20% (Fig. 4B) to 70% (Fig. 4C). Regarding
the regularized least-squares regression, its performance improved by increasing thevalue of a regularization parameter ( ! : relative importance of penalty for large weightvalues, see Eq. 3) from 1 to 104 when the size of the training dataset is small (Fig. 4A).
Anatomical selection of source currents (an anatomical constraint, see Materials andmethods) is also effective (solid line) for reconstruction accuracy when the size of the
training dataset is small. That is, for 5% training data, the anatomical constraint reduced
the error not only in the regularized least-squares regression but also slightly in thesparse regression (22.68 ± 1.21 mm with and 22.70 ± 1.50 mm without the constraint,
black arrow). These results indicate the importance of a priori anatomical knowledge in-26-

selecting lattice points when the size of the training dataset is small. Independently of
the size of the training dataset, the sparse regression was always more accurate than theregularized least-squares regression (Figs. 4A - 4C). The best performance was obtainedby the sparse regression without the anatomical constraint when the ratio of trainingdatasets was 70% (white arrow). We then closely examined the performance of thereconstruction model in this most accurate case, as described below.

Reconstruction model estimated by the sparse regression method
Figure 5 shows examples of the fingertip trajectory reconstructed by the model fromtest data (red line) in comparison with the actual trajectory recorded by the motion
tracking system (blue line) when a subject moved the fingertip toward the eight
directions indicated by black arrows (see also Supplementary Video). Figure 6A shows
positional error (Eq. 4) averaged across trials, datasets and subjects. The error was 12.7mm (SD: 1.66) for training dataset and 14.6 mm (SD: 2.55) (same as that indicated bywhite arrows in Fig. 4C) for test dataset when the movement’s path length was longerthan 200 mm.
For comparison, we also made a model that reconstructs the position directly from
MEG sensor signals (see Materials and methods). Open bars in Figure 6A indicate
errors when the position was reconstructed from sensor signals. These errors weresignificantly higher than those when positions were reconstructed from cortical currents(filled bars) for both training (cyan: F(1,4) = 72.1, P < 0.002) and test (red: F(1,4) =45.0, P < 0.005) datasets according to a two-way (datasets x type of signal used forreconstruction: cortical source current or MEG sensor signal) repeated-measures
ANOVA.

-27-

Figure 6B shows coefficients of determination averaged across trials, datasets andsubjects for the horizontal (upper panel) and vertical (bottom panel) directions. The
coefficients for the horizontal direction were significantly higher when the position wasreconstructed from source currents (filled bars) than when it was reconstructed from

sensor signals (open bars) for both training (cyan: F(1,4) = 150.70, P < 0.0005) and test
datasets (red: F(1,4) = 73.1, P < 0.001). A significant difference was also identified in
the coefficients for the vertical direction (F(1,4) = 121.14, P < 0.0005 for training
datasets; F(1,4) = 74.5, P < 0.001 for test datasets) between the source currents and thesensor signals. The coefficients were significantly higher in the vertical direction thanthe horizontal direction for both training datasets (F(1,4) = 9.72, P < 0.05) and test
datasets (F(1,4) = 9.34, P < 0.05) when the position was reconstructed from source
currents. A higher coefficient for the vertical direction was also observed when theposition was reconstructed from sensor signals (F(1,4) = 14.5, P < 0.03 for training
datasets; F(1,4) = 15.1, P < 0.03 for test datasets).



Anatomical localization of source currents contributing to reconstruction
We investigated weight values for current sources in the reconstruction model for thevertical component (Eq. 2), in which the coefficient of determination was higher thanthe other component (Fig. 6B). Figure 7A shows values averaged across time points anddatasets and projected on an inflated model of the cortical surface. Highly positive andnegative weights were concentrated in the pre- and postcentral regions and the parietal
regions. Among the parietal regions, high weight values were found in the superiorparietal lobule (SPL). Figure 7C shows absolute weight values summed within eachanatomical region and averaged across subjects. The highest values were found in theparietal and central regions of the contralateral hemisphere. Figure 7B shows weight

-28-

values for MEG sensors when the position was reconstructed from sensor signals. Highand low values were scattered around all sensors, and it was difficult to determine
which region is important for the reconstruction.

Temporal change of weight values in reconstruction modelFigure 8A shows weight values for individual current sources in representative ROIs as
a function of advanced time from the time point when a fingertip position isreconstructed (0 ms). Thin black lines indicate time courses of weight values averaged
across datasets and sources within an ROI. Each line corresponds to a subject. We fitted
an exponential curve (Eq. 7 and bottom right panel in Fig. 8A) to the time course beforeaveraging across datasets and estimated parameters (a, b and ! ) using a least-squares
method (see Materials and methods). Thick red curves show exponential functions
yielded by parameter values averaged across datasets and subjects (note that the red
curves were not directly fitted to the thin black lines). Weight values gradually
increased at an early stage in regions anterior to the central sulcus (SMA, PM and M1;upper panels in Fig. 8A), while values abruptly increased at a late stage near 0 ms inregions posterior to the sulcus (S1 and parietal regions; lower panels). Figure 8B showstime constants ( ! -values) averaged across datasets and subjects, indicating that the! -values in the anterior regions (black bars) were larger than those in the posterior
regions (white bars). This quantitatively confirms the above observation that weightvalues increased at an early stage in the anterior regions.

Examination of contribution of MEF to reconstruction
We examined whether the signals contributing to the reconstruction originated from the

-29-

MEF by shuffling the order of trials in MEG signals from sensors in the central region,
including the MEF, or in the peripheral region. For comparison, we also shuffled theorder of trials in signals of all sensors (see Materials and methods, and Fig. 9A). We
reconstructed the positions from the shuffled and intact test data and calculated

positional errors averaged across trials, datasets and subjects (Fig. 9B). We applied
two-way repeated-measures ANOVA to the error in individual trials (sets x conditions:all-, central- and peripheral-shuffling conditions and intact condition). Consequently,
the effect of the condition was significant (F(3, 12) = 41.9, P < 0.0001). Tukey’s HSD
post hoc test identified a significant difference in any comparison between theconditions, except the comparison between the central- and all-shuffling conditions (P <0.05 level). The error in the central-shuffling condition (black bar) was significantly
higher than that the peripheral-shuffling condition (gray bar), suggesting that morecrucial information for the reconstruction came from the central region including the
MEF than from the peripheral region. Moreover, the error in the central-shuffling

condition was not significantly different from that in the all-shuffling condition (leftwhite bar). This indicates the importance of information originating from sensors in thecentral region: If the information were unavailable (central-shuffling condition), errorsof reconstruction models would not be different from the “baseline” error when no
useful information can be obtained from any sensor (all-shuffling condition).

Discussion
We reconstructed the trajectory of the fingertip from the time courses of source currents
estimated on the cortical surfaces. We found that a reconstruction model estimated by asparse regression has better generalization ability than models estimated by a

-30-

regularized least-squares regression, and thus we closely examined the performance ofthe former. Consequently, the error between the predicted and actual positions for thetest dataset was 14.6 mm averaged across subjects, which was 7% of the path length of
the required movement (200 mm), and the averaged coefficient of determination was

0.68. These performances were significantly superior to that of a comparable modelusing MEG sensor signals. We believe the better performance was due to the selective
use of information directly related to sensory processing and motor control forreconstruction based on the following reasons. First, important (highly positive andnegative) weights were mainly found in the pre- and post-central regions and theparietal regions in the contra-lateral hemisphere to the moving hand when the positions
were reconstructed from cortical currents (Figs. 7A and 7C). Significant weight values
were found in not only the primary sensorimotor regions but also in the parietal region(SPL). Many studies have indicated a contribution of the SPL to the planning of action(Culham et al., 2006; Field et al., 2007; Medendorp et al., 2005), suggesting that

reconstruction was partly based on activity related to the planning of movement. By
contrast, when they were predicted from MEG sensor signals, high and low values were
scattered over many sensors (Fig. 7B), and thus it was difficult to infer which sensors
were important for the reconstruction. Second, we shuffled the order of trials in MEG
signals from sensors in the central region, including the MEF, and found a significantdecrease in reconstruction performance (Fig. 9), suggesting that reconstruction depends
on signals originating from the MEF.
Georgopoulos and colleagues (Georgopoulos et al., 2005) reconstructed a cursor
position from a linear weighted summation of MEG sensor signals when subjects
continuously moved a joystick while drawing a pentagon at a constant velocity. By

-31-

contrast, our study reconstructed a fingertip position from cortical currents. While
fingertip velocity ranged from 0 up to 700 mm/s, we succeeded in reconstructing its
position at various velocities using a single set of weights. Bradberry and colleagues(Bradberry et al., 2009) reconstructed hand velocity from MEG sensor signals duringcenter-out reaching tasks, and the mean correlation coefficients (CCs) between
reconstructed and actual velocity were 0.48 and 0.32 for horizontal and verticalcomponents, respectively. For comparison, CCs calculated from our results ofreconstructing the position using source currents were 0.74 and 0.81, and those fromreconstructing the velocity were 0.49 and 0.56.

The advantage of the hierarchical Bayesian method is to estimate a large numberof current sources densely distributed on the cortical surface. The effectiveness of thismethod has already been confirmed by the successful reconstruction of retinotopic
activities in primary visual areas (Yoshioka et al., 2008). We used this method for the
first time to estimate cortical currents evoked by a sensorimotor task and confirmed thevalidity of the estimation by comparing time courses of the estimated currents insensorimotor regions to time courses of activity reported by previous studies that usedother methods for MEG source localization (Cheyne et al., 2006; Huang et al., 2004)
and ECoG (Kunieda et al., 2000).
Our proposed method for reconstruction of movements from cortical currents

after solving the inverse problem can contribute to BMI using non-invasive techniqueswith high temporal resolution as follows. First, the accuracy of reconstruction increases
in comparison to reconstruction from sensor signals as the current study demonstrated.
A study of cursor control using a combination of amplitudes of EEG rhythms alsosuggested that the reliability of control increases when the amplitude is calculated from-32-

cortical currents in comparison to EEG sensor signals (Cincotti et al., 2008). Second,
simultaneous estimation of cortical and artifact currents can remove artifacts originatingfrom outside of the brain, such as eye and cardiac movements and muscle activity(Fujiwara et al., 2009; Morishige et al., 2009). Third, knowledge in neuroscience aboutbrain regions and their functions could be applied to the selection of feature signals. Our
results indicate that anatomical selection of lattice points for cortical currents based on
knowledge in neuroscience is an effective way to improve the ability of reconstruction
when the size of the training dataset is small (Fig. 4A). The regularized least-squares
regression with anatomical constraint could maintain a relatively constant performance(solid line in Fig. 4A) regardless of the value of the regularization parameter, whereas
the reconstruction error of the regression without the constraint increased as the value
became small (broken line). This result systematically shows that the effect ofanatomical selection is more prominent if the selected machine-learning algorithm isless efficient in the selection of feature signals, which corresponds to smaller values ofthe regularization parameter. When considering practical clinical applications of BMIs,
it might be difficult to obtain many training data from patients or to use a sophisticatedlearning algorithm such as a sparse regression for selecting feature signals due to theinsufficient time and power for computation. Results of our analysis suggest thatanatomical selection would be effective in such cases.

Several possibilities remain regarding the exact neural information that was
decoded to reconstruct movement trajectories in our experiment. Because we fixed
finger joints (Fig. 1C), the fingertip positions projected to the x-y plane (the extrinsic
coordinates) uniquely correspond to the wrist joint angles (the intrinsic coordinates) in a
one-to-one manner. Thus, these two representations cannot be behaviorally dissociated.

-33-

We speculate that information in both representations contributed to the reconstruction
because the weights in sparse linear regression were widely distributed across theprimary motor and parietal regions (Figs. 7A and 7C). Reconstruction accuracy was
higher in the vertical direction than in the horizontal direction (Fig. 6B), which isconsistent with a previous study reconstructing three-dimensional arm movements ofmonkeys from ECoG activity (Chao et al., 2010). A possible reason is that the activity
level of muscles markedly changes during vertical movement due to gravity. If so,
neuronal activity directly related to the control of muscle tensions also contributes to thereconstruction. Our analysis of behavioral data revealed a significant effect ofmovement direction on movement time (see Materials and methods). However,
movement times for only one or two targets are significantly different from that for eachtarget. Thus, information related to movement time has coarse spatial information,
which is unlikely to be crucial for our reconstruction of movements with positional
error of 14.6 mm.

Although our current study reconstructs movement in offline processing,
online/real-time reconstruction is necessary for application to BMIs and brain-computerinterfaces (BCIs). Our method has the following potential advantages for online
processing. First, we do not need to solve the inverse problem for every trial but canapply an inverse filter already estimated from previous datasets to sensor signals in eachtrial. Since the filter is a linear transformation, the time needed for calculation of thecortical currents is almost negligible. Second, reconstruction models of finger positions
from cortical currents are also simple linear transformations (Eqs. 1 and 2). However,calculation of the inverse filter and reconstruction models from the training datarequires almost a full day. Consequently, for online reconstruction, the filter and models

-34-

have to be generalized to test data on a different day. Because it is impossible to place
the subjects’ heads in the MEG scanner exactly in the same position every time, we will
have to develop a method to adjust the inverse filter within a short time period (during
which subjects can wait in the scanner). Fortunately, our hierarchical Bayesian methodestimates current variance on the cortical surface, which does not depend on the MEG
sensor position. The inverse filter corresponding to new sensor position can becalculated using previously estimated current variance and a new lead-field matrix
corresponding to the new sensor position (see Supplementary Information 1). This
calculation can be done in a short time. Furthermore, theoretically, we can use the samereconstruction models of movements even if the head position changes, since themodels are based on cortical currents instead of sensor signals. This is the third
advantage of our method. However, we still have to confirm this assumption.
Furthermore, offline preprocessing of the data such as noise filtering, normalization of
signals, and rejection of impaired sensors and trials contaminated by eye movementsand blinks (see Materials and methods) will have to be adjusted for online processing.
The generalization ability of reconstruction models to subsequently obtained test
datasets is also important for online applications to BMI/BCI. Although we randomly
assigned trials to training dataset (70%) and test dataset (30%) regardless of thechronological order in the main analysis, in an additional analysis we examined thegeneralization ability when the early 70% and the later 30% of trials were assigned tothe training and test datasets, respectively. Consequently, positional error averaged
across subjects was 14.8 mm (SD: 2.52), which was comparable to the error in the main
analysis (14.6 mm, SD: 2.55). This result further suggests the effectiveness of ourreconstruction model for online applications.

-35-

When considering application of our method to BMI/BCI, it would be useful if we
could avoid the need to acquire fMRI data. In another additional analysis, we tried toreconstruct movements using the same method as that in the main analysis but without
incorporating fMRI activity as a hierarchical prior in the estimation of the inverse filter.
The resulting positional error averaged across subjects was 13.7 mm (SD: 2.56), whichis smaller than that in the main analysis (14.6 mm, SD: 2.55). We confirmed that the
distribution and time courses of estimated cortical currents were not different fromthose in the main analysis (Fig. 3A). This result appears preferable for simplification of
our method and, to some degree, expected by a previous study indicating that ahierarchical Bayesian method can solve the inverse problem without fMRI activity
because it uses the activity as a soft constraint on the variance of source currents (Satoet al., 2004). Possible reasons for the small contribution of fMRI data to the estimationin the current study are as follows. In our pointing task, the pattern of cortical activity is
expected to dynamically change according to various phases of movement (e.g.,planning, execution and feedback) within several seconds. Our fMRI scanning with arepetition time of 3.5 s (see Materials and methods) may have been insufficient formeasurement of the rapid change in activity pattern. On the other hand, regions
activated by the task were scattered around a wide cortical region (Figs. 3A, 7A and 7C),
and thus cortical source currents could be estimated mainly from MEG data havingrelatively low spatial resolution. However, it is premature to conclude that fMRI data is
unnecessary for estimation of cortical currents by a hierarchical Bayesian method, sincea previous study using the method showed that fMRI data plays a crucial role in finespatiotemporal reconstruction of retinotopic activities in the primary visual area
(Yoshioka et al., 2008). Our preliminary data using EEG instead of MEG, which isessential for portable BMI/BCI applications, suggests that fMRI data is needed to-36-

estimate cortical currents related to hand movements.
In the current study, we aimed to reconstruct continuous trajectories during
pointing movements by using linear regression models, but we did not examine whether
our method could estimate the intended target before movement initiation. By contrast,

some non-invasive BMI/BCI methods for communication prostheses decode just theintended target before reaching it by classification of MEG and EEG activity patterns(e.g., Hammon et al., 2007; Jerbi et al., 2007). Human movements do not consist of
only discrete ones, such as rapid aiming at a static target, but also continuous ones suchas drawing a figure or tracking a moving target. Although our current study

reconstructed nearly linear trajectories toward a target, our method can be extended toreconstruction of continuous trajectories of various shapes. Prediction of the goal
(end-point) of movement by classification of activity patterns and reconstruction of theintended trajectory by regression of activity time courses could be complementaryelements for effective motor prostheses.



If we could reconstruct movements from source currents in smaller regions thanthose in the current study, it might be possible to extend our approach to examining
what information is represented in a particular region by investigating which type of
movement parameters (e.g., position, velocity and acceleration of the fingertip, andtarget direction) can be extracted from source current signals in that region. Such anextension would be a powerful tool in computational neuroscience.

Conclusion
Few previous studies have adopted reconstruction of movements from cortical currents

-37-

estimated from EEG or MEG signals, probably because calculation of inverse filtersfrom sensor signals to cortical currents requires considerable time and complicates dataprocessing. If estimations of the inverse filter and the reconstruction models consistedof linear components only, there would be little merit in calculating cortical currents.
However, our study demonstrated that non-linear selection of feature signals in corticalspaces based on sparse regression and knowledge in neuroscience can help to improvereconstruction accuracy. Such selection becomes possible after thousands of source
currents are estimated on the cortical surface. Furthermore, the best reconstruction
performance in the current study was obtained when we applied a sparse regression

(rather than a regularized least-squares regression) to cortical currents (rather than MEG
sensor signals). Since both a sparse regression and an inverse filter for calculating
cortical current are based on sparse Bayesian estimation, our results suggest that a
combination of the two sparse Bayesian methods can provide a key technology for
non-invasive BMI without depending too much on user training.

Acknowledgements
This work was partly supported by Honda Research Institute Japan and Fundamental
Technology Research Center of Honda Research and Development Corporation. We
thank Yuka Furukawa and members of the Brain Activity Imaging Center of Advanced
Telecommunications Research Institute-Promotions for technical support and Dr.

Toyama for helpful comments.

-38-

References

Akaike, H., 1974. A new look at the statistical model identification. IEEE Transactions
on Automatic Control AC-19, 716-723.
Bradberry, T.J., Gentili, R.J., Contreras-Vidal, J.L., 2010. Reconstructing
three-dimensional hand movements from noninvasive electroencephalographic
signals. J Neurosci 30, 3432-3437.
Bradberry, T.J., Rong, F., Contreras-Vidal, J.L., 2009. Decoding center-out handvelocity from MEG signals during visuomotor adaptation. Neuroimage 47,

1691-1700.
Chao, Z.C., Nagasaka, Y., Fujii, N., 2010. Long-term asynchronous decoding of arm
motion using electrocorticographic signals in monkeys. Front Neuroengineering
3, 3.
Cheyne, D., Bakhtazad, L., Gaetz, W., 2006. Spatiotemporal mapping of corticalactivity accompanying voluntary movements using an event-related
beamforming approach. Hum Brain Mapp 27, 213-229.
Cincotti, F., Mattia, D., Aloise, F., Bufalari, S., Astolfi, L., De Vico Fallani, F., Tocci,
A., Bianchi, L., Marciani, M.G., Gao, S., Millan, J., Babiloni, F., 2008.
High-resolution EEG techniques for brain-computer interface applications. JNeurosci Methods 167, 31-42.
Culham, J.C., Cavina-Pratesi, C., Singhal, A., 2006. The role of parietal cortex invisuomotor control: what have we learned from neuroimaging?
Neuropsychologia 44, 2668-2684.
Field, D.T., Wilkie, R.M., Wann, J.P., 2007. Neural systems in the visual control ofsteering. J Neurosci 27, 8002-8010.
Fujiwara, Y., Yamashita, O., Kawawaki, D., Doya, K., Kawato, M., Toyama, K., Sato,
M.A., 2009. A hierarchical Bayesian method to resolve an inverse problem of
MEG contaminated with eye movement artifacts. NeuroImage 45, 393-409.
Geman, S., Bienenstock, E., Doursat, R., 1992. Neural networks and the bias/variance-39-

dilemma. Neural Comput 4, 1-58.
Georgopoulos, A.P., Langheim, F.J., Leuthold, A.C., Merkle, A.N., 2005.
Magnetoencephalographic signals predict movement trajectory in space. Exp
Brain Res 167, 132-135.


Hammon, P.S., Making, S., Poizner, H., Todorov, E., de Sa, V.R., 2007. Predicting
reaching targets from human EEG. IEEE Signal Processing Magazine 25, 69-77.
Hochberg, L.R., Serruya, M.D., Friehs, G.M., Mukand, J.A., Saleh, M., Caplan, A.H.,
Branner, A., Chen, D., Penn, R.D., Donoghue, J.P., 2006. Neuronal ensemble
control of prosthetic devices by a human with tetraplegia. Nature 442, 164-171.



Huang, M.X., Harrington, D.L., Paulson, K.M., Weisend, M.P., Lee, R.R., 2004.
Temporal dynamics of ipsilateral and contralateral motor activity during
voluntary finger movement. Hum Brain Mapp 23, 26-39.
Jerbi, K., Lachaux, J.P., N'Diaye, K., Pantazis, D., Leahy, R.M., Garnero, L., Baillet, S.,
2007. Coherent neural representation of hand speed in humans revealed by MEGimaging. Proc Natl Acad Sci U S A 104, 7676-7681.
Kamitani, Y., Tong, F., 2005. Decoding the visual and subjective contents of the humanbrain. Nat Neurosci 8, 679-685.
Krepki, R., Blankertz, B., Curio, G., Müller, K., 2006. The Berlin Brain-Computer
Interface (BBCI) towards a new communication channel for on-line control ingaming applications. International Conference on Distributed Multimedia
System (DMS’03).
Kunieda, T., Ikeda, A., Ohara, S., Yazawa, S., Nagamine, T., Taki, W., Hashimoto, N.,
Shibasaki, H., 2000. Different activation of presupplementary motor area,
supplementary motor area proper, and primary sensorimotor area, depending on

the movement repetition rate in humans. Exp Brain Res 135, 163-172.
Maldjian, J.A., Laurienti, P.J., Kraft, R.A., Burdette, J.H., 2003. An automated method
for neuroanatomic and cytoarchitectonic atlas-based interrogation of fMRI data
sets. NeuroImage 19, 1233-1239.
Medendorp, W.P., Goltz, H.C., Crawford, J.D., Vilis, T., 2005. Integration of target andeffector information in human posterior parietal cortex for the planning of action.

-40-

J Neurophysiol 93, 954-962.
Miyawaki, Y., Uchida, H., Yamashita, O., Sato, M.A., Morito, Y., Tanabe, H.C., Sadato,N., Kamitani, Y., 2008. Visual image reconstruction from human brain activity
using a combination of multiscale local image decoders. Neuron 60, 915-929.
Morishige, K., Kawawaki, D., Yoshioka, T., Sato, M.A., Kawato, M., 2009. Artifact
removal using simultaneous current estimation of noise and cortical sources.
Lecture Notes in Computer Science: ICONIP 2008 Springer-Verlag, Berlin, pp.
335-342.
Musallam, S., Corneil, B.D., Greger, B., Scherberger, H., Andersen, R.A., 2004.



Cognitive control signals for neural prosthetics. Science 305, 258-262.
Nambu, I., Osu, R., Sato, M.A., Ando, S., Kawato, M., Naito, E., 2009. Single-trial
reconstruction of finger-pinch forces from human motor-cortical activation
measured by near-infrared spectroscopy (NIRS). Neuroimage 47, 628-637.
Nicolelis, M.A., 2001. Actions from thoughts. Nature 409, 403-407.

Sato, M.A., 2001. Online Model Selection Based on the Variational Bayes. NeuralComputation 13, 1649–1681.
Sato, M.A., Yoshioka, T., Kajihara, S., Toyama, K., Goda, N., Doya, K., Kawato, M.,
2004. Hierarchical Bayesian estimation for MEG inverse problem. NeuroImage23, 806-826.

Taylor, D.M., Tillery, S.I., Schwartz, A.B., 2002. Direct cortical control of 3D
neuroprosthetic devices. Science 296, 1829-1832.
Ting, J.A., AD’Souza, Yamamoto, K., Yoshioka, T., Hoffman, D.S., Kakei, S., Sergio,
L.E., Kalaska, J.F., Kawato, M., Strick, P.L., Schaal, S., 2005. Predicting EMG
Data from M1 Neurons with Variational Bayesian Least Squares. Advances inNeural Information Processing Systems.
Ting, J.A., D'Souza, A., Yamamoto, K., Yoshioka, T., Hoffman, D., Kakei, S., Sergio,L., Kalaska, J., Kawato, M., Strick, P., Schaal, S., 2008. Variational Bayesian
least squares: An application to brain-machine interface data. Neural Netw 21,1112-1131.

-41-

Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix,N., Mazoyer, B., Joliot, M., 2002. Automated anatomical labeling of activations
in SPM using a macroscopic anatomical parcellation of the MNI MRI
single-subject brain. NeuroImage 15, 273-289.
Waldert, S., Preissl, H., Demandt, E., Braun, C., Birbaumer, N., Aertsen, A., Mehring,C., 2008. Hand movement direction decoded from MEG and EEG. J Neurosci 28,
1000-1008.
Wolpaw, J.R., McFarland, D.J., 2004. Control of a two-dimensional movement signal
by a noninvasive brain-computer interface in humans. Proc Natl Acad Sci U S A101, 17849-17854.
Yoshioka, T., Toyama, K., Kawato, M., Yamashita, O., Nishina, S., Yamagishi, N.,
Sato, M.A., 2008. Evaluation of hierarchical Bayesian method throughretinotopic brain activities reconstruction from fMRI and MEG signals.
NeuroImage 42, 1397-1413.



-42-

Figure legends

Figure 1(A) Eight directions (gray arrows) in which subjects were required to move the index
fingertips. Terminals of the arrows correspond to target positions at which subject
aimed. (B) Posture of subjects in MEG scanner (scanner not drawn in figure). The
experimental setup in the MR scanner was the same except that subjects viewed avertical screen placed near their heads though a mirror. (C) Fixation of joints(metacarpophalangeal joint, distal and proximal interphalangeal joints) of the right

index finger. Surgical tape was used to fix a wooden brace on the back of the finger andhand so that only wrist-joint movements were allowed. (D) Timeline in a trial of MEG
experiment. Top panels show examples of images projected on the screen. Bottom
panels illustrate states of subjects' hands. (E) A part of the timeline in a session of the
fMRI experiment. Each session had 16 alternating blocks of execution and observation



periods.
Figure 2(A) Time courses of the horizontal (blue) and vertical (red) positions of the fingertip
recorded by a motion tracking system when a subject moved the fingertip toward the
top target. The time courses were aligned to the movement initiation (0 ms) and

averaged across trials and sessions. (B) Time courses of MEG-sensor signals aligned tothe movement initiation and averaged across trials and sessions. (C) Topographic field
patterns (corresponding to Fig. 2B) sampled at intervals of 50 ms ranging from 50 ms
before to 200 ms after the movement initiation. A red rectangle indicates the clearest

-43-

pattern of the movement-evoked field (MEF). (D) fMRI activity when the same subject
moved the fingertip toward eight targets in execution periods. We contrasted activity inthe execution periods to that in the observation periods when the subjects observed
replay of screen images recorded in the preceding execution period (P < 0.001
uncorrected). The activity map was rendered on inflated cortical surfaces(lateral-posterior view). Sulci (concave) and gyri (convex) are indicated by dark andlight gray shading, respectively.
Figure 3(A) Cortical source currents estimated by a hierarchical Bayesian method with the samesubject as Figure 2. Cortical current signal was transformed into current density. The
top panels show a current density map rendered on inflated cortical surfaces(lateral-posterior view). The current density was averaged within periods of 150–50 ms
before, 50–0 ms before, and 50–100 ms after the movement initiation separately for
individual current sources on the cortical surface. The bottom panel shows time coursesof cortical currents aligned to the movement initiation and averaged across sources in aregion of interest (ROI). (B) ROIs rendered on the left cortical surface. S1, primarysomatosensory cortex; M1, primary motor cortex; SMA, supplementary motor area; PM,premotor regions. (C) A linear reconstruction model of fingertip positions (horizontal orvertical) from cortical currents. Positions were reconstructed from the weighted sum of

cortical currents sampled at 11 time points (open circles) at interval of 20 ms within 200
ms preceding the currently reconstructed position.

Figure 4

-44-

Positional error between predicted and actual position averaged across trials, datasets
and subjects for test datasets (±SD). Prediction model was estimated by using different
sizes of training datasets, i.e., ratio of number of trials in training dataset to that in total
data was 5% (A), 20% (B) or 70% (C). Abscissa indicates types of regression methodused for estimation, and numbers (1 - 104) indicate the value of a regularization
parameter in regularized least-squares regression. Black arrow indicates errors ofmodels estimated by a sparse regression when the ratio was 5%. White arrow indicatesthe best performance (smallest error) in the figure.

Figure 5

Example of fingertip trajectories in the horizontal (left panels) and vertical (rightpanels) directions when a subject moved the fingertip toward eight directions as
indicated by black arrows. Red lines indicate trajectories predicted by a linear modelthat reconstructs the fingertip position from cortical source current. Blue lines indicate
an actual trajectory recorded by a motion tracking system.

Figure 6(A) Positional error between predicted and actual positions averaged across trials,
datasets and subjects for test and training data (+SD). Filled bars indicate error made bya prediction model using cortical source currents. Open bars indicate error made by amodel using MEG sensor signals. (B) Coefficient of determination for the models

averaged across trials, pairs of datasets and subjects for test and training data. Upperand lower panels show coefficients for the horizontal and vertical directions,
respectively. **: P < 0.005, ***: P < 0.001 according to a repeated-measures ANOVA.

-45-

Figure 7(A) Weight values for cortical source currents in a model predicting the fingertip
positions from the source current in the same subject as Figs. 2, 3 and 5. These values
were averaged across time points and datasets and rendered on the cortical surface. (B)

A topographic pattern of weight values for MEG sensors in a model predicting theposition from the sensor signals in the same subject. The averaging is the same as above.
(C) Absolute weight values for cortical source currents summed within an anatomical
ROI and averaged across subjects (+SD; see Materials and methods for definition of
anatomical ROIs).



Figure 8(A) Time courses of weight values for individual current sources as a function ofadvanced time from the time point when a fingertip position is predicted (0 ms). Eachthin black line corresponds to a subject and indicates a time course of weight value
averaged across datasets and sources within an ROI for the subject. Thick red curvescorrespond to exponential functions whose parameter values were estimated from timecourses of weight values before averaging across datasets and then averaged acrossdatasets and subjects. The bottom right panel illustrates how an exponential curve (w =a + b • exp(t / ! )) changes depending on each parameter. (B) ! -value averaged acrossdatasets and subjects for each ROI (+SD). Black or white bars indicate that ROIs areanterior or posterior to the central sulcus, respectively.

Figure 9(A) Definition of a central region (dark gray), including sensors near the-46-

movement-evoked magnetic field (MEF), and a peripheral region (light gray), including
the remaining sensors. The MEF is represented by a color-coded contour map. We
created this map by averaging MEG signals separately for each sensor at 50 ms after the
movement initiation across trials and subjects. (B) Error between predicted and actual

positions when the order of trials in MEG signals from sensors in the central region(black bar) or in the peripheral region (gray bar) was shuffled in relation to when theorder was not shuffled in either region (intact, right white bar) or when the order oftrials in signals from all of the sensors were shuffled (left white bar). *: P < 0.05, N.S.:
not significant according to Tukey's HSD post hoc test.



-47-



Figure 1

-48-






Figure 2

-49-
Figure 3


-50-











Figure 4





-51-


Figure 5

-52-




Figure 6

-53-







Figure 7

-54-







Figure 8

-55-







Figure 9

-56-



Supplementary Information

Re: “Reconstruction of two-dimensional movement trajectories from selectedmagnetoencephalography cortical currents by combined sparse Bayesian
methods” by Akihiko Toda, Hiroshi Imamizu, Mitsuo Kawato and Masa-aki Sato

1. Hierarchical Bayesian estimation


2. Sparse linear regression
3. Supplementary video showing examples of reconstructed trajectory
4. References

-57-

1. Hierarchical Bayesian estimation

The hierarchical Bayesian method (Sato et al., 2004; Yoshioka et al., 2008) estimatescortical activity from MEG sensor signals in a distributed source model, in which thecortical current is modeled by a number of current dipoles with fixed position andorientation. Generally, the position and direction of the dipoles are determined based on
the widely accepted assumptions that the principal sources of MEG signals are
post-synaptic currents arising from pyramidal neurons in the cortex and that theirdirection is perpendicular to the cortical surface (Dale and Sereno, 1993).
The MEG inverse problem is generally ill-posed, since the number of current

dipoles (several thousand) is often much larger than the number of MEG sensors(several hundred). To overcome this problem, fMRI information has been widely usedas prior information (Dale et al., 2000; Kajihara et al., 2004; Phillips et al., 2002). Fromthe viewpoint of Bayesian estimation, fMRI information determines the priordistribution of the variance in the current amplitude. In the hierarchical Bayesianestimation, the variance of the current amplitude is also parameterized and estimatedfrom the MEG signal. fMRI information is imposed on the prior distribution of thevariance in the current amplitude through hyperparameters rather than the variance itself
so that it can give a soft constraint on the variance.
MEG forward model

Under the quasi-static approximation assumption (Hamalainen et al., 1993), the MEG
forward model, that is, the relationship between the amplitude of current dipoles andobserved magnetic field at time point t, is given by

-58-

B(t) = G ! J(t) ,(1)where B(t) is an N-by-1 vector for the magnetic field measured at N sensors, J(t) is
an I-by-1 vector for the current amplitudes at I current dipoles, and G is an N-by-I
matrix, referred to as the lead field matrix. In this study, the inner skull surface wasapproximated as a sphere, and Sarvas's equation was used to calculate lead fields(Sarvas, 1987).
Estimating current varianceIn the hierarchical Bayesian method, the current amplitude is estimated by introducing
an Automatic Relevance Determination (ARD) hierarchical prior (Neal, 1996):' "*P ( J(t) ! , " ) # exp ) $ J%(t) & A & J(t) ,
( 2
+



P (! i ) = " (! i ! 0i , # 0 )P (! ) =

(2)

1
!,

where ! is the inverse noise variance of the observed MEG signal, A = diag(!), and !
is an I -by-1 vector whose component ! i is the inverse current variance correspondingto the i-th current dipole. " represents the Gamma distribution with mean ! 0i anddegree of freedom ! 0 . Intuitively, the hyper-parameter ! 0 represents confidence of"1
the hierarchical prior information. A prior current variance v0i = ! 0i
represents theprior information on current intensity. For large and small v0i , estimated current J i (t)tends to be large and small, respectively. These values were determined from the fMRI


information:v0i = vbase + ( m0 ! 1) " vbase " (tˆi )2 ,

(3)

-59-

where tˆi is a normalized T-value on the i-th vertex. Normalized T-values are computedby dividing the original T-values by the maximum of those T-values (thus ranging from0 to 1). vbase is a baseline of the current variance, which is estimated from thepre-movement interval (1.0 s ~ 0.5 s before the movement initiation) of the MEG databy a Bayesian minimum norm estimation. A variance magnification parameter m0 ,which is the other hyper-parameter, specifies the scaling between the current variancesin the baseline and task periods. We used m0 = 100 and ! 0 = 10.
Due to the hierarchical prior, the estimation problem becomes nonlinear andcannot be solved analytically. Therefore, we employed the Variational Bayesian (VB)

method (Attias, 1999; Sato, 2001). In the VB method, J(t) , ! and ! are iteratively
updated until convergence.
Spatial smoothness constraint
We assumed a spatial smoothness constraint on the current distribution along with thecortical surface. To do this, we employed a smoothing filter matrix



(

)

Wij ! exp "dij2 / R 2 , where dij is the distance between the i-th and j-th current
dipoles. The smoothing radius parameter, R, was set to give a full-width at
half-maximum (FWHM) of 6 mm. By introducing an auxiliary variable Z ( t ) andlettingJ(t) = W ! Z ( t ) ,



(4)

equation (1) can be replaced by
B ( t ) = Ĝ ! Z ( t ) ,(5)where Ĝ is a smoothed lead-field matrix. Therefore, the MEG inverse problem ischanged to an estimation of Z ( t ) , instead of J(t) , with the smoothed lead-field matrix

-60-

Ĝ . After estimating Z ( t ) , the actual current amplitude J(t) is calculated using
equation (4).
This method smoothes a lead-field matrix and estimates hypothesized unsmoothcortical currents that are smoothed after the estimation. Our incorporation of asmoothness constraint as a soft constraint on the estimation of source currents may
seem simple. Actually, this method was adopted in our study that originally proposed ahierarchical Bayesian method (Sato et al., 2004). However, we found that the method
could reduce the number of parameters and save computational time without degrading
estimation accuracy (Yoshioka et al., 2008).

Specifically, in our earlier paper (Sato et al., 2004), we introduced a soft constraint
for the spatial smoothness constraint:' 1*
P ( J Z, ! ) " exp ) # ( J # W $ Z )% $ & $ ( J # W $ Z ) , 
( 2
+
& 1
)
P ( Z ! ) " exp ( # Z$ % A % Z + ,
*
' 2

(6)

where ! = diag( ! ) and A = diag(#). This prior is equivalent to the following prior as

explained in the paper (Sato et al., 2004).

( 1+
P ( J ! , " ) # exp * $ J% & ' & J - 
) 2
,

(7)

where the prior covariance matrix ! is given by! = " #1 + W $ A #1 $ W% .

(8)However, this introduced additional hyperparameters, $, in addition to #, that should beestimated from the data. In order to reduce the number of hyperparameters, we assumed
the relation between hyperparameters ($ = %*#, % is a scalar parameter) in the

-61-

simulation (Sato et al., 2004).
Introducing the spatial smoothness constraint as a hard constraint ( J = W ! Z ) is
simpler than the above method and does not degrade performance. This corresponds tothe limit % ! &. The corresponding prior for J is given by


' 1
*
P ( J ! ) " exp ) # J$ % & % J , 
( 2+

! = W " A #1 " W$ .

(9)Furthermore, this is equivalent to the following prior and the forward model, bychanging the variable for J to Z using J = W ! Z & 1)
P ( Z ! ) " exp ( # Z$ % A % Z + 
' 2
*


B = G ! W ! Z .(10)This method was used in this paper and our previous paper (Yoshioka et al., 2008).

-62-

2. Sparse linear regression

In this section, we briefly explain the sparse linear regression used for reconstructing
fingertip position (Nambu et al., 2009; Sato, 2001; Sato et al., 2004; Ting et al., 2005;
Ting et al., 2008). In the following, we only consider regression for horizontal fingertip
position from estimated cortical currents for brevity (regression for vertical position andregression from MEG sensor signals can be done in the same way).
In the current problem, the number of selected currents, M (= N_source), is 1,500,
and the number of time lags, D, is 11. Consequently, the number of total input
dimensions, N, becomes very large: M x D = 16,500. Therefore, we use a Bayesiansparse linear regression with the automatic relevance determination (ARD) prior to
avoid the over-fitting problem and increase generalization ability. However, simple
application of the sparse linear regression causes computational difficulty because it
requires an inverse matrix calculation for the input covariance matrix in each iteration.
Inverse calculation of a large-sized matrix with N = 16,500 is practically intractable,since the inverse calculation requires O (N^3) computational time.
There are several sparse regression methods designed to avoid full inverse
covariance calculation. An incremental sparse method (Tipping and Faul, 2003) adds
weight parameters one by one, which increases the evaluation function (free energy)
most. The Variational Bayesian Least Squares (VBLS) (Ting et al., 2005; Ting et al.,2008) approach updates weight parameters according to the correlation between
residual errors and input variables in the manner of gradient descent by introducing
auxiliary variables, and it does not need inverse covariance calculation. Least Absolute
Shrinkage and Selection Operator (LASSO) (Efron et al., 2004; Tibshirani, 1996) is

-63-

another incremental algorithm using an L1 norm for sparse linear regression.
We evaluated the effectiveness of these methods by using various simulation data.
For medium-sized problems with fewer than 1,000 input dimensions, these methods
demonstrated good performance (not shown). On the other hand, when inputdimensions were over 10,000 and input variables were correlated with each other, theconvergence of these methods became very slow and they exhibited poor generalization
performance (not shown).
Therefore, we use another approach that considers the structure of input variables.
In our problem, we can expect sparse representation for cortical vertices (currentsources) while temporally continuous weight can be expected for the temporaldimension. Accordingly, we introduce a sparse condition only for cortical vertices but
not for the temporal dimension. By introducing an auxiliary variable for each time lag
dimension, we can derive an efficient algorithm that uses the inverse of a spatial
covariance matrix whose dimension is M, which is much smaller than the total inputdimension N.
When training datasets X = { x ( t ) | t = 1 : T } andJ = {J i ( t ) | t = 1 : T, i = 1 : M } are given, the data likelihood is given by$
M D')2 '))&& 1 T $&P ( X | J,W, ! ) ! exp && " ! # && x ( t ) " # # wij J i ( t " j ) ))) ))) , (1)&& 2 t =1 &%)( ))i =1 j =1%(
where weight parameters are represented by W = { wij | i = 1 : M, j = 1 : D } and

the bias term w0 is omitted by assuming that the means of X and J are removed
for notational simplicity. ! represents the inverse of an observation noise variance. By
introducing auxiliary variable z j ( t ) for each time lag dimension, we can rewrite Eq.

-64-(1) asP ( X | J,W, ! ) = ! dZ P ( X | Z, ! ) P ( Z | J,W, ! )%D(*2 (**'' D + 1 T %'P ( X | Z, ! ) " exp '' #! $ '' x ( t ) # $ z j ( t ) *** ***
'
''
*) **
2
j =1
t =1 &
&
)

,

(2)

2(
%
**
'' D + 1 T D! $ $ z j ( t ) # w +j , J ( t # j ) ***P ( Z | J,W, ! ) " exp '' #''**2
t =1 j =1
&
)

(

)

where Z = { z j ( t ) | t = 1 : T, j = 1 : D } . An M-dimensional current vector at timet, J ( t ) = !" J1 ( t ) ,…,J M ( t ) #$ % , and an M-dimensional weight vector for the j-th timelag, w j = !"# w1j ,…, wMj $%& ' , are also defined.
Furthermore, we assume the ARD prior for the spatial dimension:
D
$ 1 M
')
&
P (W | ! ) ! exp && " # !i # wij2 )))&% 2 i =1 j =1 )(MP ( ! ) ! * 1 / !i

,

(3)

i =1

P (") ! 1 / "where precession parameters ! = !" !1 ,…, !M #$ % are only introduced for the corticalspatial dimension and non-informative priors are assumed for ! and ! . According to

the Bayes rule, the posterior probability of W , P (W | X, J ) , can be calculated byusing the joint posterior

P (W,Z, !, ! | X, J ) :-65-

P (W | X, J ) = ! dZ d! d!P (W,Z, !, ! | X, J )P (W,Z, !, ! | X, J ) = P (W,Z, !, !, X | J ) / P ( X | J )P (W,Z, !, !, X | J ) = P ( X | Z, ! ) P ( Z | J,W, ! ) P (W | ! ) P ( ! ) P ( ! ) .

(4)

P ( X | J ) = ! dW dZ d! d! P (W,Z, !, !, X | J )The calculation of marginal probability, P ( X | J ) , cannot be done analytically.
Therefore, we use the variational Bayesian method (Attias, 1999; Sato, 2001). We
introduce a trial distribution Q (W,Z, !, ! ) for the approximation of the jointposterior P (W,Z, !, ! | X, J ) 
By defining the free energy F (Q ) , as
F (Q ) = ! dW dZ d! d! Q (W,Z, !, ! ) log ( P (W,Z, !, !, X | J ) / Q (W,Z, !, ! ) ) ,(5)
the joint posterior P (W,Z, !, ! | X, J ) can be obtained as the maximum of the freeenergy. We also introduce the factorization approximation by assuming

Q (W,Z, !, ! ) = QW (W )QZ ( Z, ! )Q" ( ! ) .

(6)

The maximum of the free energy F (Q ) is obtained by alternately maximizing F (Q )with respect to QZ , QW and Q! .


Z-stepIn the first step, we maximize F (Q ) with respect to QZ while QW and Q! arefixed. The solution is given by

-66-

(

QZ ( Z, ! ) = ! ! | !,Tz j (t ) = wj$T

) # N ( z (t ) | z (t ), ! ! )
"1

t =1

Z

D
*,
$
1 ')
,,
)
% J (t " j ) +
%
J
t
"
k
x
t
"
w
(
)
(
)
&
)
k
,,
D + 1 )(+k =1

1
!Z = eye ( D ) "ones ( D )
D +1
!

"1

,

T '
D
*,2
$
1
))
= & ) x ( t ) " & w j % J ( t " j ) ,,, + D ( D + 1 )Tr ( !W % !J ),+T t =1 )(j =1() represents the normal distribution with a mean z and acovariance matrix ! . ! ( ! | !,T ) represents the gamma distribution with a mean !

where N z | z, !

and a degree of freedom T, and eye(D) and ones(D) represent the D x D identity matrix
and a D x D matrix with one for all components, respectively. !J represents the M x
TM spatial covariance matrix of J ( t ) : !J = 1 / T # J ( t ) ! J ( t )" , and Tr represents
t =1

the matrix trace. The expectation value of the weight vector for the j-th time lag isrepresented by w j , and !W represents the M x M covariance matrix. They arecalculated in the w-step explained below.
w-stepIn the second step, QW is obtained by maximizing the free energy with respect to QW ,while QZ and Q! are fixed. The solution is given by
D(

QW (W ) = ! N w j | w j , !W
j =1

)

T

w j = ( D + 1 ) !!W # J ( t " j ) z j ( t ) ,
t =1

(

"1

( ))

!W = T ( D + 1 ) !!J + diag "( ) represents a diagonal matrix with the diagonal component ! , which

where diag !

-67-

represents the expectation value of the precision parameter ! and is calculated in the! -step explained below. Notice that !W is the same for all time lags j, since thespatial current covariance matrix !J is independent of the time lag j. This reduces thecomputational time considerably.

! -step
Finally, Q! is obtained by maximizing the free energy with respect to Q! , while QZ
and QW are fixed. The solution is given byM(

Q! ( ! ) = " ! !i | !i ,D / 2
i=1

!i

#1

D

2
1
= ( "W ) + $ wijii
D j =1)These three steps are alternatively repeated until the free energy F (Q ) is converged.


-68-

3. Supplementary video showing examples of reconstructed trajectoryThis movie continuously shows eight example trajectories ranging from 0.3 s before to
1.0 s after the movement initiation for a subject. The left panel shows the hand
movement based on the reconstructed trajectory of the fingertip. The upper-right panelshows the reconstructed (red circle) and actual (blue circle) positions that were
projected to the x-y plane as indicated in Fig. 1B. The path of each circle in the past 0.5
s is also shown. The lower-right panel shows the trajectories in the horizontal (x) andvertical (y) directions. This movie is also temporarily located on our web server(http://www.cns.atr.jp/~imamizu/Movie_01.mpg) only for review purposes, and it is notlinked from any other web page.

-69-

4. References
Attias, H., 1999. Inferring parameters and structure of latent variables models by
variational Bayes., Proceedings of 15th Conference on Uncertainty in Artificial
Intelligence, pp. 21-30.
Dale, A.M., Liu, A.K., Fischl, B.R., Buckner, R.L., Belliveau, J.W., Lewine, J.D.,
Halgren, E., 2000. Dynamic statistical parametric mapping: combining fMRI
and MEG for high-resolution imaging of cortical activity. Neuron 26, 55-67.
Dale, A.M., Sereno, M.I., 1993. Improved Localization of Cortical Activity by
Combining EEG and MEG with MRI Cortical Surface Reconstruction: ALinear Approach. J Cogn Neurosci 5, 162-176.
Efron, B., Hastie, T., Johnstone, I., Tibshirani, R., 2004. Least angle regression. The
Annals of Statistics 32, 407-499.
Hamalainen, M., Hari, R., Ilmoniemi, P., Uutila, J., Lounasmaa, O., 1993.
Magnetoencephalography theory, instrumentation, and applications tononinvasive studies of the working human brain. Reviews of Modern Physics
65, 413-497.
Kajihara, S., Ohtani, Y., Goda, N., Tanigawa, M., Ejima, Y., Toyama, K., 2004. Wiener
filter-magnetoencephalography of visual cortical activity. Brain Topogr 17,
13-25.



Nambu, I., Osu, R., Sato, M.A., Ando, S., Kawato, M., Naito, E., 2009. Single-trial
reconstruction of finger-pinch forces from human motor-cortical activation
measured by near-infrared spectroscopy (NIRS). Neuroimage 47, 628-637.
Neal, R.M., 1996. Bayesian learning for neural networks. Springer, New York.
Phillips, C., Rugg, M.D., Friston, K.J., 2002. Anatomically informed basis functions for

EEG source localization: combining functional and anatomical constraints.
NeuroImage 16, 678-695.
Sarvas, J., 1987. Basic mathematical and electromagnetic concepts of the biomagnetic
inverse problem. Phys Med Biol 32, 11-22.

-70-Sato, M.A., 2001. Online Model Selection Based on the Variational Bayes. NeuralComputation 13, 1649–1681.
Sato, M.A., Yoshioka, T., Kajihara, S., Toyama, K., Goda, N., Doya, K., Kawato, M.,
2004. Hierarchical Bayesian estimation for MEG inverse problem. NeuroImage

23, 806-826.
Tibshirani, R., 1996. Regression shrinkage and selection via the lasso. J. Royal. Statist.
Soc B. 58, 267-288.
Ting, J.A., AD’Souza, Yamamoto, K., Yoshioka, T., Hoffman, D.S., Kakei, S., Sergio,
L.E., Kalaska, J.F., Kawato, M., Strick, P.L., Schaal, S., 2005. Predicting EMG

Data from M1 Neurons with Variational Bayesian Least Squares. Advances in
Neural Information Processing Systems.
Ting, J.A., D'Souza, A., Yamamoto, K., Yoshioka, T., Hoffman, D., Kakei, S., Sergio,L., Kalaska, J., Kawato, M., Strick, P., Schaal, S., 2008. Variational Bayesian
least squares: An application to brain-machine interface data. Neural Netw 21,

1112-1131.
Tipping, M.E., Faul, A.C., 2003. Fast marginal likelihood maximization for sparse
Bayesian models. In: Bishop, C.M., Frey, A.C. (Eds.), Proceedings of the Ninth
International Workshop on Artificial Intelligence and Statistics.
Yoshioka, T., Toyama, K., Kawato, M., Yamashita, O., Nishina, S., Yamagishi, N.,

Sato, M.A., 2008. Evaluation of hierarchical Bayesian method throughretinotopic brain activities reconstruction from fMRI and MEG signals.
NeuroImage 42, 1397-1413.

