腦磁波訊號源估測與同調性造影

全文

(1)國立交通大學資訊工程學系碩士論文. 腦磁波訊號源估測與同調性造影. Neuromagnetic Source Estimation and Coherence Mapping of Brain Activities. 研究生：鄭志瑜指導教授：陳永昇. 博士. 中華民國九十四年八月.

(2) Neuromagnetic Source Estimation and Coherence Mapping of Brain Activities. A dissertation presented by. Chih-Yu Cheng. to. Department of Computer Science and Information Engineering. in partial fulfillment of the requirements for the degree of Master of Science in the subject of. Computer Science and Information Engineering. National Chiao Tung University Hsinchu, Taiwan 2005.

(3) Neuromagnetic Source Estimation and Coherence Mapping of Brain Activities. c 2005 Copyright by Chih-Yu Cheng.

(4) 摘要. 腦磁圖儀非侵入地量測腦部活化源誘發出的磁場訊號。其能夠在高時間解析度的取樣下對腦部活化源做時空造影，利用此方式將可有助於人腦功能在臨床與基礎神經生理學方面的研究。在此論文之中，我們將提出新的空間濾波技術，其利用腦磁圖儀所量測的訊號記錄來估測腦中訊號源的顯著程度以及訊號源間的相關性。利用磁圖儀所量測的紀錄來估算腦部活化源的問題稱之為逆估算問題。為了解此類的不適定問題，加入些假設與限制是需要的，例如等效電偶極的腦部訊號源模型、固定訊號源電偶極個數等假設以及結構性或是最小範數限制等。在眾多的訊號源逆估算方法中，光束構成法—一種空間濾波技術—在近十年來漸漸的引人注目。此訊號源逆估算法是藉由探測訊號源空間中一個接一個的體素，此時空間濾波器在每個位置會分別地被計算。此空間濾波器依據單位增益與最小變異量兩準則將可保留目標訊號源訊號並同時抑制其他訊號所造成的影響。然而，如何決定訊號源電偶極的方向會是個問題。在既有的文獻中主要提及三種方式，第一、此電偶極方向會垂直於該處大腦皮質表面，然而自動地且精確地重建出曲折的皮質表面是非常困難的；第二、利用搜尋的方式找出，但此方法是非常耗時的；第三、將訊號源電偶極拆解成三個互成正交的分量，但其有遺漏偵測的風險。此論文中，我們開發了一個新空間濾波技術，稱之為最大對比光束構成法，用以對神經訊號源做統計量的造影。除了如同傳統的光束構成法一樣的採用單位增益與最小變異量準則，我們所提的方法還利用了最大對比（活化與休息狀態逆估算出的神經活化訊號變異量的比值）準則。藉由此最大對比準則，我們將可用閉形式解的方式解析地決定出訊號源電偶極方向，也就是說，對於某一位置而言，能夠有效率地得出其對應的空間濾波器。一旦藉由對腦磁圖儀所量測訊號濾波所估算出各個位置的訊號源，將可計算對應整個訊號源空間的 F 統計量以視覺化在活化狀態時於皮質區的訊號源相對應於休息狀態的顯著程度。.

(5) 同調性訊號源造影是另一個在人腦功能研究上有趣的主題，其能有助於探索人腦相關功能區域間聯繫的機制。近來，在特殊頻帶下的同步震盪訊號被認為與神經網路間的溝通有著密切的關係。同調訊號源造影法利用光束構成法為基礎，對相對於一參考訊號擁有同調性的皮質區訊號源做造影。但其僅適用於同步頻帶具有穩恆的特性。在此，我們提出另一新技術—最大正規化相關光束構成法—利用最大正規化相關性的準則解析地且閉形式的決定訊號源方向以對同調訊號源造影。理論上，此方法是一般化的，也就是它能對在任一有定義自相關性與互相關性的定義域中互有相關性的訊號作造影。此論文中，我們在感興趣的 Morlet 小波域中計算腦磁波儀所量測訊號的自相關性以及其與參考訊號的互相關性，接著利用所提出的最大正規化相關光束構成法來對動態跨頻同調訊號源造影。模擬、假體以及手指抬動的實驗被採用來驗證與評估所提方法的正確性與能力。根據模擬與假體實驗所得的分析結果，我們的方法的確，第一、能既有效率地並準確地決定出訊號源電偶極的方向；第二、能準確地分別定位出擁有顯著變異量與時頻同調性的訊號源。當應用在手指抬動的實驗分析時，我們可從 F 統計量圖明顯的指出在腦部感覺運動區在活化狀態相對於休息狀態有高度的對比。.

(6) Abstract Magnetoencephalography (MEG) non-invasively measures the electromagnetic signals induced by brain activities. It can provide spatiotemporal brain activation imaging with high temporal resolution to facilitate functional brain research in both clinical and basic neuroscience fields. In this thesis, we propose novel spatial filtering techniques for statistical mapping of neuronal sources as well as cortical oscillatory coupling by using the wholehead MEG recordings. The problem of estimating the activation sources in the brain from the MEG recordings is called the inverse problem. To solve this ill-posed problem, approximations such as equivalent current dipole for source modeling, assumptions such as a fixed number of dipoles during the task, and constraints such as anatomical constraint and minimum-norm constraint are required to limit the solution space. Among the various kinds of source estimation methods, beamforming technique, a kind of spatial filtering technique, has becoming more and more attractive during the past decade. By probing the source space in a voxel-by-voxel manner, a spatial filter is individually calculated for each position. This spatial filter can reconstruct the activation magnitude of the targeted source while suppressing the contribution from other sources by imposing the unit-gain constraint and by applying the minimum-variance criterion. However, the determination of dipole orientation can be problematic. There are three major kinds of methods proposed in the literature. First, the dipole orientation is aligned to be perpendicular to the cortical surface. Unfortunately, the surface reconstruction for the convoluted cortex is very difficult and the reconstruction deviation will decrease the accuracy of the orientation. Second, the dipole orientation is determined by (exhaustive) search, which is time-consuming. The third kinds of methods decompose the dipole into three orthogonal components, which may suffer the risk of miss-detection. In this work, we develop a novel spatial filtering technique, called the maximumcontrast beamformer, for statistical mapping of neuronal sources. In addition to the unitgain constraint and the minimum-variance criterion, as in the conventional beamformers, our method exploits a maximum-contrast criterion that can maximize the discrimination between the estimated neuronal activities in the active state and those in the control (or i.

(7) resting) state. The maximum-contrast criterion helps to analytically determine the dipole orientation in a closed-form manner and the spatial filter can be obtained very efficiently for each targeted position. Once the neuronal activity waveform is estimated in the source space by spatially filtering the MEG recordings, F-statistic map can be calculated to reveal cortical regions with significant difference of activities between the control and active states. Another interesting issue in functional brain studies is the coherent source mapping for probing the binding mechanism of connected functional assemblies. Recently, oscillatory synchronization in particular frequency bands has been shown to be closely related to the communication within a neural circuit. By using the beamforming-based algorithm, DICS (Dynamic Imaging of Coherence Source) method can map the cortical sources that are statistically coherent to a specified reference at a certain frequency band. The limitation of the DICS method is that the synchronization frequency band is considered to be stationary during the task. Here, we propose another new method, maximum-normalizedcorrelation beamformer, for the mapping of the cortical oscillatory coupling. Theoretically, this method is very general that can image the sources correlated in the domain where the autocorrelation and crosscorrelation can be defined. To demonstrate the capability of this method, we compute the autocorrelation and crosscorrelation for the MEG recordings in the Morlet wavelet domain and image the dynamic coherent sources across multiple frequency bands during the task. Moreover, the dipole orientation has a closed-form solution by applying the maximum-normalized-correlation criterion. Experiments with simulation, phantom, and real data are conducted to verify the correctness and to assess the capability of the proposed methods. According to the experiments with simulation and phantom data, our methods indeed can efficiently and accurately calculate the dipole orientation. Also, our methods correctly locate the sources with significant variance and significant time-frequency coherence. When applied to a finger-lifting study, F-statistic map computed from the estimated neuronal activities on the cortical surface clearly identify the sensorimotor area with high contrast.. ii.

