img-thumbnail

sn42

R&D Job in Japan.

Topics
My Twitter

fMRI解析の基礎 (８)：接続性解析

最終更新日: 2020年7月18日
fMRI解析の基礎シリーズ

0. はじめに
1. BOLD反応
2. 実験デザイン
3. 前処理
4. GLMの理論
5. GLMと検定
6. 多重比較問題
7. 集団解析
8. 接続性解析
9. コヒーレンス解析
10. 主成分分析
11. 独立成分分析
12. デコーディング
13. エンコーディング
14. 動的因果モデリング
15. 表現類似性分析


ここでは、結合性解析とよばれる、あるタスクによって機能的に接続される神経回路を特定する解析を行います。

この解析は、GLMを用いたクラスタ解析と似ているようで違います。GLMは「タスクに関連するクラスタ」を特定するのに対し、結合性解析は「タスクに関連するネットワーク」を特定するのです。

fMRIの文脈では、機能的結合性 と 効果的結合性 を区別します。

  • 機能的結合性 (Functional Connectivity)
    • 統計的な領域間の関連性
    • 相関係数など
  • 効果的結合性 (Effective Connectivity)
    • 神経回路モデルに基づいた領域間の接続性
    • Dynamic Causal Modeling など

機能的結合性を測定する最も簡単な方法は、領域間でピアソンの相関係数を計算することです。ただし、これには以下のような問題があります。

  • ２つの領域が、異なる HRF (血流動態反応関数) に従う可能性がある。

→ 同じ神経活動でも相関係数が低くなってしまう。

  • ２つの領域が、独立したノイズを含んでいる。

→　同じ真の信号でも相関係数が低くなってしまう。

ここでは、相関係数を計算するより効果的にネットワークを検出する方法について考えていきます。

  1. 心理生理相互作用 (PPI)
  2. ベータ系列回帰
  3. グレンジャー因果関係

1. 心理生理相互作用 (PPI)

PPI分析では、「タスク時の相関係数」と「非タスク時の相関係数」の差を利用して、タスクに関連するネットワークを同定します (Friston et al. (1997), Gitelman et al. (2003))。

Psychophysiological Interactions（心理生理相互作用）の名前の由来は、タスクを識別する変数が心理学的変数と生理学的変数の相互作用であったからです。

まずは、解析する領域（seed領域）を決めます。全脳のボクセル間の組み合わせを考えると計算時間が膨大となってしまうので、関心のある領域を含むようにseed領域を選択します。seed領域の信号はひとつの信号になるように処理しておきます。

(例)：Finger-Tapping タスクの場合、一次運動野をseedとして他のボクセルとの接続性を調査する。一次運動野の信号も、GLMの結果を利用して効果的に決定する。

Block Design

  1. 神経活動 N(t)N(t) を仮定する。
  2. N(t)N(t) にHRFをたたみ込んで、予測BOLD yp(t)y_p(t) を算出する。
  3. seed領域の信号 ys(t)y_s(t) を算出する。
  4. 乗数 ym(t)=2yp(t)−1y_m(t) = 2y_p(t) - 1 を算出する。
  5. 偏差 yds(t)=ys(t)−yˉsy_{ds}(t) = y_s(t) - \bar{y}_s を算出する。
  6. 乗数と偏差から、PPI変数 yppi(t)y_{ppi}(t) を算出する。
yppi(t)=ym(t)×yds(t)y_{ppi}(t) = y_m(t) \times y_{ds}(t)
  1. 次のGLMを立てて、seed以外のボクセルに回帰させる。
y=θ1yp+θ2ys+θ3yppi+B0+nuisance+εy = \theta_1 y_p + \theta_2 y_s + \theta_3 y_{ppi} + B_0 + \text{nuisance} + \varepsilon

このモデルにおいて、yy が ypy_p と ysy_s のみで説明されるなら θ3\theta_3 は 00 となるため相互作用は無いと考えられます。そのため、θ^1,θ^2\hat{\theta}_1, \hat{\theta}_2 が有意に大きければ、心理学的・生理学的と別々にseedと相関している、つまり機能的に接続しているといえるかもしれません。しかし、これは同じ交絡要因の入力による偽相関の可能性があります。

そこでPPI変数 yppi(t)y_{ppi}(t) を利用します。乗数はタスク時に+1+1となり、非タスク時に−1-1となります。そのため、ysy_s とタスク時に正の相関があり、非タスク時に負の相関があれば、PPI変数は相関するわけです。したがって、θ^3\hat{\theta}_3 が大きい場合にタスク＞非タスクの相関があり、機能的に接続しているといえます。

