State-space solutions to the dynamic magnetoencephalography inverse problem using high performance computing

The Annals of Applie d Statistics 2011, V ol. 5, No. 2B, 1207–1228 DOI: 10.1214 /11-A OAS483 c  Institute of Mathematical Statistics , 2 011 ST A TE-SP A C E SOLUTIONS TO THE D YNAMIC MA GNETOENCE PHALOGRAPHY INVERSE PR OBLEM USING HIGH PERF ORMANCE COMPUTING 1 By Christopher J. Long 2,3 , P a t rick L. Purdon 4,5 , Simona Teme reanca 4 , Neil U . Des ai 6 , Ma t ti S. H ¨ am ¨ al ¨ ainen 4,7 and Emer y N. Bro wn 4,5,6,7 2 Imp erial Col le ge L ondo n, 3 GlaxoSmith Kline, 4 Harvar d Me dic al Scho ol, 5 Massachusetts Gener al Hospital, 6 Massachusetts Institute of T e chnolo gy (MIT), 7 Harvar d/MIT Divisi on of He alth Scienc es & T e chnolo g y Determining th e magnitude and lo cation of neural sources within the brain that are responsible for generating magnetoencephalogra- phy (MEG) signals measured on the surfac e of the h ead is a c h al- lenging problem in functional n eu roimaging. The number of p oten- tial sources within the brain exceeds by an order of magnitude the num b er of recording sites . As a consequ ence, the estimates fo r the magnitude and l ocation of the neural sources will b e ill-conditioned b ecause of the underdetermined natu re of the problem. On e wel l- known technique designed to add ress this imbalance is the minimum norm estimator (MNE). This approach imp oses an L 2 regularization constrain t that serves to stabilize and condition the source parame- ter estimates. How ever, these classes of regularizer are static in time and do not consider th e temp oral constraints in herent to the bio- physics of the MEG ex p eriment. In this pap er we prop ose a dynamic state-space mo del that acco unts for b oth spatial and temp oral cor- relations within and across candidate intracorti cal sources. In our mod el, the observ ation mo del is d erived from the steady- state so- lution to Maxw ell’s equ ations while the latent mo del representing neural dy namics is giv en by a random wa lk pro cess. W e show that the Kalman filter (KF) and the Kalman smo other [also known as Received Sep tember 2010; revised April 2011. 1 Supp orted i n part by the National Science F ound ation through T eraGrid resources pro- vided by t he T exas Adv anced Computing Center (T ACC ). This work w as also supp orted by N I H Grants NIBIB R01 EB0522 (E.N.B.), DP2-OD 006454 (P .L.P .), K25-NS 05780 (P .L.P .), DP1-OD003646 (E.N.B.), N CRR P41-RR 14075, R01-EB006385 (E.N.B., P .L.P ., M.S.H.) at Massac husetts General Hospital. Key wor ds and phr ases. Magneto en cephalograph y, source lo calization, Kalman fi lter, fixed interv al smooth er. This is an electro nic reprint of the or iginal article published by the Institute of Mathematica l Statistics in The Annals of App lie d Statistics , 2011, V ol. 5, No. 2B, 1207 –122 8 . This r eprint differs from the o riginal in pagination and typogra phic detail. 1 2 C. J. LON G ET AL. the fixed-interv al smo other (FIS) ] may b e used to solve t h e en su- ing high-d imensional state-estimation problem. U sing a well-kno wn relationship betw een Bay esian estimation and Kalma n filtering, w e sho w that the MNE estima tes carry a significan t zero bias. Calcu- lating these h igh-dimensional state estimates is a comput ationally chal lenging task th at requires High Pe rformance Computing (HPC) resources. T o this end, we employ the NSF T eragrid Sup ercomput ing Netw ork to compu te t h e source estimates. W e demonstrate improv e- ment in p erformance of t he state-space algorithm relative to MNE in analyses of simulated an d actu al somatosensory MEG exp eriments. Our fin dings establish t he b enefits of high-d imensional state-sp ace mod eling as an effective means t o solv e th e MEG source localization problem. 1. In tro d uction. Electromagnetic source imaging is a neuroimaging tec h - nique that p ermits study of neural even ts on a millisecond timescale. Th is t yp e of imaging rev eals brain dynamics that cannot b e seen w ith imaging mo dalities suc h as fun ctional magnetic r esonance imaging (fMRI) or p ositron emission tomograph y (PET) that capture brain activit y on a muc h slo wer timescale. Two of the most imp ortan t examples of electromagnetic imaging are electro encephalography (EEG) and magneto encepholograph y imaging (MEG). These mo dalities acquire multiple time series recorded at locations exterior to the skull and are generated, resp ectiv ely , b y electric and magnetic fields of n euronal cu r rent s within the cortex of the brain [H¨ am¨ al¨ ainen et al. ( 1993 ), Nunez ( 1995 )]. In these exp eriments sub jects execute a task that putativ ely activ ates m u ltiple brain areas. In the case of EEG recordings, the mo de of acquisition inv olv es affixing an arra y of electrod es to the su rface of the scalp an d measur in g the p oten tial differences relativ e to a reference no de. In contrast , MEG records extremely w eak m agnetic fields emanating fr om the brain on the order of femtoT e s la, that is, 10 − 15 T esla [Nun ez ( 1995 )]. T o put these quan tities into co n text, the magnetic field of the ea rth is around 5 × 10 − 5 T esla while that of a present- da y magnetic resonance scanner used for medical diagnostic imaging ranges b et wee n 1 . 5 and 7 T esla. Because the magnetic fields of the br ain are extremely weak, a sp ecially shielded ro om and highly sensitiv e arra y of detectors, called sup ercondu cting quan- tum interference devices (SQUIDs), are requir ed to u ndertak e the MEG recordings. Generating images of the brain’s cortical activit y from MEG time-series recordings r equires solving a h igh-dimensional ill-p osed in verse problem. In the case of MEG imaging, the nonunique nature of the inv erse problem arises b ecause th e num b er of p oten tial three-dimensional b rain sources (as a con- sequence of the biophysics of the problem) is infin ite, while the measured t wo- dimensional S QUID arr ay t ypically con tain 306 magnetomete r and gra- diometer detectors. DYNAMIC ST A TE-SP ACE METHODS F O R MEG 3 The accuracy of the inv erse solution dep ends critically on t w o main fea- tures of the p roblem, sp ecifically , the precision of the forward mo d el and the c hoice of in verse solution algorithm. The MEG forw ard mo del is d e- signed to describ e ho w activit y propagates from sour ces in the cortex to the SQUID d etectors. These forward mo dels are n ormally constructed by com binin g information ab out head and br ain geometry with th e electromag- netic p rop erties of the skull, du r a and scalp to yield a n umerical appr oxima- tion to Maxw ell’s quasistatic equations; see, for example, H¨ am¨ al¨ ainen et al. ( 1993 ). Und er this forw ard mo d el the stren gth of the magnetic field recorded at a giv en S Q UID d etector n ear the surface of th e head is appro ximately prop ortional to the recipro cal of its squared distance from a giv en corti- cal source. The seco n d co nsideration w e face relates to the c hoice of s ou r ce lo calizat ion algorithm. Solving this localization problem is one of the most tec hnically c hallenging in MEG imaging researc h and h as b een an activ e area of metho d s dev elopment for b oth EEG and MEG o v er the past t wo decades. Historically , t wo of the most p opular metho ds f or solving the EEG/MEG in verse problem h a ve b een the dip ole-based and (l inear) distributed-source metho ds. Dip ole-based metho d s mo del the u nderlying neuronal sou r ces as discrete current dip oles whose lo cation, orien tation and amplitude are unkn own quan- tities to b e estimated [Mosher, Lewis and L eah y ( 1992 ), Uutela, H¨ am¨ al¨ ainen and Salmelin 1998 ]. Cho osing the n u m b er of d ip oles to include in the m o del is problematic. Eigen v alue decomp osition metho ds ha ve b een deve lop ed to address this mod el selection issue [Mosher, Lewis and Leah y ( 1992 )]. Ho w- ev er, these app roac hes still require su b jectiv e inte rpretation when the eigen- v alue distribution do es n ot yield an ob vious distinction b et ween the signal and noise subsp aces. In addition, find ing the b est-fitting parameters for the m u ltidip ole m o del requir es nonlinear optimization. Since the measur ed fields dep end n onlinearly on the dip ole p osition parameters, t yp ical min imization routines ma y not yield the globally optima l estimates for these p arameters. Prop osed solutions to this pr oblem include the MUSIC (MUltiple SIgnal Classification) algorithm [Mosher and Leah y ( 1998 )], global optimizatio n algorithms and partly h eu ristic mo d el constructions that use in teractiv e soft- w are to ols [Uutela, H¨ am¨ al¨ ainen and Salmelin ( 1998 )]. Ho wev er, n on e of these metho ds u tilize the temp oral sequence of the s ignals, that is, if the original data p oint s are fi rst p ermuted and an in verse p ermutation is applied to the source estimat es, the r esults are the same as without p erm utation. Finally , in man y im p ortant clinical and neuroscien tific app lications, suc h as epilepsy , sleep or general anesthesia, the dip ole m o del is inappropr iate, since the gen- erating currents are m ost lik ely distrib uted across large areas on th e cortex. In th e d istributed source mo dels, eac h lo cation in th e b rain r epresen ts a p ossible site of a curr en t source. S ince the n umb er of u nkno wn sources v astly outnum b ers the num b er of EEG or MEG data r ecordings, constraint s 4 C. J. LON G ET AL. are required to obtain a unique solution. Th e m inim u m-norm estimator (MNE) emplo ys an L 2 -norm data “fidelit y” term to quantify the r elationship b et w een the obser ved magnetic field recordings and the estimated sour ce p re- dictions and an L 2 -norm regularization or p enalt y term on th e magnitude of those solutions. While the m inim u m-norm estimate (MNE) yields the so- lution with th e smallest energy (su m of s q u ares) across th e solution space [H¨ am¨ al¨ ainen et al. ( 1993 )], it tends to consistently fa v or lo w-amplitud e so- lutions th at are lo cated close to the scalp sur f ace. The relativ e cont ribution of the data term and sou r ce term can b e con tr olled or tun ed through a r eg- ularization parameter. In MNE, the data term is sp ecifically the L 2 -norm of spatially w h itened data while th e source term is the L 2 -norm of the current amplitudes. A we ll-kno wn v arian t of this approac h is the L O RET A (LOw- Resolution Electromagnetic T omography Algorithm), where the L 2 -norm source term sp ecified in MNE is r eplaced by a w eigh ted norm that appro x- imates the spatial Laplacian of the curren ts [P ascual-Marqui, Mic hel and Lehmann ( 1994 )]. I n the case of minimum-curren t estimates (MCE) [Uutela K. H¨ am¨ al¨ ainen and Somersalo ( 1999 )], the source term is tak en to b e the L 1 - norm of th e source cur ren t amplitud es. A consequence of this c hoice is that the estimates will b e sparse as opp osed to diffuse, leading to families of solu- tions that may closely r esem ble th ose un co v ered th r ough a dip olar analysis. Beamformer appr oac hes, originally inspir ed b y problems in radar and so- nar arr a y signal pro cessing [Krim and Vib erg ( 199 6 )], seek to lo calize EEG/ MEG activit y b y sp ecifying a series of sp atial filters, eac h tun ed to maxi- mize sensitivit y to a particular spatial location within the brain [V an V een et al. ( 1997 )]. These metho ds assume, ho wev er, that EEG/MEG sources are sp atially uncorr elated, and are limited by in terference crosstalk from sources outside the fo cus of th e b eamformer spatial fi lter. A more recen t related approac h by Zumer et al. ( 2007 ) features a graphical mod el frame- w ork to c h aracterize the probabilistic relationships among latent and ob- serv able temp oral sour ces for a verag ed ERP data. On ce a lo w-d imensional statistica l data subsp ace is estimated, the algorithm compu tes full p osterior distributions of the parameters and hyperp arameters with the primary ob- jectiv e of alleviating temp oral correlation b et ween spatially distinct sources (source separation). In later work [Zu m er et al. ( 2008 )] this approac h is extended somewhat to include explicit temp oral mo deling b y dr a wing on estimated w eigh tings from a predefined family of smo oth basis f u nctions. Wipf et al. ( 2010 ) explores an alternate pr o cedure to empirical Ba yesian source mod eling that attempts to s imultaneously captur e and model corre- lated source-space dynamics in the presence of u nkno wn dip ole orientat ion and bac kground interference. Ba y esian approac hes to MEG sour ce lo calizati on ha ve included Phillips et al. ( 2005 ), Mattout et al. ( 2006 ) and F r iston et al. ( 2008 ). Collectiv ely , the- se tec hniques constr u ct the source lo calization problem within an in stan ta- neous empirical Ba y es framewo rk. Th ese approac hes sp ecify a nested and DYNAMIC ST A TE-SP ACE METHODS F O R MEG 5 hierarc h ical series of spatial linear mo dels to b e s u bsequently solv ed un - der a p enalized Restricted Maximum Lik eliho o d (ReML) pro cedure. Th e m u ltiresolutional natur e of this regime carries w ith it a measure of spatial adaptation for th e estimated h yp erparameters since the s p atial clusters of activ ated dip oles are amalgamat ed across sev er al (orthogonal) scales. A t the same time the most lik ely set of spatial priors for eac h comp eting mo del is c hosen usin g information selectiv e mo del rankin g pro cedu res. Th e lat- ter approac h [F riston et al. ( 2008 )] extends the use of ReML to engender a framework aiming to m ov e from p arametric spatial pr iors to larger sets of imputed sour ces with lo calized b ut not n ecessarily contin uous supp ort. The principle adv an tage of this extension is its in trin sic abilit y to captur e source solutions ranging from sparse dip ole mo d els at one end of the lo calization sp ectrum to dense but smooth distrib uted co nfigurations of spatial dip oles at the ot her. State-space m etho ds featuring laten t temp oral mo dels ha ve recen tly b een prop osed. F or example, studies by Galk a, Y amashita and Ozaki ( 2004 ) and Y amashita et al. ( 2004 ) ha ve used a rand om-w alk dyn amical mod el w ith Laplacian spatial constrain ts to r ep resen t the d ynamics of EEG sour ce cur- ren ts, emplo ying the Kalman filter and r ecursiv e least-squares framework to p erform source localization. Daunizeau and F riston ( 2007 ) used a state - space mo del to capture laten t neuronal d ynamics by placing a first-order autorgressiv e (AR) mo del within a fu ll spatiotemp oral v ariational Ba yes pro cedur e. Eac h ensemble of sources is c haracterized by establishing an ef- fectiv e region of influence in whic h the hidden neural state-dynamics are assumed constan t. T his leads to a temp oral representati on of eac h region in terms of a sin gle a verag e dynamic. Ou, H¨ am¨ al¨ ainen and Golland ( 2009 ) ha ve devel op ed a m ixed L 1 L 2 -norm p enalized mo d el to estimate MEG/EEG sources. This metho dology allo ws for sp atially sparse solutions while ensuring physiologi cally plausible and n u merically stable time course reconstructions through the use of basis fun c- tions. Since estimation of this kind of m ixed-norm regularized p roblems rep- resen ts a compu tational c hallenge, conv ex-cone optimization metho ds are implemen ted to efficien tly compute the MEG/EEG sour ce estimates. Th ese pro cedur es are shown to yield s ignifican t improv emen t ov er conv en tional MNE appr oac hes in b oth simulate d and actual MEG data. In the ma j ority of MEG stim u lus resp onse exp eriments, the ob jectiv e is to estimate th e activit y in a give n brain region based on m u ltiple rep etitions of the same s timulus. The time int erv al of the MEG recordings after the onset of the stim u lus is usu ally aroun d 1,000 msec. C on s equen tly , the optimal estimate of the source acti vit y at a given time point shou ld b e based on all the d ata recorded in this time w indo w , as opp osed to only the recordings within a small in stan taneous sub in terv al as is the case in the MNE algorithm. 6 C. J. LON G ET AL. Th us, in th e curr en t w ork, we mo del the spatiotemp oral structure of the MEG source lo calization problem as a full-rank state-space estimat ion pro- cedure. In co ntrast to p r evious related app roac hes, for example, Gal k a, Y a- mashita and Ozaki ( 2004 ), Y amashita et al. ( 2004 ) and Ou, H¨ am¨ al¨ ainen and Golland ( 2009 ), we reco ver the un d erlying neur al state dynamics across th e en tire solution space usin g all information a v ailable in the MEG measure- men ts to compute a new source estimate at eac h lo cation an d time p oint . W e do not restrict our state-co v ariance structure to op erate eit her through reduced-rank appr o ximations or th rough “small-v olume” state-space mo d - els that op erate w ith in lo calized spatial n eighb orh o o ds [Galk a, Y amashita and Ozaki ( 2004 ), Y amashita et al. ( 2004 )]. Th e analysis computes activit y at a p articular sour ce with activit y at that source com bin ed with (p oten- tially all) neigh b oring sources in b oth preceding and subsequent time p erio ds (spatiotemp oral correlation). W e find that solution of the resulting high- dimensional s tate-space mod el requires the ap p lication of high p erformance computing (HPC) resources. W e show th at t w o solutions that are naturally suited to this type of dy- namic inv erse pr oblem are the K alman filter and the fixed -in terv al sm o other [Kitaga w a and Gersc h ( 1996 ), Kay ( 199 3 )]. These algorithms are based on a series of mathematical relationships that emp lo y prin ciples from electro- magnetic field theory to relate the observe d MEG m easuremen ts to th e underlying dynamics of the cu r rent source dip oles. In con trast to p opu - lar sour ce localization method s such as MNE where the regularizat ion co n- strain ts are formulate d in terms of a fixed-error co v ariance matrix, the state- space estimators feature time-v arying err or co v ariances and pr op agate past (and futur e) information into the cur r en t up date. W e describ e how these algorithms m a y b e used to p erform this state estimation by condu cting th e computation on the NSF T eragrid Sup ercom- puting Net w ork. In Section 3 we detail the necessary s up ercompu ting meth- o ds needed to im p lemen t the high-dimensional state-space estimation. W e demonstrate the impro v ement in p erformance of th e state-space algorithm relativ e to MNE in analyses of a sim u lated MEG exp eriment and actual somatosensory MEG experiment. 2. Theory . W e assume that in an MEG exp eriment a stim ulu s is ap- plied T times and at eac h time for a (short) window of time follo wing the application of the stim u lu s MEG activit y is measur ed sim u ltaneously in S SQUID recording sites. Let Y s,r ( k ) d enote the measurement at time k at lo- cation s on rep etition r w here k = 1 , . . . , N , s = 1 , . . . , S and r = 1 , . . . , T . W e tak e Y s ( k ) = T − 1 P T r =1 Y s,r ( k ) and we define Y k = Y ( Y 1 ( k ) , . . . , Y S ( k )) to b e the S × 1 v ector of measur emen ts recorded at time k . W e assume that there are P current sour ces in the b rain and that the relation b etw een the MEG DYNAMIC ST A TE-SP ACE METHODS F O R MEG 7 recordings and the s ou r ces is d efined by Y k = H J k + ε k , (2.1) where J k = ( J k , 1 , . . . , J k ,P ) is the 3 P × 1 vec tor of source activities at time k , eac h J k ,i is a 3 × 1 vecto r, and H is the S × 3 P lead field matrix, deriv ed as describ ed in H¨ am¨ al¨ ainen et al. ( 1993 ). This matrix is approxima ted in practice using the known conductivit y profile of an MRI-deriv ed T1 im- age to der ive a one-la yer h ead mo del using the b ound ary elemen t metho d [H¨ am¨ al¨ ainen and Sarv as ( 1989 )]. F urtherm ore, ε k is a zero mean Gauss ian noise with co v ariance matrix Σ ε . T o build on an id ea originally suggested in the EEG literature [Galk a, Y amashita and Ozaki ( 2004 ), Y amashita et al. ( 2004 )], w e assume that the J k b ehav e lik e a sto c hastic firs t-order pro cess with a p oten tially w eigh ted neigh b orho o d (six direct neighbors in the case b elo w ) constraint defined at v o xel p b y J k ,p = J k − 1 ,p + 1 6 m − 1 ℓ X ℓ ∈ M ( ℓ ) J k − 1 ,ℓ + v k , (2.2) where M ( ℓ ) is some small neighb orh o o d aroun d v oxel p and m ℓ is a mask of (estimated) autoregressive w eigh ts, and v k is a 3 P × 1 v ector of Gaussian noise with mean zero and co v ariance matrix Σ w . The linear str u cture in equation ( 2.2 ) means that we can generalize this structure into matrix format suc h that all sources in the brain (i.e., not only lo cal relationships) are encompassed in the mo del. Th at is, J k = F J k − 1 + v k . (2.3) If F is an iden tity matrix, the state dynamics reduce to a random walk pro cess. Equations ( 2.1 ) through ( 2.3 ) defi ne the observ ation and laten t time-series equations, r esp ectiv ely , for a state-space mo del form u lation that can b e used to pro vid e a dyn amic description of the MEG sour ce lo calizatio n problem. 2.1. Standar d ME G inverse solution. Th e standard approac h to the MEG source lo calization prob lem is to solv e an L 2 regularized least-squares prob- lem at ea ch time k . T hat is, k Y k − H J k k 2 Σ ε + k J k − µ k 2 C , (2.4) where k y − x k 2 Q denotes the Mahalanobis distance b et ween th e v ectors y and x with err or co v ariance matrix Q , µ is an offset p arameter (pr ior mean) and C is th e regularization co v ariance matrix [H¨ am¨ al¨ ainen and Ilmoniemi ( 1994 )]. If w e tak e µ = 0 and define the s ou r ce co v ariance matrix as C = λR , where R is a d iagonal, scaled m atrix normally compu ted in adv ance, and 8 C. J. LON G ET AL. if λ > 0 , th en at eac h time k an in stan taneous MEG sou r ce estimate (the MNE solution) is giv en as J MNE k = λR H T ( λ HRH T + Σ ε ) − 1 Y k (2.5) for k = 1 , . . . , N . The MNE estimate is a lo cal Ba y es’ estimate b ecause it uses in its compu tation only the data Y k at time k and the Ga ussian prior distribution with mean µ = 0 and co v ariance matrix λR . Hence, th e MNE estimation p ro cedure imp oses n o temp oral constrain t on the sequ ence of solutions. 2.2. Kalman filter solution. Giv en the s tate-space formulatio n of the MEG observ ation pro cess in equations ( 2.1 ) and ( 2.3 ), it follo ws that the optimal estimate of the current source at time k us in g the data Y 1 , . . . , Y N up through time N is giv en b y the Kalman filter [Kitaga wa and Gersc h ( 1996 )] J k | k − 1 = F J k − 1 | k − 1 , (2.6) W k | k − 1 = W k − 1 | k − 1 + Σ w , (2.7) K k = W k | k − 1 H T [ H W k | k − 1 H T + Σ ε ] − 1 , (2.8) J k | k = W k | k − 1 H T [ H W k | k − 1 H T + Σ ε ] − 1 Y k (2.9) + [ I − K k H ] J k | k − 1 , whic h simplifies to J k | k = J k | k − 1 + K k [ Y k − H J k | k − 1 ] , (2.10) W k | k = [ I − K k H ] W k | k − 1 (2.11) for k = 1 , . . . , N , give n initial conditions J 0 ∼ N (0 , W 0 | 0 ) , Σ w = W 0 | 0 . At eac h time k the K alman filter compu tes p ( x k | Y 1 , . . . , Y k ), wh ic h is the Gaussian distribution with mean J k | k and co v ariance matrix W k | k . In terms of the reg- ularization criterion f unction in equation ( 2.4 ), the K alman filter solution is equiv alen t to c h o osing at time k , µ = [ I − K k H ] F J k − 1 | k − 1 and C = W k | k − 1 , where we ha ve us ed th e matrix in version lemma to re-express the Kalman filter up date in equation ( 2.9 ) in order that w e can compare it directly with the MNE solution in equation ( 2.5 ). F rom this comparison w e see that the Kalman solution impro ves up on the MNE s olution in t w o im p ortant wa ys. First, in the Kalman solution, the offset parameter and the regularization matrix are different at eac h time k and are giv en , resp ectiv ely , by µ and the one-step pr ediction error co v ariance m atrix W k | k − 1 . The MNE solution at eac h time k has a fixed prior mean µ = 0 and (temp orally) fi xed regular- ization matrix C = λR . Because of this choic e of regularization constrain t, equation ( 2.9 ) sh ows that the Kalman filter estimate at time k is a linear DYNAMIC ST A TE-SP ACE METHODS F O R MEG 9 com bination of J k − 1 | k − 1 , the cu r rent s ource estimate at time k − 1, and Y k , the observ ations at time k . In con trast, the MNE solution at eac h time k is a linear combinatio n of the observ ed d ata and a fixed prior mean of µ = 0. In this regard, the MNE estimate biases the s olution at eac h time k to ward 0. T ak en together, these observ ations sh ow that the sto c h astic con tin uity as- sumption in the Kalman state mo del results in a time-v arying constraint up on the fluctuating source estimates. 2.3. Fixe d-interval smo othing algorithm solution. Because th e MEG time series are often fixed length recordings, w e can go b eyo n d the estimate J k | k pro v id ed b y the Kalman up date and compute the p osterior density p ( x k | Y 1 , . . . , Y N ) at time k giv en all th e data in the exp eriment. T o do so, w e com bine th e Kalman filter with the fi x ed -in terv al smo othing algorithm (FIS) [Kitaga w a and Gersc h ( 1996 ), Ka y ( 1993 )], that may b e computed as follo ws: A k = W k | k W − 1 k +1 | k , (2.12) J k | N = x k | k + A k [ J k | k − J k +1 | k ] , (2. 13) W k | N = W k | k + A [ W k | N − W k +1 | k ] A T k (2.14) for k = N − 1 , . . . , 1 and initial conditions J N | N and W N | N computed from the last step of th e Kalman filter N . It is w ell kno wn that p ( x k | Y 1 , . . . , Y N ) is a Gaussian distribu tion with mean J k | N and co v ariance matrix W k | N [Kitaga w a an d Ge rsc h ( 1996 ), Kay ( 1993 )]. 2.4. Pr actic al c onsider ations. Implemen tation of the Kalman filter first requires making estimates of b oth the initial state, the s tate co v ariance ma- trix Σ w and the error co v ariance matrix Σ ε . W e estimate d the noise co v ari- ance matrix Σ ε from measurement s tak en in th e scanner in the absen ce of a sub ject. W e next estimated the initial state as the MNE solution J MNE 0 and the state co v ariance matrix by firs t computing th e MNE source estimates J MNE 1 , . . . , J MNE N , subsequ en tly deriving Σ w from the sample co v ariance ma- trix u sing a differenced sequence of these sta tic estimates. 3. Sup ercomputer implementat ion. Historically the Kalman filter has found widespread u se in sev eral high-dimensional mo deling domains, in - cluding weather forecasting [F arrell and Ioannou ( 2001 )] and oceanograph y [F uku mori et al. ( 1993 ), F u kumori and Malanotte -Rizzoli ( 1995 )]. Kalman filters and fixed-inte rv al smo others are adv ant ageous in these scenarios as, under assumptions of linearit y and normalit y , they are (near) optimal esti- mators. In add ition, their desirable prop erties hold across a wide v ariet y of time-v aryin g linear (and n onlinear) m o dels. Ho wev er, in its standard form 10 C. J. LON G ET AL. the Kalman fi lter is compu tationally prohibitiv e for these classes of prob- lems. In the example applications listed ab o ve, the n umerical calculati ons are often carried out on systems of state dimension N ∼ O (10 7 ) with state co v ariance matrices of size N 2 ∼ O (10 14 ) [F arrell and Ioannou ( 2001 )]. Com- putationally , the most intense asp ects of the Kalman algorithm stem from the p rediction u p date in the co v ariance matrices th at requir e costly linear algebraic up dates at ea c h time step. Since the dynamical error s tr ucture o f these systems is often well und ersto o d, many n umerical solutions to these kinds of paradigms, for exa mple, F arrell and Ioannou ( 2001 ) and F uku m ori and Malanotte-Riz zoli ( 1995 ), emplo y a r ange of mo del-redu ction techniques in their form ulations to ac hiev e computational trac tabilit y . In the case of the MEG and EEG inv erse problem, the solution sp ace is ∼ O (10 3 − 10 4 ), le ading to error co v ariances of dimension P 2 ∼ O (10 6 –10 8 ). Th us, the compu tational problem is n ot quite s o intensiv e as in the fore- casting app lications, meaning that we can feasibly emplo y full-rank state- estimation metho ds to compu te th eir solution on an HPC system. When p erformin g Kalman filtering, th e computations dominating eac h time s tep in volv e three high-dimensional full-rank matrix multiplica tions, leading to a total app ro ximate computational cost of 3 P 3 (excluding auxiliary lo we r di- mensional linear op erations). When taking in to accoun t an cillary v ariables, the amount of memory requir ed is at least 24+ gigab ytes for eac h up date op eration. In addition, th e FIS r equires one add itional such m ultiplication at eac h time step follo we d by a full-rank P × P m atrix inv ersion (2 P 3 ). Also, FIS r equires appro ximately three gigab ytes of storage sp ace p er time p oin t in order to sav e the prediction and up date co v ariances (in 32-bit storage format) estimated du ring the f orw ard pass of the Kalman filter. The scale of these compu tations ren d ers th em infeasible on nearly an y state-of-the-art standalone computing resource. T o addr ess this limitation, we parallellized the Kalman and FIS filtering computations suc h that the data-in tensive parts of the algorithm w ere d istributed across m u ltiple no des of a High P er- formance Computin g system [Blac kford et al. ( 1997 )]. F or this pur p ose w e utilized the NS F T eragrid resource at the T A C C (T exas Adv anced Com- puting Cen ter). Th is resource co mprised a 1024-processor Cray/D ell Xeon- based Linux cluster with a total of 6.4 T eraflops co mputing ca pacit y . 4. Results. 4.1. Simulate d MEG exp eriment. W e designed a set of sim ulation stud- ies to compare the p erformance of the dynamic lo calizatio n metho ds against MNE. Prior to constructing the dynamic simulation, we first computed th e spatial conductance pr ofile across the head and sp ecified the spatial r eso- lution of the discretized solution space (source lo cations). T o r estrict this source space to the cortical su r face, we emp loy ed anatomical MRI data DYNAMIC ST A TE-SP ACE METHODS F O R MEG 11 obtained from a sin gle sub ject with a high-resolution T1-w eight ed 3D se- quence (TR/TE/flip = 2,530 ms/3.49 ms/7 ◦ , p artition thic kness = 1.33 mm , matrix = 256 × 256 × 128, field of view = 21 cm × 21 cm) in a 3-T MRI scan- ner (Siemens Medical Solutions, Erlangen, German y). The geometry of the gra y–white matter surf ace was subsequen tly computed u sing an automatic segmen tation algo rithm to yield a triangulated mo del with ap p ro ximately 340,00 0 vertic es [Dale, Fisc hl and Ser en o ( 1999 )]. Finally , we utilized the top ology of a recurs ively sub d ivided icosahedron with appr o ximately 5 mm spacing b et w een the sour ce n o des to giv e a cortical sampling of appro xi- mately 10 4 lo cations across b oth hemisph eres. Using MNE, we c hose as the region of in terest a section of the left h emi- sphere ov er the p rimary somatosensory and motor cortices. W e computed a single la yer homogeneous forward mo del or lead field matrix H o ver all the sampled vo xels in the left and right hemisp h eres. The region of in terest con tained 125 activ e v oxels f rom the more than 10,000 gra y matter vo xels that could b e p oten tial sources f or an observ ed magnetic fi eld un der this parcellation. F or this simulatio n , therefore, the num b er of activ e sources w as P activ e = 125, and the num b er of measurement channels wa s S = 306. W e n ext constructed H with dimension 306 × 5,120 (left-hemisphere), cor- resp ond ing to a spatial sampling of 5 mm. Note that we c hose to estimate only the d ominan t n ormal MEG comp onent at eac h ve rtex (i.e., z -direction) as opp osed to esti mating th e triplet of x, y and z con tributions from which b oth magnitude and directional information could b e resolv ed. In eac h of the c h osen 125 v oxel s in the sour ce sp ace, we generated an “activ ation” signal consisting of a mixture of 10 and 20 Hz f requencies mo dulated by a 0.4 Hz en v elop e to sim u late typical signal patterns commonly observ ed in the motor and somatosensory cortices; see, for example, Lounasmaa et al. ( 1995 ). Next w e ad d ed Gaussian noise to all P = 5,120 v ertices to generate an image-wide signal-to-noise ratio (SNR) that ranged b etw een 0.1 and 2 (in line with typical SNRs encountered in MEG studies). Sp ecifically , w e defined SNR as k H J 2 k S P σ 2 . (4.1) When computing the s ou r ce reconstructions w e set the observ ation co v ari- ance matrix as the empt y ro om cov ariance, that is, that wh ic h is generated from bac kground noise measurements tak en from the system while the sub- ject is absent from the scanner. W e computed three in v erse solutions: the MNE solution, the KF solution, and the FIS solution. In this s im u lation study , w e r eco v ered the source estimate for eac h method an d for eac h v alue of SNR using three different c hoices of tun in g parameter ( λ = 0 . 5 , 1 , 3). The time-course length w as 12 0 time p oin ts and with an assumed sampling rate of 600 Hz, equated to a data seg men t of ab out 200 msec. 12 C. J. LON G ET AL. T o examine the b enefits of the dyn amic state-space pro cedu res in relation to the MNE sour ce localization algorithm, we computed the Kalman and Fixed-In terv al Sm o othing fi lters, generating estimates of the source lo cations and amplitudes with in the en tire cortex. When applied to the 200 ms window of MEG data, th is analysis a v eraged ab out one hour for eac h simulated record when d istributed across 16 of the CRA Y-DELL cores. F or eac h c h oice of SNR (thr ough a swee p of 20 SNRs c h oices equisp aced b et ween 0.1 and 2), the FIS algorithm to ok around t wo hours on 24 C RA Y-DELL computational no des, leading to a total CP U time of around 1,280 hours for the wh ole sim u lation. Storage of the prediction and up date co v ariances for eac h SNR required approximate ly 110 gigab ytes of disk space, r esu lting in a total of 2 T erab ytes of storage for eac h of the th ree c h oices of tuning parameter ( λ ), that is, 6 T erabytes in total. Figure 1 illustrates the spatial extent of the s imulated signal (a) tak en as a snapshot at the p eak of the damp ed sinusoidal signal. Y ello w to r ed color map v alues indicate progressiv ely smaller current sour ces ( J ) represent ing ground truth, panel (a), or th e reconstructed estimates (b)–(d). Figure 1 panels (b)–(d) show the estimates as computed with (b) the min im u m norm estimate, (c) th e Kalman filter and (d) the Fixed In terv al Smo other. F or eac h metho d, the SNR as defined by equation ( 4.1 ) and the tuning parameter go v erning the strength of the estimate smo othness were set to SNR = λ = 1. Figure 2 depicts temp oral reconstructions f or the th ree metho ds with the c hoice of SNR = λ = 1. F or eac h region A, B or C, the time series at the v oxe l corresp onding to the p eak MNE r esp onse w as plotted and o verlaid on to the “ground -truth” simulate d signal. The state-space reconstructions sho w considerably less v ariabilit y than the static MNE time course, ho w ever, the KF appr oac h disp la ys elev ated b ias as compared to the FIS metho d. F ur thermore, the KF metho d con tains a temp oral shift in relation to the sim u lated s ignal that is n ot apparent in MNE and FIS. These fu ndamenta l differences in th e temp oral b eha vior of the three methods are consistent in space, across S NR, and are in v arian t to choic e of λ . Figures 3 and 4 illustrate the b ehavior of MSE (Mean-Squared Error) as a function of SNR for the three c hoices of tu ning parameter: panel (a) λ = 0 . 5, panel (b) λ = 1 , and panel (c) λ = 3 b et w een the tru e s ignal and the estimated signal for eac h estimatio n metho d. In Figure 3 MSE is calculated at all v ertices in th e r econstructed cortical surface an d av eraged sp atially for eac h c h oice of SNR. The fi xed-in terv al smo other shows the lo west MSE w hic h indicates that the estimated source cu rv es are closer to the ground tru th sim u lations compared to the alternativ e metho ds. By con trast, the MNE appro ximations sh o w the least f av orable reco very of the sim ulated sources. In Figure 4 w e computed the MSE at eac h v ertex for eac h source r econstruction metho d in the same m an n er but calculate d the a v erage MSE only w ithin the three activ ated regions. These plots, th erefore, examine the b eha vior DYNAMIC ST A TE-SP ACE METHODS F O R MEG 13 Fig. 1. A ctivation maps in thr e e r e gions, A, B and C at the p e ak r esp onse time of 126 ms for: (a) T rue si gnal; and sour c e estimates c ompute d by (b) MNE; (c) the Kalman filter; and (d) the Fixe d Interval Smo other. These images wer e c ompute d using a r e gularization p ar ameter ( λ = 1 ). of the source reconstructions by zo oming into areas that con tain signal, th u s decouplin g the m o del fi ts to areas of true signal fr om those vertice s con taining bac kground noise. Th e three metho d s show similar b ehavio r as in Figure 3 in the act iv ated regio ns, with the fixed-in terv al smo other again sho w in g sup erior p erforman ce to the other t wo m etho ds. As a consequence of th e r elativ e MSE errors exhibiting s ignifican tly smaller v alues in these signal-ric h regions, we displa y their effects on a log linear scale to b etter delineate the salien t differences b et ween metho ds. 5. Somatosensory MEG exp eriment. T o illustrate the state-space and MNE methods on an actual MEG experiment, w e estimat ed the sources of the alpha (8–13 Hz) and mu (10 Hz and 20 Hz) sp onta n eous b rain activit y 14 C. J. LON G ET AL. Fig. 2. Time c ourse of sour c e estimates fr om the simulate d M EG exp eriment for the maximal ly r esp onding voxel r elative to the MNE solution, wher e e ach p anel (fr om top to b ottom) c orr esp onds dir e ctly to e ach of the thr e e r e gions, A, B and C denote d i n Fi gur e 1 (a) . Gr e en denotes the MNE solution, r e d the KF, dark blue is the FIS solution, whi le cyan r epr esents the true signal. DYNAMIC ST A TE-SP ACE METHODS F O R MEG 15 Fig. 3. Aver age p er c entage r elative change in MSE c ompute d as a f unction of SNR fr om the simulate d MEG exp eriment for al l voxels on the c ortic al surfac e for e ach lo c alization metho d: MNE (li ght gr ay), Kalman filter (dark gr ay), and Fixe d Interval Smo other (char- c o al). (a) λ = 0 . 5 , (b) λ = 1 , (c) λ = 3 . 16 C. J. LON G ET AL. Fig. 4. Av er age p er c entage r elative change i n MSE c om pute d as a function of SNR in the simulate d MEG exp eriment acr oss al l active sour c es f or e ach lo c al i zation metho d: MNE (light gr ay), Kalman filter (dark gr ay), and Fixe d I nterval Smo other (char c o al). (a) λ = 0 . 5 , (b) λ = (c) , λ = 3 . DYNAMIC ST A TE-SP ACE METHODS F O R MEG 17 using anatomicall y-constrained whole-head MEG. Syn c hr onous co rtical rhyt hmic activit y of large p opu lations of neurons is asso ciated with distinct brain states and has b een the sub ject of extensiv e inv estigati on u sing unit electroph ysiology , as well as EEG, MEG and fMRI tec hn iques [Lounasmaa et al. ( 1995 ), Salmelin and Hari ( 1994 )]. The alpha rhythm is r ecorded ov er the p osterior p arts of the br ain, b eing strongest when the sub ject h as closed ey es and reduced wh en the sub ject h as open ey es. Sour ces of alpha rhythm ha ve b een iden tified with MEG and EEG to lie in th e parieto o ccipital su lcus and calcarine fissu r e [F reeman, Ah lfors and Menon ( 2009 )]. The m u rhythm represent s activit y close to 10 and 20 Hz in the somatomotor cortex, and is kno w n to reflect resting states c haracterized b y lac k of mo vemen t and somatosensory in put. The 20 Hz comp onent tends to b e lo calized in the mo- tor cortex, ante rior to lo cation of the 10 Hz comp onen t [Lounasmaa et al. ( 1995 )], b eing sp ecifically suppressed or elev ated at top ographic lo cations represent ing moving or r esting ind ivid ual lim bs, r esp ectiv ely . By compari- son, th e 10 Hz comp onent arises in th e somatosensory cortex, and is thought to reflect the absence of sensory input from the upp er lim bs [Liljestr¨ om et al. ( 2005 ), Salmelin et al. ( 2005 )]. F ollo win g approv al b y the Massac husett s General Hospital Hu man Re- searc h Committee, the MEG signals were acquired in a h ealth y male su b- ject in a m agnetical ly and electrically shielded r o om at the Martinos Center for Biomedical Imaging at the Massac husetts General Hosp ital. MEG signals w ere recorded from th e en tire head using a 306-c hannel d c-SQUID Neuromag V ectorview system (Elekta Neuromag, Elekta Oy , Helsinki, Finland) while the sub ject sat with his head inside the helmet-shap ed dewar con taining the sensors. The mag netic fields w ere recorded simultaneously at 102 locations, eac h with 2 planar gradiometers and 1 magnet ometer at a sampling r ate of 600 Hz, minimally filtered to (0.1– 200 Hz). Th e p ositions of the electrod es in add ition to fiduciary p oin ts, such as the n ose, n asion and preauricular p oints, w ere digitized u sing the 3Space Isotrak I I Sy s tem for subsequent precise co-registrat ion with MRI images. T he p osition of th e head with re- sp ect to the helium-co oled dew ar con taining th e measuremen t SQUIDS wa s determined by digitizing the p ositions of f ou r coils that w ere attac hed to the head. These four coils are subsequently emp lo y ed by the MEG sensors to assist in co-reg istration. Ongoing magnetic activit y wa s thus recorded un d er the follo wing three conditions: (1) at rest w ith ey es closed; (2) at rest with ey es op en; and (3 ) during sustained fin ger mo vemen t with ey es op en. As describ ed in previous studies, examination of raw data rev ealed strong alpha activit y wh en the sub j ect had ey es closed, and damp ened or no alph a rhyt hm when the sub ject had eyes op en (i.e., at rest or dur ing sustained fi n ger mov emen t conditions). The mu activit y w as strongest at rest (i.e., with ey es op en or closed) and suppr essed dur ing sus tained fingers mov emen t condition. 2 s long p erio ds of 18 C. J. LON G ET AL. ra w d ata with either large amplitude alpha wa v es (e.g., dur ing ey es-closed conditions) or large amplitude m u wa v es (e.g., during hand/finger resting conditions) w ere used as th e inpu t signal for MEG source lo calizatio n. In ad- dition, the empt y MEG r o om n oise w as recorded ju st b efore the s tart of the exp eriment and subsequently used as measurement noise in the estimation of the alph a and m u activit y generators. Figure 5 depicts th e net estimated current d istr ibutions c haracterized in eac h panel by the length of eac h cu rrent triplet. In the inv erse calculations, the orien tation of the estimated cur ren ts w as fixed such that they la y p er- p endicular to the cortical mesh . Eac h panel rev eals the relativ e effect of the state-space app roac hes as compared to those gained from th e static MNE tec hnique as the dynamics of the mu-rh ythm time course unfold (Figure 6 ). In all cases, we can observ e regions in th e cen tral sulcus sho win g strong activ ation, bu t the FIS and KF solutions sho w stronger activ ations than those obtained through the MNE metho d. Moreo v er, Figure 6 shows th at the reconstructed time cours es of the KF and FIS s olutions are muc h less v ariable than those of the MNE solution. In su mmary , analysis of the somatosensory exp eriment sho wed app ro xi- mate agreement b et ween MNE and the state-space mo d els. As in the simu- lated MEG exp eriment , w e n ote that the magnitudes of the estimate sour ces w ere greater for the stat e-space models (Figure 5 ). The temp oral dynamics of the FIS agreed m ore closely with the MNE than with the KF (Figure 5 ). Not surp risingly , b ecause th e state-space mo d els imp osed a temp oral con- strain t, their estimate s w ere smo other than the MNE estimates (Figure 6 ). Also, as in the simulated MEG exp erimen t, the KF estimates of the sou r ce amplitudes were larger than those computed b y eit her MNE or the FIS. 6. Discussion. W e ha v e dev elop ed and imp lemen ted a state-space mo del solution to the MEG in verse problem. W e examined its p erformance rel- ativ e to the MNE algorithm on simulated and actual MEG exp eriment al data. A simple analytic analysis sh ow ed that the state-space approac h of- fers tw o im p ortant impro vemen ts ov er MNE. Namely , the KF/FIS estimate equation ( 2.9 ) uses a dyn amic L 2 -norm regularizati on m ean equati on ( 2.6 ) and co v ariance m atrix equ ation ( 2.7 ) up d ate at eac h step, whereas the MNE estimate equation ( 2.5 ) uses a static L 2 -norm regularization zero mean and static co v ariance matrix to compute its instant aneous estimate. Making the prior mean zero at eac h up date b iases eve r y solution to wa r d zero. This mak es it difficu lt to id entify sources that are nonzero b ut w eak (low amplitude). In cont rast, the KF/FIS estimate is more likel y to iden tify a w eak source b ecause if the s ource estimate at the pr evious time p oin t, that is, 1 msec b e- fore, w as similarly we ak, it w ould b e com bined with the current observ ation to iden tify wh at would most likel y b e an atten uated source at the current time p oint . Second, the KF sour ce estimates u se all of the observ ations from DYNAMIC ST A TE-SP ACE METHODS F O R MEG 19 Fig. 5. Sour c e activation maps showing instantane ous r e c onstructions for the sensorimo- tor task. Each c olumn shows the temp or al evolution at the five p e ak m o des of the mu-rhythm for e ach of the thr e e m etho ds: M inimal Norm Estimator (MNE); Kalm an filter (KF); and Fixe d Interval Smo other (FIS). The r e d and blue c olor maps c orr esp ond to inwar d and out- war d curr ent flow, r esp e ctively. This r e c onstruction was c arrie d out using a r e gularization p ar ameter value of λ = 1 . 20 C. J. LON G ET AL. Fig. 6. R epr esent ative tim e c ourse r e c onstructions extr acte d fr om l o c ations in the so- matosensory c ortex at 492 ms in Figur e 5 : MNE (blue), Kalman filter (gr e en) and Fi xe d Interval Smo other (r e d). the start of the exp eriment to the cur ren t time, then th e FIS equation ( 2.13 ) uses the KF sour ce estimates to compu te new estimates based on all of the exp erimenta l observ ations. In con tr ast, the MNE source estimat e uses only the observ ations recorded at the current time. Giv en that the MEG up dates are computed ev ery m illisecond, the solutio ns at adjacen t time p erio ds are lik ely to b e strongly correlated. Our state-space approac h imp oses a d y- namic L 2 -norm regularization constraint to use this temp oral information, whereas MNE imp oses the same static L 2 -norm regularizatio n constrain t at eac h time p oin t and thus do es not consider the temp oral s tr ucture. O ur sim u lation stu dies d emonstrated th at the KF and FIS estimates are m ore accurate in terms of MSE (Figures 3 and 4 ) compared with MNE. Although the sp atial maps of the state-space and MNE h ad comparable spatial extent, the KF and FIS estimates captur ed m ore accurately the magnitude of the ac- tiv atio ns (Figure 1 ). This accoun ts for the low er MSE for the FIS and the KF relativ e to MNE despite small num b ers of erroneous activ ations in th e case of the state-space mo dels outsid e the areas of actual activ ations (Figure 1 ). In add ition, the tw o state-space s olutions were considerably smo other and in closer agreemen t with the simulat ed signal than the MNE solution. Th e KF whic h compu tes its estimates based exclusiv ely on data up to the current time had a temp oral lag with r esp ect to the tr u e signal (Figure 2 ). That is, the lo cations of the signal p eaks and troughs w ere offset with resp ect to their true temp oral lo cations. In addition, the KF estimates o v erestimated the true signal amplitud e (Figure 2 ). The backw ard pass of the FIS corrected b oth of th ese p roblems, y ieldin g a closer agreemen t with the true signal and DYNAMIC ST A TE-SP ACE METHODS F O R MEG 21 a lo wer MSE (Figure 2 ). Although in the analysis of the actual MEG exp eri- men t MNE and the state-space m o dels source estimates were in appro ximate agreemen t (Figure 5 ), the magnitudes of th e FIS estimates w ere more con- sisten t w ith those of the KF estimates, whereas th e temp oral dyn amics of the FIS estimates follo w ed more closely those of the MNE estimates (Fig- ure 5 ). The KF estimates of the sour ce amplitudes were larger than those computed by the either MNE or the FIS (Figure 6 ). Although our fi nd- ings establish the feasibilit y of using high-d im en sional state-space mo dels for solving the MEG in v er s e problem, they also suggest several exte nsions. In our analysis we compu ted the K F and FIS estimates b y testing differen t v alues of the regularizat ion paramete r. T h is parameter can b e estimated in an empirical Ba yesian [Lam us et al. ( 2007 )] or a fully Ba y esian f ramew ork. The s p atial comp on ent of the mo del can b e imp ro ved by u sing the structure of the particular exp erimental d esign to pr op osed sp ecific forms of the F state-transitio n matrix. There are a br oad r ange of tec hn iqu es that h a ve b een used to accele rate computations in high-d im en sional state-space m o d- els. The forward mo del can b e extended to include sub cortical sources. The MEG and EEG r ecordings are usually recorded simulta neously . Ou r state- space paradigm can b e applied to th e problem of estimating the sources f rom these tw o s ou r ces. These extensions will b e the topics of futur e rep orts. REFERENCES Black f ord, L. S. , Choi, J. , C lear y, A. , D’Azeve do, E. , Demme l, J. , Dhillon, I. , Dongarra, J. , H ammarling, S . , H enr y, G. , Petitet, A . , S t anley, K. , W al- ker, D. and Whaley, R. C . (1997). Sc aLAP A CK Users’ Guide . SIAM, Philadelphia, P A. Dale, A. M. , Fischl, B. an d Sere no, M. I. (1999). Cortical surface-based analysis. I. Segmentation and surface reconstruction. Neur oImage 9 179–19 4. Da unize au, J. and Friston, K. J. (2007). A mesostate-space mo del for EEG and MEG. Neur oImage 38 67–81. F arrell, B. F. and Ioannou, P. J. (2001). State estimation using a redu ced order Kalman filter. J. Atmosph eric Sci. 58 3666–368 0. Freeman, W. , A hlf ors, S. and Menon, V. (2009). Combining fMRI with EEG and MEG in order t o relate patterns of brain activity to cognition. Int. J. Psychophys iol. 73 43–52. Friston, K. J. , Harrison, L. , D aunizeau, J. , Ki ebel, S. , Phillips, C. , T r ujillo- Barreto, N. , Henson, R. , Flandin, G. and Ma ttout, J. (2008). Multiple sparse priors for the M/EEG inv erse problem. Neur oI mage 39 1104–11 20. Fukumori, I. and Malanotte-Ri zzoli, P. (1995). A n appro x imate Kalma n filter for ocean data assimilatio n: An example with an idealized Gulf stream mo del. J. Ge ophys. R es. 100 6777–679 3. Fukumori, I. , Be nvensite, J. , Wunsch, C. and Haid-vogel, B. D. (1993). Assimi- lation of sea surface topography into an o cean circulation mo del using a steady-state smoother. J. Phys. Oc e ano gr. 23 2162–21 81. Galka, A. , Y amashit a, O. Ozaki, T. et al. ( 2004). A solution to the dy namical inv erse problem of EEG generation using spatiotemp oral Kalman filtering. Neur oImage 23 435– 453. 22 C. J. LON G ET AL. H ¨ am ¨ al ¨ ainen, M. S. and Ilmoniemi, R. J. (1994). I nterpreting magnetic fields of the brain: Minimum norm estimates. Me d. Biol. Eng. Comput. 32 35–42. H ¨ am ¨ al ¨ ainen, M. S. and Sar v a s, J. (1989). Realistic conductivity geometry mod el of the human head for interpretation of neuromagnetic data. IEEE T r ans. Biome d. Eng. 36 165–171 . H ¨ am ¨ al ¨ ainen, M. S. , H ari, R. , Ilmoniemi, R. J. , Knuu tila, J. and Lounasmaa, O. V. (1993). Magnetoenceph alography—Theory , instrumentation, and app lication to nonin- v asive studies of th e wo rking human brain. Re v. Mo dern Phys. 65 413–497. Ka y, S. (1993). F undamentals of Statistic al Signal Pr o c essing, V ol . I—Estimation The ory . Prenti ce Hall, Upp er Sadd le River, NJ. Kit aga w a, G . and Gersch, W. (1996). Smo othness Priors Analysis of Time Series . L e ctur e Notes in Statistics 116 . Springer, N ew Y ork. MR1441074 Krim, H. and Vi berg, M. (1996). Tw o decades of arra y signal pro cessing research: The parametric app roac h . IEEE Signal Pr o c essing Magazine 13 67–94. Lamus, C. , Long, C. J. , Hamalaine n, M . S. , Bro wn, E. N. and Purdon, P. L. (2007). P arameter estimatio n and d ynamic source localization for th e magneto encephalogra- phy (MEG) inv erse problem. I EEE International Symp osium on Biome dic al Im aging, Washington, DC, April 12–15, 2007 . Liljestr ¨ om, M. , Ku jala, J. , Jensen , O. and Salme lin, R. (2005). Neuromagnetic localization of rhythmic activity in the human brain: A comparison of three metho ds. Neur oImage 25 734–745 . Lounasmaa, O. , H ¨ am ¨ al ¨ ainen, M. , Hari, R. and Salmelin, R. (1995). Information processing in the human brain: Magneto encephalographic approach. Pr o c. N atl. A c ad. Sci. USA 93 8809–8 815. Ma ttout, J. , Phillips, C. , Pen ny, W. D. , Rugg, M. D. and Friston, K. J. (2006). MEG source localization und er multiple constraints: An extended Bay esian framew ork. Neur oImage 30 753–768 . Mosher, J. C. and Leahy, R. M. ( 1998). R ecursive MUS IC: A framewo rk for EEG and MEG source localization. IEEE T r ans. Biome d. Eng. 45 1342–1354. Mosher, J. C. , Lewis, P. S. and Leahy , R. M. (1992). Multiple dip ole mo deling and localization from spatio-temporal MEG data. IEEE T r ans. Biome d. Eng. 39 541–55 7. Nunez, P. L. (1995). Ne o c ortic al Dynamics and Human EEG Rhythms . Ox ford Un iv. Press, New Y ork. Ou, W. , H ¨ am ¨ al ¨ ainen, M. S. and Golland, P. (2009). A d istributed spatio-temporal EEG/MEG invers e solv er. Neur oIm age 44 932–94 6. P ascual-Marq ui, R . D. , M ichel, C . M. and Lehm ann, D . (1994). Lo w resolution electromagnetic tomography: A new meth o d for lo calizing electrical activit y in the brain. Int. J. Psychophysiol. 18 49–65. Phillips, C. , Ma ttout, J. , Rugg, M. D. , Maquet, P. and Friston, K. J. (2005). An empirical Ba yes ian solution to th e source reconstruction problem in EEG. Neur oImage 24 997–1011. Salmelin, R. and Hari, R. (1994). Characterization of sp ontaneous MEG rhythms in healthy adu lts. Ele ctr o enc ephalo gr. Cl in. Neur ophysiol. 91 237–248. Salmelin, R. , Forss, N. , Knuutila, J. and Hari, R. (2005). Bilateral activ ation of the human somatomotor cortex by distal hand mov ements. Ele ctr o enc eph. Clin. Neu- r ophysiol. 95 444–452. Uutela, K. , H ¨ am ¨ al ¨ ainen, M. a nd Salmelin, R. (1998). Global optimization in th e localization of neuromagnetic sources. IEEE T r ans. Biome d. Eng. 45 716–723. Uutela K. H ¨ am ¨ al ¨ ainen, M. and Somersalo, S. (1999). Visualizatio n of magnetoen- cephalographic d ata using minim u m current estimates. Neur oIm age 10 173–18 0. DYNAMIC ST A TE-SP ACE METHODS F O R MEG 23 V an Veen , B. D. , v a n Dronglen, W . , Yuchtman, M. and Suz uki, A. (1997). Lo cal- ization of brain electrical activity via linearly constrained minimum v ariance spatial filtering. I EEE T r ans. Biome d. Eng. 44 867–880. Wipf, D. P. , O wen, J. P. , A ttias, H. T. , Sekiha ra, K. and Nagarajan, S. S. (2010). Robust Ba yesian estimation of the lo cation, orientation, and time course of multiple correlated neu ral sources using MEG. Neur oImage 49 641–655. Y amashit a, O. , Galka, A. , Ozaki, T. , Bisca y, R. and V aldes-Sosa, P. (2004). Recur- sive p enalized least squares solution for dy namical inv erse problems of EEG generation. Hum. Br ain Mapp. 21 221–235. Zumer, J. M. , A ttias, H. T. , Sekihara , K. and Nagarajan, S. S. (2007). A p roba- bilistic algorithm fo rm integrating source lo calization and noise suppression for MEG and EEG data. Neur oImage 37 102–115. Zumer, J. M. , A tti as, H. T. , Sekiha ra, K. and Nagarajan, S. S. (2008). Probabilistic algorithms for MEG/EEG source reconstruction using temp oral b asis funct ions learned from data. Neur oImage 41 924–940 . C. J. Long GlaxoSmithKline Clinical Imaging Center Hammersmith Hospit al London W12 0NN UK E-mail: cjlong1@gmail.com P. L. Purdon A thinoula A. Mar tinos Center for Biomedical Imaging Massachusetts General Hospit al Har v ard Medical School Charlestown, M assachusetts USA S. Temereanca A thinoula A. Mar tinos Center for Biomedical Imaging Massachusetts General Hospit al Har v ard Medical School Charlestown, M assachusetts USA N. U. Desai Dep ar tment of Electrical En gineering & Computer S cience Massachusetts Institute of Technology Cambridge, Massachusetts USA M. S. H ¨ am ¨ al ¨ ainen A thinoula A. Mar tinos Center for Biomedical Imaging Massachusetts General Hospit al Har v ard Medical School Charlestown, M assachusetts USA and Har v ard/MIT Division of Heal th Sciences & Technology Massachusetts Institute of Technology Cambridge, Massachusetts USA E. N. Brown Dep ar tment of Brain And Cognitive S ciences Massachusetts Institute of Technology Cambridge, Massachusetts USA and Dep ar tment of Anesthesia & Critical Care Massachusetts General Hospit al Boston, Massachusetts USA and Har v ard/MIT Division of Heal th S ciences & Technology Massachusetts Institute of Technology Cambridge, Massachusetts USA E-mail: en br o wn1@mit.edu