(8) 誌謝. 本篇論文得以完成，首先要誠摯的感謝自小到大花費無盡心血栽培我的父母親。有他們在背後的支持，我才得以完成此碩士學位。再則，感謝我的指導教授陳永昇老師以及陳麗芬老師，在這兩年多來引領我如何做研究，並對於我的論文提出了許多寶貴看法與意見。他們所給予的助益，是本篇論文完成的直接重要因素。最後，感謝實驗室的同儕們，在這段日子以來，對我直間或間接的幫助。. iii.

(9) iv.

(10) Contents List of Figures. vii. List of Tables. ix. 1. Introduction 1 1.1 Backgrounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Thesis Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10. 2. MEG Forward Prediction with Spherical Head Model 13 2.1 General Forward Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Source Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Forward Solution for Spherical Head Model . . . . . . . . . . . . . . . . . 17. 3. Beamforming Methods for Functional Brain Imaging 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Source Estimation . . . . . . . . . . . . . . . . . . . . . 3.2.1 Data Model . . . . . . . . . . . . . . . . . . . . 3.2.2 Spatial Filter Design . . . . . . . . . . . . . . . 3.2.3 Source Space Scanning Approaches . . . . . . . 3.2.4 Maximum Contrast Beamformer . . . . . . . . . 3.3 Coherence Mapping . . . . . . . . . . . . . . . . . . . . 3.3.1 Signal Similarity . . . . . . . . . . . . . . . . . 3.3.2 DICS . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Maximum Normalized Correlation Beamformer . 3.3.4 Time-frequency Coherent Sources Mapping . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. 23 24 27 27 30 33 37 43 43 44 46 47. Experiment Results 4.1 Simulation and Phantom Experiment . . . . . . . . . . . 4.1.1 Forward Model Accuracy . . . . . . . . . . . . 4.1.2 Regularization Term β . . . . . . . . . . . . . . 4.1.3 Imaging of Source Power and Coherent Sources .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. 53 57 57 58 60. 4. v.

(11) 4.1.4 4.1.5. 4.2 5. Estimation of Source Orientation by Maximum Contrast Beamfomer Estimation of Source Orientation by Maximum Normalized Correlation Beamformer . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.6 Computational Cost of the Proposed Algorithms . . . . . . . . . . Experiment of Self-paced Finger Movement . . . . . . . . . . . . . . . . .. Conclusions. 65 66 67 72 79. Bibliography. 83. vi.

(12) List of Figures 1.1 1.2 1.3 1.4. A finger lifting example about brain functional mapping. . . . . . . . . . . Large scale integration problem. . . . . . . . . . . . . . . . . . . . . . . . Acquisition of the electric and induced magnetic signal by EEG and MEG. . Forward and inverse problem and their relationship. . . . . . . . . . . . . .. 2.1 2.2. MEG forward model with spherical head model. . . . . . . . . . . . . . . . 19 MEG forward example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21. 3.1 3.2. Basic concept of spatial filter . . . . . . . . . . . . . . . . . . . . . . . . Example of the application of the beamforming technique: radar array system and MEG inverse problem. . . . . . . . . . . . . . . . . . . . . . . . Spatial response of the unit-gain constraint and minimum variance spatial filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Determination of the position and orientation probed by beamformer when given the geometry model of the cortical surface. . . . . . . . . . . . . . Effect of miss-detection . . . . . . . . . . . . . . . . . . . . . . . . . . . Determination of the position and orientation probed by beamformer in general case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Examples of choosing the interesting region on time-frequency map . . . Procedure of imaging (event-related) time-frequency . . . . . . . . . . .. . 24. Concept of signal space projection (SSP). . . . . . . . . . . . . . . . . . Usual preprocessing for MEG recordings. . . . . . . . . . . . . . . . . . Configuration of phantom. . . . . . . . . . . . . . . . . . . . . . . . . . Dipole fitting error for checking the accuracy of Sarvas forward model by phantom experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Computer generated sources for simulation. . . . . . . . . . . . . . . . . F-statistic map of the estimated sources by maximum contrast beamformer with the simulated recordings. . . . . . . . . . . . . . . . . . . . . . . . Time-frequency coherence map of the estimated sources by maximum normalized correlation beamformer with the simulated recordings. . . . . . .. . 55 . 56 . 58. 3.3 3.4 3.5 3.6 3.7 3.8 4.1 4.2 4.3 4.4 4.5 4.6 4.7. vii. 3 5 6 7. . 26 . 31 . 35 . 36 . 37 . 49 . 51. . 59 . 62 . 63 . 64.

(13) 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15. 4.16. Error of maximum contrast beamformer to determine the orientation of tangential source with simulated recordings. . . . . . . . . . . . . . . . . Error of maximum contrast beamformer to determine the orientation of arbitrary source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Error of maximum contrast beamformer to determine the orientation of arbitrary source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The effect of the radial component of dipole source upon estimated dipole orientation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Error of the maximum normalized correlation beamformer to determine the orientation of arbitrary source. . . . . . . . . . . . . . . . . . . . . . . . F-statistic map of the estimated sources by maximum contrast beamformer with the recordings from left index finger movements experiment. . . . . F-statistic map of the estimated sources by maximum contrast beamformer with the recordings from right index finger movements experiment. . . . . Time-frequency coherence map of the estimated sources by maximum normalized correlation beamformer with the recordings from left index finger movements experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . Time-frequency coherence map of the estimated sources by maximum normalized correlation beamformer with the recordings from right index finger movements experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . .. viii. . 67 . 68 . 69 . 70 . 71 . 75 . 76. . 77. . 78.

(14) List of Tables 4.1 4.2. The effect of the regularization term β . . . . . . . . . . . . . . . . . . . . 61 Comparison of SAM and proposed method in computing time . . . . . . . 73. 5.1. Comparison of SAM, LCMV and our methods. . . . . . . . . . . . . . . . 81. ix.

(15) x.

(16) Chapter 1 Introduction.

(17) 2. Introduction. 1.1. Backgrounds. In the past three centuries, the development of technology is more and more flourishing and the world rapidly change in this era. For instance, light has replaced candle to light in night, car has replaced horse-drawn vehicle to transport goods, abacus has been replaced by electric calculator to do accounts and even the computer equipping artificial intelligence could beat all humans on chess table. Almost every thing around our life is greatly effected by this technology progressing. Human knows the principles of nature to create those mentioned above. From the theory of mechanism, heats to electromagnetism, human beings create textile, steam machine and motor. Medical science brings the better prevention of and treatment of disease to extend humans’ lift. Agriculture science brings more quantity of output to supply the increasing population. Aeronautical science make the dream of fling come true by the appearance of airplane. Many people donate on each field to think and develop new methods or tools to push the era forward. All ideas of those principles and invents come from a powerful device, human brain. Human brain is often analogized with the most powerful computer. That brain controls hands, feet, eyes and most organs is like that the control unit in the central processing unit (CPU) coordinates memory, peripheral devices in the whole computer system. Brain processes the information receiving by eyes, ears, nose resembles that the arithmetic/logic unit (ALU) in CPU processes the information which inputs through keyboard and mouse. Also, brain supports parallel processing so that we can read the news on web, press the mouse button, drag it to scroll the page, listen music and follow the melody to hum the tune simultaneously. Besides, brain wins the computer in several aspects. The most advantage of human brain is that it owns the originality being able to cerate infinite possible, whereas the computer is just a deterministic machine. Moreover, it could perform so many functions but its size is magically small. ”Deep blue” that is match for the best human player on chess weights 1.4 ton about 1000 times of the matured human brain, whereas it is just able to chess. Human brain is what a dedicate device! Of course that there is no reason that no one has interest to know the mechanism of the operations in this mystical and attractive device. From many centuries ago, the studies.