問題点は、非タスク時にseedとtargetに負の相関がある場合に結合性を特定しやすいといった、余計な仮定を含んでいる点です。そのため、非タスク時に相関が無い場合の検出力は下がってしまいます。

Rapid Event-Related Design

Rapid Event-Related Design では、 神経活動 N(t)N(t) と 予測BOLD yp(t)y_p(t) の差が大きいため、PPI変数を生成する乗数に時間的なずれが発生する可能性があります。これを緩和するために、実際のBOLD ys(t)y_s(t) を利用してずれを補正することができます (Gitelman et al. (2003))。

ただし、Block Design より大きく検出力が下がります。これは、BOLDが飽和することが少ないために、相関を上手く検出できないからです。

2. ベータ系列回帰

ベータ系列回帰では、全てのイベントを別々の説明変数として、GLMで重み β\beta を推定します。さらに、２つの領域間で一連の β\beta の相関係数を測定し、機能的結合性を評価します。

つまり、イベントごとに神経活動の度合いが協調して変動する領域間が機能的に接続していると考えます。

β⃗seed=(βs,event1,⋯ ,βs,eventN)β⃗target=(βt,event1,⋯ ,βt,eventN) FC=corr(β⃗seed,β⃗target)\vec{\beta}_{\text{seed}} = (\beta_{\text{s,event1}}, \cdots, \beta_{\text{s,eventN}}) \\ \vec{\beta}_{\text{target}} = (\beta_{\text{t,event1}}, \cdots, \beta_{t,\text{eventN}}) \\ \ \\ \text{FC} = \text{corr}(\vec{\beta}_{\text{seed}}, \vec{\beta}_{\text{target}})

ただし、ベータ系列は各イベントについて独立に回帰係数を測定するため、時系列情報が失われることになります。そのため、毎回同様な β\beta を示す領域間はその値が高くても相関係数が低くなります。この場合は、FBR法を用いたGLMを利用することで時系列情報を取り入れることができます。

(Example)

Correlation GLM:β⃗=(3, 3, 3)FBR-based GLM:β⃗=(1,3,2, 1,3,2, 1,3,2)\begin{array}{rl} \text{Correlation GLM:}& \vec{\beta} = (3,\ 3,\ 3) \\ \text{FBR-based GLM:} & \vec{\beta} = (1,3,2,\ 1,3,2,\ 1,3,2) \end{array}

ベータ系列回帰の利点は、イベントの種類ごとに分けて機能的結合性が測定できる点です。課題としては、β\beta の推定が正確でなければならないことが考えられます。

3. グレンジャー因果関係

経済学による因果推論方法を応用して、脳領域間の因果関係を推定します。中でもグレンジャー因果関係は、機能的結合性を評価するために広く利用されています。

領域A→領域B\boxed{\text{領域A}} \to \boxed{\text{領域B}}

グレンジャー因果推定では、領域Aが領域Bを活性化させているかどうかを推定できます。このとき、領域Bの信号の未来予測が、領域Aの信号を利用することで予測精度が向上する場合、A→Bの因果があると考えます。

if the causality A→B :(A + B⟶predictB)>(B⟶predictB)\begin{array}{l} \text{if the causality A} \to \text{B :} \\ (\text{A + B} {\overset{\text{predict}}{\longrightarrow}} \text{B}) \gt (\text{B} {\overset{\text{predict}}{\longrightarrow}} \text{B}) \end{array}

これを実行するためには、「予測」という意味を詳細に決める必要があります。グレンジャー因果推定では、線形自己回帰モデル (Linear AR model) を利用します。以下の AR(p) 過程では、pステップ前の時間までの信号を利用して予測を行います。

y(t)=a1y(t−1)+a2y(t−2)+⋯+apy(t−p)+ε(t)y(t) = a_1y(t-1) + a_2y(t-2) + \cdots + a_py(t-p) + \varepsilon(t)

誤差項 ε(t)\varepsilon(t) が平均 00 分散 σ2\sigma^2 の正規分布に従うとしたとき、誤差の分散 σ2\sigma^2 が予測精度と解釈できます。もし、任意の時間 tt の信号 y(t)y(t) を完全に予測できるのであれば、誤差の分散 σ2\sigma^2 は 00 となるからです。

3.1. 因果関係の定量的測定

領域A→領域Bの因果関係を調べるために、以下の２つのARモデルを作ります。

