img-thumbnail

sn42

R&D Job in Japan.

Topics
My Twitter

ボルテラモデルでBOLD反応をモデル化する

最終更新日: 2020年7月26日

Jupyter Notebookはこちら


非線形なBOLD反応

神経活動と共に血流流量が増加することで、BOLD反応が重ね合わせの原理を破ることがあります。重ね合わせの原理を破るということは、線形関係の成り立たない非線形な反応が生じるということです。

非線形な応答をモデル化した有名な例は バルーンモデル (Balloon model; Buxton et al. (1998)) であり、神経活動からBOLD反応が生成する過程を詳細にモデル化することで非線形な反応をモデル化しました。

しかし、バルーンモデルは複雑であるために計算量が膨大となってしまいます。そこで、非線形性だけに注目し、非線形系理論を応用すると計算量を削減できます。ここでは、ボルテラモデル (Volterra model) という非線形系のモデルを利用してBOLD反応をモデル化します。

ボルテラモデル (Volterra model)

時間不変なシステムは、入出力の関係をボルテラ級数で表現できます。ここで、hkh_k を kk 次ボルテラカーネルといいます。

y(t)= h0+∫h1(τ1)x(t−τ1)dτ1+∫∫h2(τ1,τ2)x(t−τ1)x(t−τ2)dτ1dτ2+⋯\begin{alignedat}{2} y(t) = \ & h_0 + \int h_1(\tau_1) x(t-\tau_1) d\tau_1 \\ &+ \int \int h_2(\tau_1, \tau_2) x(t-\tau_1) x(t-\tau_2) d\tau_1 d\tau_2 + \cdots \end{alignedat}

ボルテラ級数を用いて、BOLD反応 y(t)y(t) と神経活動 N(t)N(t) の関係を表現します。

y(t)=∑i=1∞Hi[N(t)]whereHi[N(t)]=∫⋯∫hi(τ1,…,τ2)N(t−τ1)⋯N(t−τi)dτ1…dτi\begin{alignedat}{2} y(t) &= \sum_{i=1}^\infty H_i[N(t)] \\ \text{where} \\ & H_i[N(t)] = \int \cdots \int h_i(\tau_1, \dots, \tau_2) N(t-\tau_1) \cdots N(t-\tau_i) d\tau_1 \dots d\tau_i \end{alignedat}

Friston et al. (1998) は、非線形なBOLD応答は次の式で近似できるとしました。

y(t)= H1[N(t)]+H2[N(t)]= ∫0tN(τ1)h1(t−τ1)dτ1+∫0t∫0th2(τ1,τ2)N(t−τ1)N(t−τ2)dτ1dτ2　\begin{alignedat}{2} y(t) = \ & H_1[N(t)] + H_2[N(t)] \\ = \ & \int_0^t N(\tau_1) h_1(t- \tau_1) d\tau_1 \\ & + \int_0^t \int_0^t h_2 (\tau_1, \tau_2) N(t-\tau_1) N(t-\tau_2) d\tau_1 d\tau_2　\\ \end{alignedat}

さらに、Friston et al. (1998) では、22 次ボルテラカーネルを次のように設計しました。

h2(τ1,τ2)=∑i=13∑j=13βijbi(τ1)bj(τ2)h_2 (\tau_1, \tau_2)=\sum_{i=1}^3 \sum_{j=1}^3 \beta_{ij}b_i(\tau_1)b_j(\tau_2)

ここで、b(τ)b(\tau) は次のようなガンマ基底関数です。

b1(τ)=13!τ3e−τb2(τ)=17!τ7e−τb3(τ)=115!τ15e−τ\begin{alignedat}{2} b_1(\tau) &= \frac{1}{3!} \tau^3 e^{-\tau} \\ b_2(\tau) &= \frac{1}{7!} \tau^7 e^{-\tau} \\ b_3(\tau) &= \frac{1}{15!} \tau^{15} e^{-\tau} \\ \end{alignedat}

したがって、BOLD反応の 22 次ボルテラモデルは次のようになります。