(18) 1.1 Backgrounds. Figure 1.1: A finger lifting example about brain functional mapping. When we perform a certain functional task, the corresponding cortex areas in the brain will activate. From sagittal, coronal and transverse view, the figure shows significance of the average power of the estimated signal by beamforming-based method using MEG recordings. It evidently shows that the sensorimotor area activates during the period of finger lifting.. on brain started form the anatomy to explore the structure of it. Then, in the later of nineteen century, the electric brain signal was found. Through the researches in early and mid twenty century, it was more and more clear that when we perform a certain task, the neurons at the corresponding cortex area in the brain will activate. The activation is the electrical signals resulting from the irons flowing in and out the cells of the neurons to perform the information interchange between the neurons. For instance, when we look at something, the neurons in the visual cortex area, at the occipital part of the brain will activate. Figure 1.1 gives another example that sensory and motor areas activate when we lift our finger up. However, comparing with other technologies, the functionality of the human brain is still a virgin land to scientists. We still have few idea about our brain yet. It is very difficult to explore the human brain because of the extraordinary sophistication and complexity of it, which consists of billions of micron-sized neurons connecting with each other to knit a very complex neural network for performing many works. Besides, due to morality, there are many limitations of the researches on living bodies. For instance, it can not force an ordinary one to be incised his scalp to denoting his brain for studies. Hence, the progress. 3.

(19) 4. Introduction. of the studies on the field of functionality of the human brain is comparatively slow. Today, by the progress of the technologies in both hardware and software, there are many useful instrumentations can help us to explore the functionality of human brain more accurate and more convenient, even some can under the non-invasive condition, that is, during the experiment, there is no physical harm to the subject. Now, for non-invasively studying on brain functionality, there are three primary modalities, which are functional magnetic resonance imaging (fMRI), electronic encephalograph (EEG) and magnetic encephalograph (MEG). All of these three receives the brain signal outside the head by sensors and can be think that they take the picture for the phenomena caused by the brain electrical activities. fMRI is a widely-used device measuring the changes in blood oxygenation which is related to the neural activities. It has the advantage of high spatial resolution, but with the disadvantage of low temporal resolution. EEG and MEG directly measure the electrical potential and induced magnetic filed caused by the iron flow respectively. Besides, they own the ability to measure the signal with high temporal resolution up to the order of millisecond. Temporal resolution of the instructions for the studies of brain functionality is important. Because, the the brain activities could dynamically change up to about 70 ∼ 80 Hz, that is, by the Nyquist sampling theorem [50], the sampling frequency at least needs about 160 Hz to reconstruct the origin signals. In addition, recently, more and more studies direct that the large-scale integration in the brain is done by the information involving in the frequency domain of the signals [3, 4, 56, 64, 67, 71]. For a functional task, there are not only the primary but other regions associating with it. It is the ensemble functionality of the involved neural network that achieve this task. For example, when we lift our finger up, the regions corresponding visual, auditory, dorsolateral prefrontal association, somatosensory and motor cortex may be involved. How those spatially-distributed regions know that they cooperate in the same job is so called large-scale integration problem, or say binding problem. Although there is still no concrete answers to the problem, the most plausible candidate is by the modulation of the wave with certain frequency bands. The concept is illustrated in the Figure 1.2. For these reasons, EEG and MEG plays indispensable roles for the study of brain connectivity because they are with the ability of high temporal resolution. Figure 1.3 illustrate the acquisition of EEG and MEG recordings..

(20) 1.1 Backgrounds. Figure 1.2: Large scale integration problem. When performing a certain task, several cortex areas (circle) may be involved. Cooperation mechanism of these involved cortex areas is still in debate. However, more and more studies consider that the most possible candidate is the rhythm which modulates the signals generated by the cooperated cortex ares. Thus, those signals are correlated in the certain frequency band(s) and are called as coherent sources.. Nevertheless, not only low signal to noise ratio (SNR) but also ill-posed inverse problem are two major difficulties in the studies of brain functionality by the modality of EEG or MEG. First, the scale of the brain electrical signals is very small compared to the environmental noise and confoundings. For example, the magnitude of the induced magnetic field is around 10−13 ∼ 10−15 Tesla while the omnipresent earth’s static filed is about 10−5 Tesla (Tesla is the unit of magnetic field). Hence, both a proper environment for experiment and the signal preprocessing before analyzing or further processing are necessary. Second, the recordings on EEG/MEG sensors are not induced by a single source but are the superpositions of the contributions induced by sources distributed all over the brain. It can be considered as that we set some microphones for recording the voice in a clamorous cocktail party and hence each microphone will record the mixture of the voice from all the people in the party rather than an individual one. Although there are many researches based on EEG or MEG modalities directly analyze the recordings, they may have bias due to the mixed effect. The problem of estimating the individual signal of each source from the EEG or MEG recordings measured outside the scalp is called the inverse problem. Shown by Helmholtz (1853), the difficulty in solving the problem is that it mathematically is an ill-posed problem which has no unique solution but infinite. That is, we can not solve it easily.. 5.

(21) 6. Introduction. (www.neuromag.com). Figure 1.3: Acquisition of the electric and induced magnetic signals by EEG and MEG. The upper and lower parts respectively are for EEG and MEG. Both of they acquire the time course recordings by several sensors locating outside the scalp with the sampling rate up to the order of kilo-Hz. Due to such high temporal resolution property, they are adaptive to be the modalities for the studies on brain functionality. Before describing the individual algorithms for solving the inverse problem, it needs to know that the concepts of most algorithms involve the Bayesian statistic framework: P (x|y) =. P (x)P (y|x) . P (y). (1.1). P (x|y) denotes the conditional probability of x given y while P (x) and P (y) are the probability of x and y respectively. P (y|x) denotes the conditional probability of y given x. Applying to MEG inverse problem, let x be the parameters of the sources in the head and y be the recordings on the sensors such that P (x|y) describes the inverse problem. The Bayesian equation shows that if we want to know the source parameters by the recordings, we can transform the problem to the form at the right side of the equal sign. Thus, knowing P (y|x) is a necessary condition. P (y|x) describes the probability of the recordings when given the source parameters, that is, the forward problem, which has the unique solution when we know the electric source and head (conductor) model. The problem and its solution will be briefly describes in Chapter 2. Figure 1.4 illustrates forward, inverse problem.

(22) 1.1 Backgrounds. 7. Inverse P roblem. F orward S olution. F orward P roblem Figure 1.4: Forward and inverse problem and their relationship. Given a putative source in the brain, what the magnetic field distributes outside the head is the forward problems (right to left). Oppositely, the inverse problem describes how we reconstruct the brain sources by the recordings on sensors outside the head (left to right). To solving the inverse problem is always involving the solution of the forward problem which has the unique solution when we knows the electric source (light green arrow) and head (conductor) model.. and their relationship. As aforementioned, the inverse problem has infinite solutions; hence, several algorithms based on various constraints were been developed in the literatures to derive a unique solution for the problem. These algorithms can be roughly divided into primary two categories that think a source as a focal current or distributed currents respectively. These two assumptions are both reasonable, because we can find a group of distributed currents that effect the same contribution to the EEG/MEG recordings with a specific focal current, and vice versa. The representative algorithms in first category is based on dipole modeling and spatial filter. The dipole modeling method is a nonlinear searching method to determine the parameters of the given number sources, the least square dipole fitting algorithm [27, 30, 49, 58] is an example. The spatial filter method estimates the source by linear combining the recordings measured by sensors. The coefficients for the linear combination are adaptively determined by different criteria in different algorithms. Beamforming [22, 54, 59, 72] and multiple signal classification (MUSIC) [46, 48] algorithms are the examples. The other category assumes the sources is the distributed currents and the minimum norm estimation (MNE) [26, 43] method belongs to it..

