VBMEG Tutorial

VBMEG Tutorial

はじめに

VBMEGは、MEGやEEGのデータから皮質電流を推定するためのMatlabツールボックスである。 本チュートリアルは、サンプルデータの解析を通じて、VBMEGの処理の流れを把握することを目的としている。

準備するファイル

http://vbmeg.atr.jp/?lang=ja#DOWNLOADよりダウンロードできます。

  • MEGデータファイル(A008a_BC_0.5HPF_100LPF_trig437_label4.ave)
  • センサ位置合わせファイル(marker1.pos.mat)
  • MRI構造画像ファイル(3D.hdr, 3D.img)
  • 脳モデルファイル(lh.(curv/inflated/smoothwm).asc, rh.(curv.inflated/smoothwm).asc)
MEGデータファイル
  • 横河電機MEGデータファイル
  • 聴覚刺激課題
  • 3.2kHzの純音を、左耳に提示したときの脳活動を加算平均したもの。
センサ位置合わせファイルの詳細
  • MEGのセンサ座標を被験者のMRI構造画像座標系に位置合わせするための情報が入ったファイル。
  • 位置合わせプログラムにより作成する。
MRI構造画像ファイルの詳細
脳モデルファイルの詳細

ディレクトリ

  • 以下のようにプログラムおよび入力データが格納されていることを前提にしている。
D:\vbmeg(VBMEGプログラムディレクトリ)
D:\data(入力データディレクトリ)
│   3D.hdr
│   3D.img
│
├───FS
│       lh.curv.asc
│       lh.inflated.asc
│       lh.smoothwm.asc
│       rh.curv.asc
│       rh.inflated.asc
│       rh.smoothwm.asc
│
└───Yokogawa
        A008a_BC_0.5HPF_100LPF_trig437_label4.ave
        marker1.pos.mat

作業手順

パスを設定する

  • MATLABを起動し、実行準備をする。

setpath

プロジェクトを作成する

  1. コマンドラインからproject_mgrを起動する。

    >> project_mgr

  2. [File]->[New project]を選択する。

    Create new project

  3. プロジェクト名を入力し、Selectボタンを押す。

    Create new project

  4. カレントディレクトリをD:\に変更し、make dirボタンを押す。表示されたダイアログにAuditoryと入力し、OKボタンを押す。

    Create working directory

  5. 作業ディレクトリ(D:\Auditory)を選択し、APPLYボタンを押す。

    Create new project

  6. OKボタンを押し、ダイアログを閉じる。

    Create new project

  7. 確認ダイアログが出るので、「はい」を選択する。

    Confirm

  8. project_mgrの準備が整った。

    Ready for VBMEG

脳モデルをインポートする

  1. [Data Import]->[Import cortical model]を選択する。

    Launch import brain gui

  2. MRI構造画像ファイル, 脳モデルファイルをセットする(最初にラジオボタン:FreeSurferを押し、各Selectボタンを押して、ファイルを選択する)。

    FreeSurfer files

  3. 出力先を作成する(D:\Auditory\brain)
    • Selectボタンを押す。

    Make directory for output

    • make dirボタンを押すと、作成するディレクトリ名を入力するダイアログが表示されるので、brainと入力し、OKボタンを押す。

    Make directory for output

    • brainを選択してから、APPLYを押す。

    Make directory for output

    • 保存ファイル名が更新される。

    Save files

  4. Execボタンを押して、20分ぐらい待つと、脳モデルのインポートが完了し、結果が表示される。

    Import brain

    結果画面
    Import result

MEGデータをインポートする

  1. [Data Import]->[Import MEG data]->[Yokogawa]を選択する。

    Launch import meg data gui

  2. Yokogawa MEG fileを指定する。

    MEG data import gui

  3. ファイルダイアログの拡張子を.aveに変更する。

    MEG data import dialog

  4. Yokogawa MEGファイルを選択し、OKを押す(D:\data\Yokogawa\A008a_BC_0.5HPF_100LPF_trig437_label4.ave)

    MEG data import dialog

  5. Selectボタンを押し、センサ位置合わせファイルを選択する(D:\data\Yokogawa\marker1.pos.mat)。

    MEG data import dialog

  6. 出力先を作成する(D:\Auditory\meg)
    • Selectボタンを押す。

    Make directory for output

    • make dirボタンを押すと、作成するディレクトリ名を入力するダイアログが表示されるので、megと入力し、OKボタンを押す。

    Make directory for output

    • megを選択してから、APPLYを押す。

    Make directory for output

    • 保存ファイル名が更新される。

    Save files

  7. Execボタンを押す。10秒程度で、MEGデータのインポートが完了する。

    Import MEG data

    <<MATLAB Command window>>
    Import result