y(t)= ∑i=13αi[bi(t)∗N(t)]+∑i=13∑j=13βij[bi(t)∗N(t)][bj(t)∗N(t)]= ∑i=13αixi(t)+∑i=13∑j=13βijxi(t)xj(t)\begin{alignedat}{2} y(t) = \ & \sum_{i=1}^3 \alpha_i [b_i(t) \ast N(t)] \\ & + \sum_{i=1}^3 \sum_{j=1}^3 \beta_{ij} [b_i(t) \ast N(t)] [b_j(t) \ast N(t)] \\ = \ & \sum_{i=1}^3 \alpha_i x_i(t) + \sum_{i=1}^3 \sum_{j=1}^3 \beta_{ij} x_i(t) x_j(t) \end{alignedat}

シミュレーション

ボルテラモデルを Python を使用してシミュレーションします。

次の例では、α1=α2=α3=0.2, β11=0.1, β13=−0.2\alpha_1 = \alpha_2 = \alpha_3 = 0.2, \ \beta_{11}=0.1, \ \beta_{13} = -0.2 とし、他のパラメーターを 00 としました。

まずは、１回の神経活動に対してシミュレーションしましょう。

Python 
import numpy as np
import math
import matplotlib.pyplot as plt

# Set the experimental condition
time = np.arange(100)
boxcar = np.zeros(100); boxcar[1:21] = np.ones(20)

# Calculate a Volterra series.
def volterra_bold(time, boxcar, alpha = [0.2, 0.2, 0.2], beta = [[0.1, 0, -0.2], [0, 0, 0], [0, 0, 0]]):
    b1 = ((np.power(time, 3)) * np.exp(-time) ) / math.factorial(3)
    b2 = ((np.power(time, 7)) * np.exp(-time) ) / math.factorial(7)
    b3 = ((np.power(time, 15)) * np.exp(-time) ) / math.factorial(15)
    x1 = np.convolve(boxcar, b1)
    x2 = np.convolve(boxcar, b2)
    x3 = np.convolve(boxcar, b3)
    x = np.array([x1, x2, x3])

    alpha = np.array(alpha)
    H1 = np.dot(alpha, x)

    beta = np.array(beta)
    H2 = np.dot(beta.flatten(),
                (np.array([x,x,x]) * np.array([x,x,x])).reshape(9,-1))
    return H1 + H2

# Plot the results.
time_for_display = [0,1] + list(range(1,16)) + list(range(15,100))
boxcar_for_display = [0]*2 + [1]*15 + [0]*85
f, ax = plt.subplots(1,1,dpi=150)
ax.plot(time_for_display, boxcar_for_display,
        label='Neural activation', alpha=0.5)
ax.plot(volterra_bold(time, boxcar),
        label='BOLD response')
ax.set_xlim(0,50); ax.set_ylim(-0.2,1.2)
ax.set_xlabel('time'); ax.set_ylabel('value')
ax.legend(loc='best', frameon=False)
plt.savefig('volterra_hrf1.png')
plt.show()

volterra_hrf1.png

次に、複数回の神経活動に対してもシミュレーションしましょう。

ここでは、全ての β=0\beta=0 とした線形モデルも一緒にシミュレーションしました。

Python 
# Set the experimental condition
time = np.arange(100)
boxcar = np.zeros(100); boxcar[1:16] = np.ones(15); boxcar[21:36] = np.ones(15); 

# Plot the results.
time_for_display = [0,1] + list(range(1,16)) + list(range(15,21)) + list(range(20, 35)) + list(range(34, 100))
boxcar_for_display = [0]*2 + [1]*15 + [0]*6 + [1]*15 + [0]*66
f, ax = plt.subplots(1,1,dpi=150)
ax.plot(time_for_display, boxcar_for_display,
        label='Neural activation', alpha=0.5)
ax.plot(volterra_bold(time, boxcar),
        label='BOLD (nonlinear)')
ax.plot(volterra_bold(time, boxcar, beta=[[0]*3,[0]*3,[0]*3]),
        label='BOLD (linear)')
ax.set_xlim(0,50); ax.set_ylim(-0.2,1.2)
ax.set_xlabel('time'); ax.set_ylabel('value')
ax.legend(loc='upper right', frameon=False)
plt.savefig('volterra_hrf2.png')
plt.show()

volterra_hrf2.png

Reference
  • Ashby, F. G. (2019). Statistical analysis of fMRI data. MIT press.
  • Buxton, R. B., Wong, E. C., & Frank, L. R. (1998). Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magnetic resonance in medicine, 39(6), 855-864.
  • Friston, K. J., Josephs, O., Rees, G., & Turner, R. (1998). Nonlinear event‐related responses in fMRI. Magnetic resonance in medicine, 39(1), 41-52.

Topics
My Twitter