(23) 8. Introduction. Solving the inverse problem and visualizing the results can be think as the imaging process for the brain signals. That is, after applying the inverse procedure on the recordings, we can ”see” the (estimated) activities in the brain. For the studies of the basic functional brain mapping, imaging where the source current magnitude (power) is significant gives the direct information reflecting the activation of the neurons there. Most inverse algorithms have proposed in the literatures belongs to these category. For another interesting and mystical problem, the binding problem, studies often exploit a measurement, ”coherence”, to quantify the connective level between different neural groups distributed in the brain. It results from that, as aforementioned, several studies reveals that the cooperative neural groups may share the modulation in certain frequency bands, such that the statistic in these bands should be similar between those groups and therefore coherence which is a measurement to judge the similarity of the signals in frequency domain is a proper media to implicitly point out the significance of the connective likelihood of those neural groups. Recently, some inverse algorithms, like [13, 20, 22, 31], were developed to image the coherent level of the (estimated) brain sources under a narrow frequency band. Besides, Sun et al. [66] and Sekihara et al. [61, 63] respectively presented methods based on least square dipole fitting approach and MUSIC involving the time-frequency representation to estimate the brain sources which best match the specified the time-frequency behaviors.. 1.2. Thesis Scope. In this thesis, we focus on beamforming-based algorithm for solving inverse problem with MEG modality. The advantage of beamforming-based algorithm will be described at the first section of Chapter 3. Besides, though most algorithms or procedures associating MEG are still suitable for EEG, we choose MEG be the modality here due to first, the SNR of the recordings usually is higher than EEG and, second, the forward solution can be derived by a simpler way with sufficient accuracy simultaneously. Beamforming-based inverse algorithm, or say beamformer, is the spatial filter, which can be think as putting a virtual sensor on a specific position with a specific direction to probe the source there according to the criteria of minimizing total variance of the esti-.

(24) 1.2 Thesis Scope. mated source and retaining unit gain to the specific position and direction. Theoretically, everywhere and every orientation in the brain is need to be probed for imaging the activities distributing over the whole brain. However, it is impossible due to the computational limitation. For the end of overcoming the problem, sampling the position space is adopted and the sampling rate is dependent on the desired spatial resolution. However, at one position, it just needs to determine a direction matching the orientation of the source there to be probed. Linearly constrained minimum variance spatial filter proposed by Van Veen et al. [72] samples the orientation solution space with three probing directions which are orthogonal with each other. Consequently, it sums of the results probed on the three directions and views the summation as the estimated result. However, because the intrinsic property, unit gain constraint and minimum variance, of the beamforming-based algorithm, sampling will flat the source distribution, even miss the detection of the source in the extreme condition. Robinson and Vrba proposed a method, called ”synthetic aperture magnetometery” (SAM), [54] which searches the optimal orientation according to the maximum z-statistic criterion. They also show the resultant distribution of the estimated brain signal is more focal by SAM than LCMV especially when the SNR is high [73]. The major drawback of SAM is that the computational cost is heavy. Besides, almost studies about sources coherence for the binding problem analyzed the MEG/EEG recordings directly, like [1,32,35,38,39,57], but Gross et al. [22] proposed ”dynamic imaging of coherent sources (DICS)” method which is yet the only beamformingbased inverse algorithm, to image the coherent sources. It performs in a narrow frequency band when given a reference signal. The method determines the source orientation by find the principle axis (maximum correlation) in this orientation solution space. For the sake of overcoming the tradeoff between computational cost and resultant accuracy in the issue of determining the source orientation , in this thesis, we propose two new beamforming-based algorithms for imaging the source significance and coherent sources respectively. Also, we apply the later method to image the time-frequency coherent sources for an extension of DICS. The term, time-frequency coherence, describes that the coherent frequency of two sources are not limited in a narrow band but may cross two or several bands with time, that is, in time-frequency domain. In the following, we briefly describe the two beamforming-based algorithms and the imaging method of the time-frequency co-. 9.

(25) 10. Introduction. herent sources proposed in this thesis. I. Maximum contrast beamformer Maximum contrast beamformer is a beamformingbased algorithm to image the source significance by MEG/EEG recordings with the properties both of efficiency and no probability of miss-detection. It still follows the minimum variance and unit gain constraint criteria as conventional beamformingbased algorithms to solve the MEG inverse problem, but the source orientation is determined by maximum contrast criterion with a deterministic computational steps. II. Maximum normalized correlation beamformer Maximum normalized correlation beamformer is a beamforming-based algorithm to image the correlated sources with a given reference signal by MEG/EEG recordings with the properties both of efficiency and no probability to miss-detection. It still follows the minimum variance and unit gain constraint criteria as conventional beamforming-based algorithms to solve the MEG inverse problem, but the source orientation is determined by maximum normalized correlation criterion with a deterministic computational steps. Theoretically, it is able to image the correlated sources not only in the frequency domain but also in any domain that the auto- and cross-correlation have been defined. III. Imaging the time-frequency coherence We apply the proposed beamforming methods to image the coherent sources cross different frequency bands with time. by computing the auto- and cross-correlation of the selected wavelet coefficients both of one (estimated) targeted signal and one reference signal. In this work, we implemented the proposed methods in C/C++ language on Win32 platform. The simulation and phantom experiments are conducted to input theirs recordings into the proposed algorithm for verification and demonstration. We also apply our methods to the recordings recorded from real finger lifting experiments.. 1.3. Thesis Organization. After this chapter, we start to introduce the forward problem and its solution for the preparation for the inverse algorithms described later. Consequently, in Chapter 3, af-.

(26) 1.3 Thesis Organization. 11. ter introducing some typical inverse algorithms first and mentioning the advantage of the beamforming-base method, we elaborate the basic of beamformer for imaging the source significance. Next, the aforementioned drawbacks of LCMV and SAM are pointed out more detail respectively before we propose the maximum contrast beamforming method with the detail deduction. Then, the principles of the DICS and maximum normalized beamforming method which image the coherent sources are expounded in order. At the end of that chapter, we describe the procedure for imaging the time-frequency coherent sources. Chapter 4 consists of the results of the experiment for the verification and demonstration of the proposed methods. Finally, a short conclusion for the thesis is given at the last chapter..

(27) 12. Introduction.

(28) Chapter 2 MEG Forward Prediction with Spherical Head Model.

(29) 14. MEG Forward Prediction with Spherical Head Model. According to the Bayesian statistic framework mentioned in last chapter, the inverse algorithms have to involve the forward solution to estimate the parameters of the brain sources when given the recordings measured by MEG sensors. Therefore, we need to, at first, construct the MEG forward model which describes the distribution of the the magnetic field outside the head when giving putative current sources in the brain. The general forward model can be deduced from quasi-static Maxwell equation and Biot-Sarvas law; however, in one hand, to analyze the brain behavior to the neuron-sized level is very difficulty; in the other hand, there is no analytical solution for the general form forward model, but needs to be computed numerically. In the MEG-related studies, it usually uses equivalent current dipole (ECD) to approximate the total electromagnetic effect caused by a group of neighboring neurons such that we can microscopically view the neurons’ activities for simplification. Besides, the spherical head model is often adapted for deriving an analytical solution of the forward model to reduce the computational cost. In the following, we will briefly introduce the forward problem from the general model, ECD source model, the spherical head forward model to, finally, the formulation of deriving the MEG recordings when giving sources.. 2.1. General Forward Model. The frequency band of almost interesting bio-electromagnetic signals are under 100Hz. At this condition, the physics of MEG can be described by quasi-static Maxwell equation [27]. Under this assumption, when a current density η(rη ) locates at rη in a closed volume G, the magnetic field b(r) at location r outside G can be described by the well-known Biot-Sarvas law r − rη µ0 Z η(rη ) × drη , b(r) = 4π G kr − rη k3. (2.1). where µ0 is the permittivity of free space. The overall current density can be view as the sum of two components, passive and primary parts. In our case, the primary currents ηp (rη ) associates the original neural activity, and the passive part, or say volume current ηv (rη ), is the result of the the induced potential of the conducting medium. The relation can be.