リードフィールドを計算する

  1. [Analysis]->[Calculate leadfield]を選択する。

    Launch leadfield calculation gui

  2. MEG fileを指定する(D:\Auditory\meg\A008a_BC_0.5HPF_100LPF_trig437_label4.meg.mat)

    meg file

  3. Cortical model fileを指定する(D:\Auditory\brain\3D.brain.mat)。

    cortical model file

  4. Cortical area fileを指定する(D:\Auditory\brain\3D.area.mat)。

    cortical area file

  5. 出力先を作成する(D:\Auditory\leadfield)。
    • Selectボタンを押す。

    Make directory for output

    • make dirボタンを押すと、作成するディレクトリ名を入力するダイアログが表示されるので、leadfieldと入力し、OKボタンを押す。

    Make directory for output

    • leadfieldを選択してから、APPLYを押す。

    Make directory for output

    • 保存ファイル名が更新される。

    Save files

  6. Execボタンを押す。20~30秒程度で、リードフィールドの計算が完了する。

    Calculate leadfield

    <<MATLAB Command window>>
    Calculate leadfield

電流分散を推定する

  1. [Analysis]->[Estimate current variance]を選択する。

    Launch VB estimation gui

  2. Cortical model fileを指定する(D:\Auditory\brain\3D.brain.mat)。

    Cortical model file

  3. Leadfield fileを指定する(D:\Auditory\leadfield\A008a_BC_0.5HPF_100LPF_trig437_label4.basis.mat)。

    Leadfield file

  4. MEG fileを指定する(D:\Auditory\meg\A008a_BC_0.5HPF_100LPF_trig437_label4.meg.mat)。

    MEG file

  5. Cortical area fileを指定する(D:\Auditory\brain\3D.area.mat)。

    Cortical area file

  6. Cortical activity fileを指定する(D:\Auditory\brain\3D.act.mat)。

    Cortical activity file

  7. Variance magnification parameterを10にする。

    Variance magnification parameter

  8. Dipole reduction ratioを0.2にする。

    Dipole reduction ratio

  9. 出力先を作成する(D:\Auditory\bayes)。
    • Selectボタンを押す。

    Make directory for output

    • make dirボタンを押すと、作成するディレクトリ名を入力するダイアログが表示されるので、bayesと入力し、OKボタンを押す。

    Make directory for output

    • bayesを選択してから、APPLYを押す。

    Make directory for output

  10. 保存ファイル名を入力する。

    auditory_left[ENTER]で、拡張子: .bayes.matは自動的に付く
    Save files

  11. Execボタンを押す。20分前後で推定が完了する。

    Current Variance Estimation

    <<MATLAB Command Window>>

    Noise model: full covariance
    Load MEG data for noise estimate [D:\Auditory\.\meg\A008a_BC_0.5HPF_100LPF_trig437_label4.meg.mat]
    Session 1: -499.2 - 0.0[ms]
    Number of basis files (= Number of sessions): 1
    Area ID: Cortex
    Number of vertices: 20004
    --- Reduce cortex
    Number of reduced vertices: 4004
    --- Spatial smoothing filter calculation
    R = 6.00e-003, Rmax= 1.20e-002
    Number of vertices = 4004
    Number of vertices in expanded area = 20002
    Basis file for session 1: D:\Auditory\.\leadfield\A008a_BC_0.5HPF_100LPF_trig437_label4.basis.mat
    Number of basis files (= Number of sessions): 1
    Area ID: Cortex
    Number of vertices: 20004
    --- Reduce cortex
    Number of reduced vertices: 4004
    --- Spatial smoothing filter calculation
    R = 6.00e-003, Rmax= 1.20e-002
    Number of vertices = 4004
    Number of vertices in expanded area = 20002
    Basis file for session 1: D:\Auditory\.\leadfield\A008a_BC_0.5HPF_100LPF_trig437_label4.basis.mat
    Number of sessions: 1
    Time window: [1 1250]
    MEG data file for session 1: D:\Auditory\.\meg\A008a_BC_0.5HPF_100LPF_trig437_label4.meg.mat
    Number of sensors           : 400
    Number of trials            : 1
    --- Check variable consistency is OK
    Noise model: spherical
    Load MEG data for noise estimate [D:\Auditory\.\meg\A008a_BC_0.5HPF_100LPF_trig437_label4.meg.mat]
    Session 1: -499.2 - 0.0[ms]
    Number of sessions: 1
    Time window: [1 625]
    MEG data file for session 1: D:\Auditory\.\meg\A008a_BC_0.5HPF_100LPF_trig437_label4.meg.mat
    Number of sensors           : 400
    Number of trials            : 1
    --- Initial VB-update iteration = 100
    --- Total update iteration      = 100
    Sensor noise variance is estimated
    Background current variance is estimated
    Set prior fMRI activity pattern

    --- New VBMEG estimation program ---

    --- Initial VB-update iteration = 1000
    --- Total update iteration      = 1000
    Tn=  1, Iter=  50, FE=16183629.600099, Error=5.920164e-002
    Tn=  1, Iter= 100, FE=16184064.896437, Error=5.940531e-002
    Tn=  1, Iter= 150, FE=16184112.004743, Error=5.944172e-002
    Tn=  1, Iter= 200, FE=16184140.401725, Error=5.944495e-002
    Tn=  1, Iter= 250, FE=16184143.926441, Error=5.944783e-002
    Tn=  1, Iter= 300, FE=16184144.942685, Error=5.944801e-002
    Tn=  1, Iter= 350, FE=16184145.363067, Error=5.944747e-002
    Tn=  1, Iter= 400, FE=16184145.588821, Error=5.944667e-002
    Tn=  1, Iter= 450, FE=16184145.749733, Error=5.944572e-002
    Tn=  1, Iter= 500, FE=16184145.890986, Error=5.944465e-002
    Tn=  1, Iter= 550, FE=16184146.025482, Error=5.944347e-002
    Tn=  1, Iter= 600, FE=16184146.152109, Error=5.944222e-002
    Tn=  1, Iter= 650, FE=16184146.263429, Error=5.944096e-002
    Tn=  1, Iter= 700, FE=16184146.351531, Error=5.943978e-002
    Tn=  1, Iter= 750, FE=16184146.413128, Error=5.943876e-002
    Tn=  1, Iter= 800, FE=16184146.451154, Error=5.943794e-002
    Tn=  1, Iter= 850, FE=16184146.472196, Error=5.943733e-002
    Tn=  1, Iter= 900, FE=16184146.482881, Error=5.943690e-002
    Tn=  1, Iter= 950, FE=16184146.487982, Error=5.943661e-002
    Tn=  1, Iter=1000, FE=16184146.490319, Error=5.943642e-002
    Alpha is scaled back by bsnorm
    ----- Save estimation result in D:\Auditory\.\bayes\auditory_left.bayes.mat