yB(t)=∑k=1paBk yB(t−k)+εB(t)yB(t)=∑k=1p[dAk yA(t−k)+dBk yB(t−k)]+εB|AB(t)\begin{array}{l} y_\text{B}(t) = \sum_{k=1}^p a_{\text{B}k} \ y_\text{B}(t-k) + \varepsilon_\text{B} (t) \\ y_\text{B}(t) = \sum_{k=1}^p \left[ d_{\text{A}k} \ y_\text{A}(t-k) + d_{\text{B}k} \ y_\text{B}(t-k) \right] + \varepsilon_\text{B|AB} (t) \\ \end{array}

ここで、

εB(t)∼N(0,σB2)εB|AB(t)∼N(0,σB|AB2)\varepsilon_\text{B} (t) \sim \mathcal{N}(0,\sigma_\text{B}^2) \\ \varepsilon_\text{B|AB} (t) \sim \mathcal{N}(0,\sigma_\text{B|AB}^2)

です。予測精度を比較するためには σB2\sigma_\text{B}^2 と σB|AB2\sigma_\text{B|AB}^2 を比較すればよく、σB2>σB|AB2\sigma_\text{B}^2 \gt \sigma_\text{B|AB}^2 のとき、因果関係があるといえます。

統計的に複雑なのは、２つのモデルでパラメータの数（パラメータの自由度）が異なることです。そのため、自由度が低いモデルが正しいという帰無仮説を検定することになります。

また、相互に活性化するかどうかを検証するために、次のモデルも必要となります。

yA(t)=∑k=1pbAk yA(t−k)+εA(t)yA(t)=∑k=1p[cAk yA(t−k)+cBk yB(t−k)]+εA|AB(t)\begin{array}{l} y_\text{A}(t) = \sum_{k=1}^p b_{\text{A}k} \ y_\text{A}(t-k) + \varepsilon_\text{A} (t) \\ y_\text{A}(t) = \sum_{k=1}^p \left[ c_{\text{A}k} \ y_\text{A}(t-k) + c_{\text{B}k} \ y_\text{B}(t-k) \right] + \varepsilon_\text{A|AB} (t) \\ \end{array}

ここで、

εA(t)∼N(0,σA2)εA|AB(t)∼N(0,σA|AB2)\varepsilon_\text{A} (t) \sim \mathcal{N}(0,\sigma_\text{A}^2) \\ \varepsilon_\text{A|AB} (t) \sim \mathcal{N}(0,\sigma_\text{A|AB}^2)

です。もし、領域Aが領域Bを活性化するが逆は起こらないならば、σB2>σB|AB2\sigma_\text{B}^2 \gt \sigma_\text{B|AB}^2 かつ σA2≯σA|AB2\sigma_\text{A}^2 \ngtr \sigma_\text{A|AB}^2 となります。

次に、これらの適合度の比較を考えます。Geweke (1982) によれば、次の統計量 FA,BF_{\text{A,B}} によって領域Aと領域Bの線形関係を捉えることができます。

FA,B=FA→B+FB→A+FA⋅Bwhere:   Fi→j=ln⁡(σj2/σj∣ij2)  Fi⋅j=ln⁡[(σi∣ij2σj∣ij2)/(σi∣ij2,σj∣ij2−cov(εi∣ij,εj∣ij))]\begin{array}{l} F_{\text{A,B}} = F_{\text{A} \to \text{B}} + F_{\text{B} \to \text{A}} + F_{\text{A} \cdot \text{B}} \\ \text{where: }\\ \ \ F_{i \to j} = \ln (\sigma_j^2/\sigma_{j|ij}^2) \\ \ \ F_{i \cdot j} = \ln \left[(\sigma_{i|ij}^2 \sigma_{j|ij}^2)/(\sigma_{i|ij}^2, \sigma_{j|ij}^2 - \text{cov}(\varepsilon_{i|ij}, \varepsilon_{j|ij})) \right] \end{array}

3.2. パラメータ推定

全ての AR(p) 過程は以下のように行列形式で記述できます。