(30) 2.1 General Forward Model. 15. formalized as η(rη ) = ηp (rη ) + ηv (rη ) = ηp (rη ) + σ(rη )e(rη ) = ηp (rη ) − σ(rη )∇v(rη ),. (2.2). where e(rη ) is the electric field that equals the negative gradient of the electric potetial v(rη ), and σ(rη ) is the conductivity of the mediums which indicates the brain tissues in our case. Typically, the head can be think include three to five regions, for example, scalp, skull, cerebrospinal fluid, gray matter and white matter [47]. Moreover, the conductivity of each are considered constant and isotropic such that the Biot-Sarvas law can be rewritten as the sum contributed by the primary and volume currents: Z r − rη µ0 X (σi − σj ) dS0ij , v(rη ) × b(r) = b0 (r) + 3 4π ij kr − rη k Sij. (2.3). where σi and σj are the conductivities of the successive mediums. The second term at the right of the equal sign is the magnetic contribution from the volume currents formed by summing the integrating over the cross-medium surface S0ij , eg. scalp-skull, gray-white matter, and so on. The first term is the primary magnetic field: r − rη µ0 Z ηp (rη ) × drη (2.4) b0 (r) = 4π G kr − rη k3 which is the magnetic field contributed by primary source only and the corresponding primary potential field is 1 Z r − rη ηp (rη ) · drη (2.5) 4πσ0 G kr − rη k3 We need to know that the potential on all surfaces of mediums for deriving a general v0 (r) =. equation describing the magnetic field. By Green’s theorem, the relation of the potential at r and the potential of the source can be formulated [58] [47] as Z 1 X r − rη (σi + σj )v(r) = 2σ0 v0 (r) − (σi − σj ) v(rη ) × dS0ij . 3 2π ij kr − rη k Sij. (2.6). Up to now, if the parameters of head is known, we can substitute Eq. (2.4) and Eq. (2.6) into Eq. (2.3) to solve the forward problem, that is, when given a source inside the head, we can derive the magnetic field anywhere outside the scalp, under the assumption of quasistatic magnetic field computation and that of constant isotropic conductivity of the homogeneous regions in the brain. However, the solution, in general, is not analytical and need to be solved numerically..

(31) 16. MEG Forward Prediction with Spherical Head Model. 2.2. Source Model. According to the knowledge of neural-physiology, the active large pyramidal neurons locating at gray matter gives the primary contribution to the MEG recordings owing to two primary reasons. First, these neurons arranged at a regular direction which is perpendicular to the cortical surface so the contribution from each neuron will not be canceled out like other disorderly-arranged neurons. Second, the duration of the electrical signal at this place is more longer than others [2] such that the integration of their power is larger enough to be measured. Although the actual electric activities take place within each neuron, we always model the activities from a macro view. It is owing to the fact that, first, a large number of neighbor neurons activate together and that, second, it is very difficult to study the neuron activities finely to the neural level with the device like MEG which equips rather few number of sensors compared with the number of the neurons. Hence, for convenience and simplicity, a small patch of activated cortex centralized at r0 is considered and the primary current density can be approximately represented by the current dipole representation which is called equivalent current dipole (ECD) with the parameters of position, orientation and magnitude. Thus, we can view the large numbers of activated neurons as a point of dipole source: ηp (rη ) = qδ(rη − r0 ),. (2.7). where δ is the Dirac delta function and the q is the dipole moment defined by q=. Z. ηp (rη ) drη .. (2.8). With the dipole moment representation, Eq. (2.4) and Eq. (2.5) can be respectively simplified as b0 (r) =. r − r0 µ0 q× , 4π kr − r0 k3. (2.9). v0 (r) =. 1 r − r0 q· . 4πσ0 kr − r0 k3. (2.10).

(32) 2.3 Forward Solution for Spherical Head Model. 2.3. 17. Forward Solution for Spherical Head Model. For the sake of computing the solution for the forward problem formulated as Eq. (2.3), the information of the head is needed in advance. Realistic head model can be derived by segmenting the surface of each tissue on the volume image scanned by X-ray computed tomography (CT) imaging or magnetic resonance (MR) imaging instruments. Then, the forward solution can be calculated by finite element method (FEM) or boundary element method (BEM) [15]. However, calculating the solution by both of these methods are timeconsuming because they have to solve the problem numerically. A spherical head model is a simple model of the head with the advantage of fast computing for forward problem and also retains an acceptable accuracy to the MEG applications [2, 47]. The model assumes the head be comprised by a set of nested concentric spheres. At this condition, an analytical forward solution can be derived and the accuracy has be verified by the phantom dipole fitting experiments which shows there is no great effect to the result of source localization by MEG recordings [37]. At the following paragraph, it will describe the forward solution under the spherical head assumption. Involving the sphere head model, it divides the magnetic field b(r0 ) into two parts, radial and tangential components, the direct of which are radial and tangential to the sphere surface respectively. The radial component of the magnetic field at r can be compute as br (r) = b(r) ·. r , krk. (2.11). where krk is the magnitude of r. Because r/krk is the normal direction of the sphere surface, it can be found that the term in Eq. (2.3) caused by volume current will varnish. That is, the br (r) are just effect by the primary current: br (r) = b0 (r) ·. r µ0 r − r0 µ0 q r = (q × r × r0 · · )= . 0 3 krk 4π kr − r k |r| 4π krkkr − r0 k3. (2.12). The last equality above is due to (a · b) × a = a × b. This shows if the normal direction of the sensors is radial to the sphere surface, Eq. (2.12) is the solution. Oppositely, when the direction is not radial, tangential components which includes the effects of volume current must be considered. However, Sarvas [58] shows it has no necessary to explicitly compute the term associating volume currents. It exploits the radial magnetic field br (r) to derive a.

(33) 18. MEG Forward Prediction with Spherical Head Model. magnetic potential u(r) and then derive the full magnetic potential outside the head volume by the relation b(r) = −µ0 5 u(r). Summarizing form of the magnetic field outside the head sphere according to the Sarvas model is b(r) =. µ0 F (r, r)q × r0 − (q × r0 · r∇F (r, r)), 4πF 2 (r, r). (2.13). where ∇F is the gradient of the scalar function F which is F (r, r0 ) = kr − r0 k(krkkr − r0 k + krk2 − (r0 · r)).. (2.14). Mosher et al. [47] rearrange the solution of forward problem which utilizes the spherical head model to a pretty form by separating the linear and nonlinear parts of the dipole parameters. This reduces the computational cost for solving the inverse problem in practice, because the linear parameter has no necessary to find numerically. This form can be represented as b(r) = K(r, r0 ) · q,. (2.15). where q concludes the amplitude and orientation information of the dipole, whereas the parameter of dipole position is nonlinearly involved in the kernel matrix K(r, r0 ). For the case of Sarvas model in Cartesian coordinate system, K(r, r0 ) = [. µ0 ∇F r − F I Cr0 ], 4π F2. (2.16). where Crq equals .  0   r0  z . −rz0. ry0. 0. −rx0. rx0. 0. −ry0. .      . . and r0 = [ rx0 ry0 rz0 ]T . Figure 2.1 gives a depiction of the elementary elements and their relationship in the MEG forward model with spherical head model. When there are more than one sources, the resultant magnetic field is the sum of the effect by each individual source. That is, b(r) =. X ∀i. where i denotes the number of the dipoles.. K(r, r0 i ) · qi ,. (2.17).