電流を計算する

  1. [Analysis]->[Estimate current]を選択する。

    Launch Cortical current calculation gui

  2. Current variance fileを指定する(D:\Auditory\bayes\auditory_left.bayes.mat)。

    自動的にMEG fileはセットされる。
    Current variance file

  3. 出力先を作成する(D:\Auditory\current)。
    • Selectボタンを押す。

    Make directory for output

    • make dirボタンを押すと、作成するディレクトリ名を入力するダイアログが表示されるので、currentと入力し、OKボタンを押す。

    Make directory for output

    • currentを選択してから、APPLYを押す。

    Make directory for output

    • 保存ファイル名が更新される。

    Save files

  4. Execボタンを押して、1分程度で電流計算が完了する。

    Current calculation

    <<MATLAB Command Window>>

    ----- New VBMEG -----
    Start current estimation
    Number of sessions: 1
    Time window: [1 1250]
    MEG data file for session 1: D:\Auditory\.\meg\A008a_BC_0.5HPF_100LPF_trig437_label4.meg.mat
    Number of sensors           : 400
    Number of trials            : 1
    --- Lead field matrix of focal window 
    Number of basis files (= Number of sessions): 1
    Area ID: Cortex
    Number of vertices: 20004
    --- Reduce cortex
    Number of reduced vertices: 4004
    --- Spatial smoothing filter calculation
    R = 6.00e-003, Rmax= 1.20e-002
    Number of vertices = 4004
    Number of vertices in expanded area = 20002
    Basis file for session 1: D:\Auditory\.\leadfield\A008a_BC_0.5HPF_100LPF_trig437_label4.basis.mat
     --- Save estimated current: D:\Auditory\.\current/auditory_left.curr.mat