[yA(NTR)⋮yA(p+1)]=[yA(NTR−1)⋯yA(NTR−p)⋮⋱⋮yA(p)⋯yA(1)][aA1⋮aAp]+[εA(NTR)⋮εA(p+1)] ∴ yA=YAaA+εA\begin{array}{rcl} \begin{bmatrix} y_\text{A}(N_{TR}) \\ \vdots \\ y_\text{A}(p+1) \end{bmatrix} &=& \begin{bmatrix} y_\text{A}(N_{TR}-1) & \cdots & y_\text{A}(N_{TR}-p) \\ \vdots & \ddots& \vdots \\ y_\text{A}(p) & \cdots & y_\text{A}(1) \end{bmatrix} \begin{bmatrix} a_{\text{A}1} \\ \vdots \\ a_{\text{A}p} \end{bmatrix} + \begin{bmatrix} \varepsilon_\text{A}(N_{TR}) \\ \vdots \\ \varepsilon_\text{A}(p+1) \end{bmatrix} \\ \ \\ \therefore \ \boldsymbol{y}_\text{A} &=& Y_\text{A} \boldsymbol{a}_\text{A} + \boldsymbol{\varepsilon}_\text{A} \end{array}

aA\boldsymbol{a}_\text{A} は、ガウスマルコフ定理を用いて推定できます。

a^A=(YA⊤YA)−1YA⊤yA\hat{\boldsymbol{a}}_\text{A} = (Y_{\text{A}}^{\top} Y_\text{A})^{-1} Y_{\text{A}}^{\top} \boldsymbol{y}_\text{A}

この推定値 a^A\hat{\boldsymbol{a}}_\text{A} を利用して、誤差項とその分散を推定します。

ε^A=yA−YAa^Aσ^A2=(ε^A⊤ε^A)/(NTR−p)\begin{array}{l} \hat{\boldsymbol{\varepsilon}}_\text{A} = \boldsymbol{y}_\text{A} - Y_{\text{A}} \hat{\boldsymbol{a}}_\text{A} \\ \hat{\sigma}_\text{A}^2 = (\hat{\boldsymbol{\varepsilon}}_\text{A}^{\top} \hat{\boldsymbol{\varepsilon}}_\text{A})/(N_{TR}-p) \end{array}

また、領域AとBを組み合わせた場合は以下のようになります。

yA=[YA,YB][cAcB]+εA|AB\boldsymbol{y}_\text{A} = \left[ Y_\text{A}, Y_\text{B} \right] \begin{bmatrix} \boldsymbol{c}_\text{A} \\\boldsymbol{c}_\text{B} \end{bmatrix} + \boldsymbol{\varepsilon}_\text{A|AB}

これも同様に、分散・共分散を推定できます。

ε^A|AB=yA−[YA,YB][c^Ac^B]σ^A|AB2=(ε^A|AB⊤ε^A|AB)/(NTR−p)cov^AB=(ε^A|AB⊤ε^B|AB)/(NTR−p)\begin{array}{l} \hat{\boldsymbol{\varepsilon}}_\text{A|AB} = \boldsymbol{y}_\text{A} - \left[ Y_\text{A}, Y_\text{B} \right] \begin{bmatrix} \hat{\boldsymbol{c}}_\text{A} \\ \hat{\boldsymbol{c}}_\text{B} \end{bmatrix} \\ \hat{\sigma}_\text{A|AB}^2 = (\hat{\boldsymbol{\varepsilon}}_\text{A|AB}^{\top} \hat{\boldsymbol{\varepsilon}}_\text{A|AB})/(N_{TR}-p) \\ \hat{\text{cov}}_{\text{AB}} = (\hat{\boldsymbol{\varepsilon}}_\text{A|AB}^{\top} \hat{\boldsymbol{\varepsilon}}_\text{B|AB})/(N_{TR}-p) \end{array}

ppはBICを最小化するように決定できます。

BIC(p)=(NTR−p)ln⁡σ^2+pln⁡(NTR−p)BIC(p) = (N_{TR}-p)\ln \hat{\sigma}^2 + p \ln(N_{TR}-p)

3.3. 推論

領域Aが領域Bを活性化させるかどうかを検定するためには、帰無仮説 H0\text{H}_0 の検定をします。

H0:FA→B=0H0:FA→B>0\text{H}_0: F_{\text{A} \to \text{B}} = 0 \\ \text{H}_0: F_{\text{A} \to \text{B}} \gt 0

Geweke (1928) は、帰無仮説が真であるならば、統計量 NTRF^A→BN_{TR} \hat{F}_{\text{A} \to \text{B}} は 自由度 pp のカイ二乗分布に漸近的に従うことを示しました。

NTRF^A→B∼χ2(p)NTRF^B→A∼χ2(p)NTRF^A⋅B∼χ2(1)∴ NTRF^A,B∼χ2(2p+1)\begin{array}{l} N_{TR} \hat{F}_{\text{A} \to \text{B}} \sim \chi^2(p) \\ N_{TR} \hat{F}_{\text{B} \to \text{A}} \sim \chi^2(p) \\ N_{TR} \hat{F}_{\text{A} \cdot \text{B}} \sim \chi^2(1) \\ \therefore \ N_{TR} \hat{F}_{\text{A}, \text{B}} \sim \chi^2(2p+1) \end{array}