(34) 2.3 Forward Solution for Spherical Head Model. 19. b (r , r ' , q ). K (r , r ' ). r q (j) r'. Figure 2.1: MEG forward model with spherical head model. (Left) When given an arbitrary source moment q locating at r0 with orientation j, it can derive the magnetic field at r outside the scalp by b = Kq. The head is modeled as the nested spheres with concentric, symmetric and homogenous property. (Right) It often exploits a equivalent current dipole (ECD) (big arrow) to represent a group of neighboring pyramidal neurons’ electrical activities (tiny arrows) in the brain. These neurons locate at the cortex area with the orientation perpendicular to the cortical surface. A MEG sensor is a coil to measure the magnetic flux through it. The measured magnetic flux is the integration of the partial component, parallel to the normal direction of the coil, of the magnetic field inside the circled area of the coil. Thus, the measurement of a sensor locates at r is m(r) =. Z. n·. m(r) =. ∀i. (2.18). K(r, r0 i ) dr ) · qi ),. (2.19). ∀i. A. X. K(r, r0 i ) · qi dr,. X. (( n ·. Z A. where n and A are the the normal vector and the area circled by of the coil respectively. We can concatenate the equations for each sensor to represent the relationship of the dipole source and the recordings measured by all MEG sensors, that is,  . . m  . 1    m2   .  .  . . mN. . ( n1 · K(r1 , r0 i ) dr1 ) · qi   A1.     X    =    ∀i     . R. ( n2 ·. R A2. ( nN ·. R. 0. K(r2 , r i ) dr2 ) · qi .. . K(rN , r0 i ) drN ) · qi.     ,     . AN. where N is the total number of the sensor. Because the parameters of the sensors is constant on the same MEG instrument, we can know these in advanced and can be omit they in the.

(35) 20. MEG Forward Prediction with Spherical Head Model. equation. Therefore, we rearrange the equation to be m=. X. G(r0 i ) · qi .. (2.20). ∀i. G is the gain matrix performing the transformation from a source to the MEG measurements. In practice, the integral part, for each sensor, in G can be approximated by the summation of the weighted magnetic fields at some discrete points inside the coil [49]. There are two issues about the forward model need to be concerned. First, it needs to be noted that the parameters of the positions and orientations of both sensors and source dipoles are all in the same coordinate system here, and so do they used in the following chapters. In practice, before the experiment, we need to derive the transformation from the coordinate system of head to that of MEG device and vice versa via the step of head position identification. Therefore, we can do the computation in either the spherical head or MEG device coordinate system. Second, we need to decide the center of the sphere, that is, the origin of the spherical head coordinate system. Usually, a good sphere center of the spherical head model can be determined by finding a sphere with most fit for the head geometry derived from MR image. So far, when giving ECD sources in the brain, we can explicitly describes the measurement on all MEG sensors with an analytical formulation. Figure 2.2 gives an example of a forward solution when one dipole is or two dipoles are given. With the forward solution, we can utilize it for solving the inverse problem. Next chapter, we start to introduce a kind of the inverse algorithm and which is the most major part in this thesis, that is, beamforming methods..

(36) 2.3 Forward Solution for Spherical Head Model. 21. +. Forw ard Prediction. Forw ard Prediction. Figure 2.2: An example of the forward solution. We simulate the putative dipole(s) (yellow mark) in the head and then the resultant magnetic field recorded by sensors can be derives by the calculation of forward solution. If there are more than one source, the solution is the addition of the individual’s. The magnetic filed topography is drawn by the interpolation of the sensors recordings estimated by Sarvas forward modal. The position of the sensors are unfolded to a plane for better visualization..

(37) 22. MEG Forward Prediction with Spherical Head Model.

(38) Chapter 3 Beamforming Methods for Functional Brain Imaging.

(39) 24. Beamforming Methods for Functional Brain Imaging. w1 w2. + Estimated signal. wN. S ensors. R ecordings. Figure 3.1: Basic concept of spatial filter. Spatial filter utilizes the recordings measured by sensors cites at distinct location to estimate the target source signal. It estimates the target source signal by summing the weighed recordings on each sensor. The kernel of the filter is the weighting coefficients designed by various criteria for different applications.. 3.1. Introduction. Beamforming method [12, 33] is a well-known technique to localize the signal source in the array signal processing field and its developments can be traced from mid of twenty century. Today it is widely applied in sonar, radar, and astronomical telescope array systems. It is a spatial filter which weights different values, called beamforming coefficients, to each recording measured by each sensor at distinct cites. It then sums the weighted recordings as the estimated signal of the target source. Figure 3.1 illustrates the basic concept of spatial filter. We can think beamforming method utilizes the information recorded by physical sensors to create a virtual sensor at a specified position in the source space to measure the signals with a specified orientation. The kernel of the method is how to determine the beamforming coefficients, which is based on different assumptions in different algorithms for different applications. The MEG inverse problem is a source localization problem with recordings measured by a set of sensors distributed outside the scalp. This is similar with the radar array system.

(40) 3.1 Introduction. 25. and can be modeled as the problem solved by beamforming technique. Recently, algorithms [21, 54, 59, 72] based on this technique are developed to apply on solving MEG inverse problem and there are more and more MEG-related studies [7, 16, 25, 60, 62, 65, 68] exploit these methods for their applications which include the analyses of motor-, somatosensory- and auditory-related experiments. Figure 3.2 depicts that the beamformer applies to radar array system and to MEG inverse problem. For detecting the brain source at a position with an orientation, the beamforming coefficients applying on each recordings of different MEG sensors is mainly determined by two criteria: retaining the targeted source exactly and, second, attenuating the effect cause by other sources. From the view of designing the filter, the spatially unit pass band of the filter is exactly at the specific position and orientation of the targeted source in one hand and, in the other hand, the filter tries its best to stop the contribution from other positions and orientations. These derives the filter with a unit gain constraint and minimum variance property. We will deduce the algorithm with much detail at the latter sections. By modifying the specified position and orientation, we can construct many virtual sensors distributing in the source space (brain) to image the electric activities in the brain. Beside the beamforming method, other common methods that we briefly mentioned in Chapter 1 for solving MEG/EEG inverse problem includes least square dipole fitting, MNE and MUSIC. Least square dipole fitting method nonlinearly search the optimal dipole sources’ parameters ,r0 i and qi according to minkm −. X. G(r0 i ) · qi k,. (3.1). ∀i. where m is the N × 1 vector to represent the measured data on the N sensors at an instant time. This method does well and fast when the dipole number is small, but if the number is large, not only the computational cost of the nonlinear searching will dramatically increase but the objective function may much far from unimodal, that is, there will be many local minimum. Moreover, the major drawback of the method is that it need to determine the source dipole number in advance, which is unsuitable for real applications. The MNE method assumes that the sources are distributed currents. For deriving the source parameters, it minimizes not only the error between the real recordings and the estimated one, like least square dipole fitting, but the size (norm) of the distributed currents. The drawback.

(41) 26. Beamforming Methods for Functional Brain Imaging. Figure 3.2: Example of the application of the beamforming technique: radar array system and MEG inverse problem. The upper figure shows that the radar array system uses the beamforming technique to put a virtue sensor (the small radar circled by the dash square) at a specify position in the source space (sky) for detecting if there is object (airplane) there. It can scan whole source space by shifting the virtue sensor. The same idea is adapted for solving the EEG/MEG inverse problem. That is, we put a virtue sensor (blue dash square) in the head to measure the signals at the specific position with the specific orientation and imaging the whole brain by moving the sensor through everywhere of the source space (brain). In both cases, the estimated signals (virtually measured by the virtual sensors) are constructed by summing the weighted recordings on physical sensors..