推定電流を見る

  1. project_mgr GUIを表示し、vb_job_currentを選択する。

    すると、右上にポップアップメニューが表示されるので、View estimated currentを選択する。
    Launch Cortical current viewer

  2. 推定電流を表示するためのGUIが表示される。画面左上の空間分布表示パネルには、選択されている時間窓(time of interest; TOI)で平均された電流強度分布が皮質モデル(inflated model)上に表示される。画面下部の電流時系列表示パネルには、選択されている関心領域(region of interest; ROI)に含まれる頂点に渡り平均された電流強度時系列が表示される。ROIのサイズはデフォルトで半径4mmである。

    Cortical current viewer

  3. TOIは"Time window"エディットテキストで設定できる。ここではTOIを、左右両側に表れる聴覚関連磁場(AEF)の潜時である、80~120msecに設定する。

    Set timewindow

    TOIの設定に伴い、空間分布表示パネル上の電流強度空間分布が更新される。TOIは時系列表示パネル上で次のように灰色で表示される。
    Set timewindow

  4. Leftボタンを押すと皮質モデルを左側から見ることができる。聴覚野付近に脳活動が見られる。

    View point left

  5. ROIの中心を選択するために、"Select vertex"ボタンを押してから皮質モデル上の1点を選択する。ここでは左脳聴覚野付近の脳活動源を選択する。

    Choose activity

    "Vertex index"エディットボックスに選択された頂点のインデックスが表示される。
    Vertex index

    このエディットボックスに頂点番号を直接入力する事でもROIの中心を変更できる。
    ROIの中心は皮質モデル上で緑の×印で示される。電流時系列表示パネルには選択されたROIの電流時系列が表示される。Y軸のスケールは電流時系列に合わせて変わる。
    Select vertex

  6. より電流強度の高い頂点を探し出すために、"Spatial peak"ボタンを押す。するとROIの中で電流強度の時間平均が最大の頂点が探索され、ROIの中心がその頂点に設定される。

    Search spatial peak

    同時に"Vertex index"エディットボックスの頂点番号が更新される。
    Vertex index

    電流時系列表示パネルに更新されたROIの電流時系列が表示される。
    Current timecourse

  7. 同様の手順で右脳聴覚野付近の推定電流を見ることができる。

    Select vertex

    Current timecourse

推定電流を見る(補足)

  • 右脳聴覚野付近の脳活動源の前方に、もう一つ脳活動源が見られる。

    Select vertex

    Current timecourse

  • これら二つの脳活動源の時系列を比較すると、一方で電流が内向き(時系列表示上で正)のとき、もう一方では外向き(負)になっている。これは、一つの脳活 動源が二つの脳活動源に分離されたものと考えられる。このような推定脳活動源の分離は、二つの脳活動源(頂点)のリードフィールドが互いに反転したパター ンを持つときに起こり得る。リードフィールドの吸い込みと湧き出しの空間位置は頂点位置によって決まるので、二つの頂点は互いに近い位置にあると予想され る。さらに、リードフィールドの正負、すなわち吸い込みと湧き出しの向きは頂点位置での皮質法線方向で決まるので、二つの頂点位置ではそれらの法線方向が 互いに逆向きになっていると予想される。

    二つの要因のうち、頂点位置についてGUI上で定性的に確認する。"Inflated model"表示では元の皮質モデル(白質/灰白質境界)の位相情報しか保たれていない。元の皮質モデルを見るために、表示を"Folded model"に変更する。
    Original model

    すると元の皮質モデル上が表示される。
    Original model

    右脳聴覚野付近の脳活動源は皺に隠れている。表示角度を変える事でこの脳活動を見ることができる。
    Original model

  • リードフィールドを決定する二つの要因である頂点位置と法線方向について定量的に調べるため、MATLABコンソール上で作業を行う。次のコマンドによって、頂点位置と法線方向が変数"V"と"xx"にそれぞれ代入される。
    >> [V,F,xx] = vb_load_cortex('3D.brain.mat'); 
    
  • 二つの頂点間の距離と法線ベクトル方向余弦は次のようにして得られる。
    sqrt(sum((V(19805,:)-V(19617,:)).^2))
    
    ans =
    
        0.0141
    
    >> sum(xx(19805,:).*xx(19617,:))             
    
    ans =
    
       -0.7452
    
    >> 
    

    以上の通り、頂点間距離は14.1cm、方向余弦は-0.745であった。これより二つの脳活動源に対応するリードフィールドが互いに反転したパターンを持つ事が示唆される。