したがって、カイ二乗検定を行って領域Aと領域Bの因果関係を推定します。

3.4. 条件付きグレンジャー因果関係

もしある領域３が領域１・２を活性化する場合、領域１と２の間の因果関係をどのように推定すればよいでしょうか？

領域３の影響を取り除いて因果推論するためには、以下のような条件付きグレンジャー因果関係 F1,2∣3F_{1,2|3} を考えます。

F1→2∣3=ln⁡(σ2∣232/σ2∣1232)F2→1∣3=ln⁡(σ1∣132/σ1∣1232)F1⋅2∣3=ln⁡[(σ1∣1232σ2∣1232)/(σ1∣1232,σ2∣1232−cov(ε1∣123,ε2∣123))]∴ F1,2∣3=F1→2∣3+F2→1∣3+F1⋅2∣3\begin{array}{l} F_{1 \to 2|3} = \ln (\sigma_{2|23}^2 / \sigma_{2|123}^2) \\ F_{2 \to 1|3} = \ln (\sigma_{1|13}^2 / \sigma_{1|123}^2) \\ F_{1 \cdot 2|3} = \ln \left[(\sigma_{1|123}^2 \sigma_{2|123}^2)/(\sigma_{1|123}^2, \sigma_{2|123}^2 - \text{cov}(\varepsilon_{1|123}, \varepsilon_{2|123})) \right] \\ \therefore \ F_{1,2|3} = F_{1 \to 2|3} + F_{2 \to 1|3} + F_{1 \cdot 2|3} \end{array}

また次の式が成立します。

F1,2∣3=F13→2−F3→2where F13→2=ln⁡(σ22/σ2∣1232)\begin{array}{l} F_{1,2|3} = F_{13 \to 2} - F_{3 \to 2} \\ \text{where } F_{13 \to 2} = \ln (\sigma_{2}^2 / \sigma_{2|123}^2) \end{array}

この結果から、「１・３ →\to ２の作用から３ →\to ２の作用を引けば、１と２の線形関係が分かる」と解釈できます。

3.5. 理論的拡張

グレンジャー因果推論の拡張としては、次のようなものがあります。

3.6. 評価

fMRIではBOLD信号は利用できますが、神経活動自体は利用できません。次のような疑問があると思います。

  1. 時間解像度は十分か？

→ TR = 250 ms で行った研究では、左右視覚野の112 msの刺激差による因果推論が可能でした (Rogers et al. (2010))。これほど高い時間解像度は無いことが多いですが、離れた領域間であれば、神経活動の遅れを考慮してグレンジャー因果推論が十分できると考えられています。

  1. HRFの違いが因果推論に影響するのか？

→ HRFの違いは完全にグレンジャー因果推論に影響します。解決策としては、集団解析としてHRFの変動を抑えることや、コントロール条件を上手く利用すること、逆たたみ込みを利用することが挙げられます。

Reference

  • Ashby, F. G. (2019). Statistical analysis of fMRI data. MIT press. url
  • 多次元神経イメージングによる脳内ネットワーク解析 url
  • Friston, K. J., Buechel, C., Fink, G. R., Morris, J., Rolls, E., & Dolan, R. J. (1997). Psychophysiological and modulatory interactions in neuroimaging. Neuroimage, 6(3), 218-229. url
  • Gitelman, D. R., Penny, W. D., Ashburner, J., & Friston, K. J. (2003). Modeling regional and psychophysiologic interactions in fMRI: the importance of hemodynamic deconvolution. Neuroimage, 19(1), 200-207. url
  • Marinazzo, D., Liao, W., Chen, H., & Stramaglia, S. (2011). Nonlinear connectivity by Granger causality. Neuroimage, 58(2), 330-338. url
  • Chen, Y., Bressler, S. L., & Ding, M. (2006). Frequency decomposition of conditional Granger causality and application to multivariate neural field potential data. Journal of neuroscience methods, 150(2), 228-237. url
  • Valdes-Sosa, P. A. (2004). Spatio-temporal autoregressive models defined over brain manifolds. Neuroinformatics, 2(2), 239-250. url
  • Rogers, B. P., Katwal, S. B., Morgan, V. L., Asplund, C. L., & Gore, J. C. (2010). Functional MRI and multivariate autoregressive models. Magnetic resonance imaging, 28(8), 1058-1065. url
Topics
My Twitter