(42) 3.2 Source Estimation. 27. of the method is the solution will trend to the surface of the head because the effect of minimizing the norm of the source. MUSIC is a kind of spatial filter, which is designed by the criteria that will incline the estimated source orthogonal to the noise space which need to be estimated in advance. Using beamforming method to solve MEG inverse problem has no need to know to the number of the sources in advanced and the resultant estimated sources have no tendency to the surface of the head which is suffered by MNE approach. Furthermore, unlike MUSIC, it has no need to estimate the noise space beforehand. Robinson and Vrba [55] also mentioned that the beamforming-based method is better than MUSIC for the study based on unaveraged MEG recordings. In this chapter, we start from the deduction of the basics of beamforming-based method to imaging the brain source significance at section 3.2.1 and 3.2.2. Next, we point out the difficulty about the determination of the source orientation for probe and the disadvantage of others’ solution for it in the section 3.2.3. Consequently, at section 3.2.4, we propose a new method, maximum contrast beaformer, for solving the dilemma. After that, the DICS method which images the coherent sources is described at section 3.3.1 and 3.3.2 and its limitation. For solving it, we also proposed another new method, maximum normalized correlation beaformer, at section 3.3.3. Theoretically, the maximum normalized correlation beaformer can not only image the coherent sources in frequency domain but also in other domains if the auto- and cross-correlation are defined. We address the procedure of imaging the time-frequency coherent sources using proposed beamformer at the end of this chapter, section 3.3.4.. 3.2 3.2.1. Source Estimation Data Model. In time domain, bio-medical signal is often modeled as a random, or say chaotic, signal [5]. Random signal is the signal that we can not precisely predict its next features by its past history up to a significant level. The brain signal and the corresponding measured data by MEG are often modeled as the signal in this category. It needs to note that the model here is different with the source model mentioned in last chapter. They describe the temporal.

(43) 28. Beamforming Methods for Functional Brain Imaging. and electromagnetical characteristics of the brain signals respectively. In spatial domain, if the number of MEG sensors is N , the measured data can be view as a vector in a N dimension vector space that spanned by the basis of each sensor. The beamforming method solving MEG inverse problem spatially apples to the measured data of each sensor. That is, it performs a linear transform from the spatial domain and that is the reason why it is called ”spatial” filter. The unit-gain constraint and minimum variance beamformer uses the first and second order temporal statistic to feature both the measured data in the sensor-spanned space which equals the domain of the filter and the estimated source signal in the range of the filter. In the following, we first introduce the related notation for further deduction later. First, involving the Sarvas forward model, if there are k dipoles in head, at a certain instance, the recordings m equals the summation of the productions of each gain matrix and its corresponding dipole moment: m=. k X. G(r0 i ) qi ,. (3.2). i=1. where qi locating at r0 i is the dipole moment of the i-th in all k sources and G(r0i ) is the corresponding gain matrix. Then, m is the summation of the effects by every source dipoles in the brain. Without of losing generality, we represent the source dipole in the Cartesian coordinate system such that qi is a 3 × 1 vector and its elements represent x,y,z components of qi respectively: . .  m1      qix  m2     m =  .  , qi =   qiy  .    .    qiz mN.    ,  . and therefore, evidently, G(r0 i ) is a N × 3 matrix. In real environment, it still needs to consider the noise so we add the noise n to the above equation such that m=. k X. G(r0 i ) qi + n.. (3.3). i=1. We then separate the dipole moment into orientation and magnitude parts, and combine the.

(44) 3.2 Source Estimation. 29. former part and the gain matrix into `, the lead field, that is, m=. k X. G(r0 i ) ji qi + n =. k X. `(θi ) qi. (3.4). i=1. i=1. where ji = qi /kqi k is the orientation of the dipole source and we use θ to denotes both r0 and j. The subscript represents the number of the dipole as well. Consequently, we consider the signals as a sequence of realization rather than just a value at an instant time, that is, we adds the temporal information into the notations. The above equation needn’t modify, but the elements both in qi and in m is no longer deterministic scalars but random variables which represent the random processes of the measured and source signals respectively. Then, we can define the first and second order statistics both of these two kind signals. It denotes mean and variance of the magnitude of a dipole as q¯i = E{qi }. (3.5). cqi = E{ [qi − q¯i ][qi − q¯i ]T } ,. (3.6). respectively [17], where E performs the expectation. Similarly, the mean and covariance matrix of the measured data are ¯ = E{m} = m. k X. `(r0 i ) q¯i. (3.7). {`(θi ) cqi `T (θi )} + Cn ,. (3.8). i=1. ¯ ¯ T} = C = E{ [m − m][m − m]. k X i=1. respectively, where Cn is the covariance of the noise under a usual assumption of zero mean. In practice, it is online digitized for MEG recordings and which can be represented by a N × T matrix M, where T is the sampling numbers. With out of generality, we assumes the weight of each sample point is identical. Hence, the covariance matrix of the recordings can be estimated by 1 MM DF MTM DF , (3.9) N −1 is the mean-deviation form [36] of M, that is, every element in MM DF is C=. where MM DF. the corresponding element in M subtracting the arithmetic mean of the row at which the element lies..

(45) 30. Beamforming Methods for Functional Brain Imaging. 3.2.2. Spatial Filter Design. As mentioned in section 3.1, the virtual sensor created by spatial filter is constructed by means of summing the weighted recordings of each sensors and this can be written as y0 = wT (θ0 ) m,. (3.10). where y0 is the estimated scalar random signal at a specified position r00 and with the orientation j0 . w is the filter, the elements of which are the beamforming coefficients for weighting. The goal of the filter is that it can extract the targeted source with the parameter θ = θ0 and suppress other non-targeted (θ 6= θ0 ) interference. For achieving the goal, a unit-gain constraint and minimum variance spatial filter had be introduced, and its mathematical formula is min cy0 s.t. wT (θ0 ) `(θ0 ) = 1. w(θ0 ). (3.11). where cy0 is the variance of the estimated signal. The constraint in the formula performs the unit gain to targeted source for probe. It is because that if there is no noise and the forward model is exactly correct, then y0 = wT (θ0 ) m = wT (θ0 ) `(θ0 ) q0 = 1 × q0 = q0 ,. (3.12). where q0 is the true source signal at the target position and orientation. The minimization operator tries its best to reduce the gain from other sources. The concept of unit-gain constraint and minimum variance is depicted in Figure 3.3. After substituting Eq. (3.8) into Eq. (3.11), we can get min {wT (θ0 )Cw(θ0 )} s.t. wT (θ0 ) `(θ0 ) = 1. w(θ0 ). (3.13). The formulation is an equality-constrained optimization problem and can be solved by involving Langrange’s theorem [8]. We omit the arguments of the functions for clarity and than the theorem tells that the optimal solution of Eq. (3.13) satisfies ∇(wT Cw) + λ∇(wT ` − 1) = 0. (3.14). (wT ` − 1) = 0,. (3.15).

(46) 3.2 Source Estimation. 31. w T l=1. r', j. w T l=0. Figure 3.3: Spatial response of the unit-gain constraint and minimum variance spatial filter. For a specified position r0 with orientation j, the unit-gain and minimum variance spatial filter, first, is unit-gain at the position with the orientation and, second, tries its best to minimize the total gain from the source space. When the filter is ideal (real line), the gain is 1 (wT ` = 1) exactly at r0 with j, otherwise it is 0 (wT ` = 1). where ∇ denotes the gradient operator and after it operates, Eq. (3.14) becomes 2wT C + λ`T = 0.. (3.16). Rearranging and taking the transpose, CT w = Cw =. −1 λ`. 2. (3.17). The first equality in Eq. (3.17) is because of the symmetric property of covariance matrix. Next, multiplying C−1 on both sides of the last equality yields w=. −1 −1 λC ` 2. (3.18). Then, after multiplying `T on both sides of equal sign and according to the constraint in Eq. (3.13), we get −1 T −1 λ` C ` = 1 2 λ = −2 ( `T C−1 ` )−1. `T w =. (3.19) (3.20). Finally, substituting Eq. (3.20) into Eq. (3.18) derives the solution w = C−1 ` (`T C−1 ` )−1. (3.21). This solution will try its best to approach to the ideal spatial filter with perfect spatial specificity for θ0 , that is, the gain is zero as θ 6= θ0 , otherwise it is 1. However, at this.

(47) 32. Beamforming Methods for Functional Brain Imaging. time, the norm of the filter may be large such that it is very unstable and sensitive to the noise [72]. To comprise the spatial specificity and the noise sensitivity, a regularization tern is introduced to this optimizations problems [21, 54] and hence Eq. (3.13) and its corresponding solution becomes min {wT (θ0 )Cw(θ0 ) + αwT (θ0 )w(θ0 )} s.t. wT (θ0 ) `(θ0 ) = 1. (3.22). w = (C + αI)−1 ` (`T (C + αI)−1 ` )−1. (3.23). w(θ0 ). respectively, where α is the regularization parameter. A statistical mapping of the neural activities are often used to find the source significance [54, 72]. The reason is that, first, the norm of the lead field is inversely-squaredproportional to the dipole depth, which is relative to the center of head model [58], such that the norm of the filter is intrinsically larger when the source is deeper, and so does the estimated source and that, second, both the distribution of the noise and the sensitivity of the sensors are not uniform spatially. Therefore, to obtain a meaningful statistical brain activation, F-statistic F is calculated as F =. wT Cact w , wT Cctr w. (3.24). where Cact and Cctr are the covariance matrices estimated by the samples from active and control states respectively. The active and the control states are the periods of the brain activities in which we are interested and we are not. Here, we discuss the choice of C, Cact and Cctr . It is no doubt that we can estimate Cact directly from the recordings in the period when the interesting brain activities cite. Oppositely, Cctr can be estimated from the recordings in the period when the interesting brain activities do not cite or can exploit the empty room recordings for the estimation. If assuming the noise is uncorrelated, we can eliminate the elements in the non-diagonal part of Cctr . Moreover, If assuming the noise is independent identical distributed and the sensitivities of sensors are identical, we can use the identity matrix I as the control state covariance matrix. The period used to estimate C is the part to be suppressed its variance. If we set C the same with Cact , the meaning is that the filter will suppress the contribution from θ 6= θ0 according the statistic characteristic of the recordings in this period. Different.

(48) 3.2 Source Estimation. 33. combination of C, Cact and Cctr gives different meaning of the results. However, by our experience and other literatures, except Cact , we cannot concludes what is best for the choice of C and Cctr , that is, they still are options of the method. Another parameter in the method, α, which perform the regularization between noise sensitivity and spatial specificity is difficult to be automatically determined because the value is dependent on the SNR. Actually, we always cannot compute the SNR in advance because we have no idea about the signal which just can be estimated after applying the inverse method. In practice, it often uses recordings to guide a rough rule for setting the value. By our experience and other literatures, we often set the value vary from 1 to 0.00001 times the maximum eigenvalue of Cact . It needs to keep in mind that utilizing the maximum eigenvalue of Cact to derive an α is dominant by the maximum contribution of a certain source signal. That is, it is not able to concern the SNR of their sources, although there are no more proper ways yet.. 3.2.3. Source Space Scanning Approaches. So far, when given a position r0 and orientation j, it can use the filter to estimate the source at the specified position and the specified orientation. Ideally, it should be applied at everywhere with every orientation for imaging the activities in whole brain. Nevertheless, infinitely computation is impossible, so sampling is adapted. It consequently encounters a problem: how much is the enough sampling rate? If we leave the computational ability and the limitation of the resolution caused by the number of MEG sensors aside, according to the Nyquest sampling theorem [50], the spatial sampling period of the position to apply the beamforming filter should be smaller than the half of the distance of the neighboring pyramidal neurons so that we can see the electrical activities of the individual pyramidal neurons. Although such extremely high sampling rate is almost impossible to be adapted when practice, this still tells us that sampling the space for the position to be probed for the filter is the finer, the better because the cortex area is tightly full of the neurons from microscopic view. However, it is not true for the issue of the orientation, because, at one position, the pyramidal neuron(s) from which the most measurements of MEG are contributed is (are) just (equivalently) one orientation, perpendicular to the cortical surface. In.

(49) 34. Beamforming Methods for Functional Brain Imaging. the following, we discuss the issue about scanning the sources space by beamformer with the division of two categories, with and without the information of cortical surface geometry. Under each category, we discuss how it determines the position and orientation to be probed by beamformer respectively. Some studies [29,41] utilize CT or MR image to construct the cortical surface geometry and solve the MEG inverse problem subject to the constraint that the source site at the position under the surface about 2mm with the orientation perpendicular to it. In practice, it often models the surface of a 3D object by triangle meshes, so does the cortical surface model. Continuously, they applied the beamforming filter in accordance with the position of and the normal vector of the vertices on the triangle meshes to do the sampling both for position and orientation, which is illustrated in Figure 3.4. Although this approach seems good for the reason that it fits the neurophysiology knowledge, an accurate enough cortical surface model is necessary [29]. However, segmentation, one of the most difficult task in image processing [19], of the convoluted cortical surface from the 3D volume image data is not trivial [18]. It usually needs a great effort to achieve an accurate enough surface model. Without the cortical surface geometry, the other way simply samples the whole 3D space of the brain uniformly [21, 54, 59, 72] and there are different ways to determine the source orientation to probe in different studies. Van Veen et al. [72] applied the filter on three orthogonal axes, calculated the estimated source power on each axis and summed the results to be the final result as the source variance. They called the method linearly constrained minimum variance (LCMV) spatial filter and it was formulated as min {trace(WT (r0 0 )CW(r0 0 ))} s.t., WT (r0 0 ) G(r0 0 ) = I. W(r0 0 ). (3.25). where W is the a N × 3 matrix with that each column represents the beamforming coefficients for each axis. Probing the orientation solution space with few targeted direction like [72] is not proper because it is incorrect to think the source (fully) project to the axes. It has to sacrifice the spatial specificity by increasing the weight of regularization term α to avoid miss-detection, because the beamforming method does its best to suppress the signal other than the specified orientation. That is, when the filter approaches to ideal, we will find no contribution.

(50) 3.2 Source Estimation. 35. Figure 3.4: Determination of the position and orientation probed by beamformer when given the geometry model of the cortical surface. By the knowledge of neurophysiology, it is known that the electrical activities measured by MEG occur at cortical surface with the orientation perpendicular to it. After segmenting the cortical surface from the MR or CT image, it can construct and represent the cortical surface geometry by triangle meshes. Some studies apply the beamformer to the source space by sampling the position space at the position of the vertices on the mesh and, consequently, it determines the orientation as the vertex normal direction usually calculated by averaging the normal vectors of the connective triangle planes (right of the figure). The left of the figure is the resultant source power reconstructed by this approach for an finger lifting experiment [41]. The major difficulty of the approach is that accurately segmenting the convoluted cortical surface is not easy. from a source except that it exactly lies on one of the sampling axes. Granted that the α is not approach to zero, the result image map is flatter than the one with correct source orientation because some parts of the contribution by the source are still lost due to the power annulation tendency of the filter. This predicament occurs as well in the approaches using cortical surface geometry model when the model is not accurate enough. By a simulation, Figure 3.5 gives an example of miss-detection when the targeted orientation of beamformer is far aberrant from the orientation of the true source. Robinson and Vrba [54] proposed a method, ”synthetic aperture magnetometery,” (SAM). SAM searched the orientation on the plane tangential to the sphere of the head modal such that the value of z-deviate (performing the similar role to the F-statistic that mentioned in last section) is maximum, after sampled the orientation solution space. The method finds a meaningful orientation for probing and have no risk of missing the source if it indeed find the global maximum. Also, they [73] proved by deduction and showed by simulation that the result derived by SAM is more focal (opposite to flatter) than that derived by LCMV.