指数関数型シナプスモデル

シナプスのモデルは複数あるが, 良く用いられるのが指数関数型シナプスモデル(exponential synapse model)である．このモデルは生理学的な過程を無視した現象論的モデルであることに注意しよう．指数関数型シナプスモデルには2つの種類, 単一指数関数型モデル (single exponential model)と二重指数関数型モデル (double exponential model)がある．

数式の説明の前にモデルの挙動を示す．次図は2種類のモデルにおいて$t=0$でスパイクが生じてからのシナプス後電流の変化を示している．ただし, 実際のシナプス後電流はこれにシナプス強度 (Synaptic strength)1を乗じて総和を取ったものとなる．

using PyPlot
rc("axes.spines", top=false, right=false)

td, tr = 2e-2, 2e-3 # synaptic decay time, synaptic rise time (sec)
dt, T = 5e-5, 0.1 # タイムステップ, シミュレーション時間 (sec)
nt = Int(T/dt) # シミュレーションの総ステップ

# 単一指数関数型シナプス
r_single = zeros(nt)
for t in 1:nt-1
    spike = ifelse(t == 1, 1, 0)
    r_single[t+1] = r_single[t]*(1-dt/td) + spike/td
    #r_single[t+1] = r_single[t]*exp(-dt/td) + spike/td
end

# 二重指数関数型シナプス
r_double, hr = zeros(nt), zeros(nt)
for t in 1:nt-1
    spike = ifelse(t == 1, 1, 0)
    r_double[t+1] = r_double[t]*(1-dt/tr) + hr[t]*dt
    hr[t+1] = hr[t]*(1-dt/td) + spike/(tr*td)
    #r_double[t+1] = r_double[t]*exp(-dt/tr) + hr[t]*dt
    #hr[t+1] = hr[t]*exp(-dt/td) + spike/(tr*td)
end   

time = (1:nt)*dt
figure(figsize=(4, 3))
plot(time, r_single, linestyle="dashed", label="single exponential")
plot(time, r_double, label="double exponential")
xlabel("Time (s)"); ylabel("Post-synaptic current (pA)")
legend(); tight_layout()

2種類の指数関数型シナプスの動態．破線は単一指数関数型シナプスで, 実線は二重指数関数型シナプスである．

単一指数関数型モデル(Single exponential model)

シナプス前ニューロンにおいてスパイクが生じてからのシナプス後電流の変化はおおよそ指数関数的に減少する, というのが単一指数関数型モデルである 2. 式は次のようになる．

$$f(t)=\frac{1}{{\tau }_s}\exp (-\frac{t}{{\tau }_s})$$

この関数を時間的なフィルターとして, 過去の全てのスパイクについての総和を取る．

$$r(t)=\sum _{t_k<t}f(t-t_k)$$

ここで$r(t)$は前節におけるシナプス動態($s_{\text{syn}}$)で, $t_k$はあるニューロンの$k$番目のスパイクの発生時刻である．$t_k<t$の意味は現在の時刻$t$までに発生したスパイクについての和を取るという意味である．なお，スパイクが生じてから, ある程度の時間が経過した後はそのスパイクの影響はないと見なせるので, 一定の時間までの総和を取るのがよい．

別の表記法としてスパイク列に対する畳み込みを行うというものもある．畳み込み演算子を$\ast$とし, シナプス前細胞のスパイク列を$S(t)=\sum _{t_k<t}\delta (t-t_k)$とする (ただし, $\delta$はDiracのdelta関数において$\delta (0)=1$とした関数)．このとき, $r(t)=f\ast S(t)$と表すことができる．畳み込み演算子を用いると簡略な表記ができるが，実装上は他と同じ手法を用いる．

微分方程式による表現

上の手法ではニューロンの発火時刻を記憶し, 時間毎に全てのスパイクについての和を取る必要がある．そこで, 実装する場合は次の等価な微分方程式を用いる．

$$\frac{dr}{dt}=-\frac{r}{{\tau }_s}+\frac{1}{{\tau }_s}\sum _{t_k<t}\delta (t-t_k)$$

ここで${\tau }_s$はシナプスの時定数(synaptic time constant)である． また, $\delta (\cdot )$はDiracのdelta関数です(ただし$\delta (0)=1$です). これをEuler法で差分化すると

$$r(t+\mathrm{\Delta }t)=(1-\frac{\mathrm{\Delta }t}{{\tau }_s})r(t)+\frac{1}{{\tau }_s}{\delta }_{t,t_k}$$

となる．ここで${\delta }_{t,t_k}$はKroneckerのdelta関数で, $t=t_k$のときに1, それ以外は0となる．また減衰度として$(1-\mathrm{\Delta }t/{\tau }_d)$の代わりに$\exp (-\mathrm{\Delta }t/{\tau }_d)$を用いる場合もある．

二重指数関数型モデル(Double exponential model)

2重の指数関数によりシナプス後電流の立ち上がりも考慮するのが, 二重指数関数型モデル(Double exponential model)である3．$t=0$にシナプス前細胞が発火したときのシナプス後電流の時間変化の関数は次のようになる．

$$f(t)=A[\exp (-\frac{t}{{\tau }_d})-\exp (-\frac{t}{{\tau }_r})]$$

ただし, ${\tau }_r$は立ち上がり時定数(synaptic rise time constant), ${\tau }_d$は減衰時定数(synaptic decay time constant)である．${\tau }_d$は${\tau }_s$と同じく神経伝達物質の減少速度を決定している．$A$は規格化定数で次のように表される．

$$A=\frac{{\tau }_d}{{\tau }_d-{\tau }_r}\cdot {\left(\frac{{\tau }_r}{{\tau }_d}\right)}^{\frac{{\tau }_r}{{\tau }_r-{\tau }_d}}$$

規格化定数$A$を乗じることで最大値が1となる．ただし, シミュレーションをする上で実際に規格化をする場合は少ない．

$\alpha$関数

上記の式において, $\tau ={\tau }_r={\tau }_d$の場合は $\mathit{\alpha }$ 関数 (alpha function, alpha synapse)と呼ぶ (Rall, 1967)．式としては次のようになる．

$$\alpha (t)=\frac{t}{\tau }\exp (1-\frac{t}{\tau })$$

この式は二重指数関数型シナプスの式に単に代入するだけでは導出できない．これらの式の対応については後述する．

微分方程式による表現

ここで, 二重指数関数型シナプスの式に対応する, 補助変数$h$を用いた微分方程式を導入する．

$$\begin{array}{r}\begin{array}{rl}\frac{dr}{dt}& =-\frac{r}{{\tau }_d}+h\\ \frac{dh}{dt}& =-\frac{h}{{\tau }_r}+\frac{1}{{\tau }_r{\tau }_d}\sum _{t_k<t}\delta (t-t_k)\end{array}\end{array}$$

単一指数関数型シナプスの場合と同様にEuler法で差分化すると

$$\begin{array}{r}\begin{array}{rl}r(t+\mathrm{\Delta }t)& =(1-\frac{\mathrm{\Delta }t}{{\tau }_d})r(t)+h(t)\cdot \mathrm{\Delta }t\\ h(t+\mathrm{\Delta }t)& =(1-\frac{\mathrm{\Delta }t}{{\tau }_r})h(t)+\frac{1}{{\tau }_r{\tau }_d}{\delta }_{t,t_{jk}}\end{array}\end{array}$$

となる．

念のため, 微分方程式と元の式が一致することを確認しておこう．$t=0$のときにシナプス前細胞が発火したとし, それ以降の発火はないとする．このとき, $h(0)=1/{\tau }_r{\tau }_d,r(0)=0$ である．$h$についての微分方程式の解は

$$h(t)=h(0)\cdot \exp (-\frac{t}{{\tau }_r})$$

となるので, これを$r$についての式に代入して

$$\frac{dr}{dt}=-\frac{r}{{\tau }_d}+h(0)\cdot \exp (-\frac{t}{{\tau }_r})$$

となる．これを解くには両辺に積分因子$\exp (t/{\tau }_d)$をかけてから積分をするかLaplace変換をするかである．今回はLaplace変換を用いる．右辺一項目を移行した後に両辺をLaplace変換すると以下のようになる．

$$\begin{array}{r}\begin{array}{rl}\mathcal{L}[\frac{dr}{dt}+r/{\tau }_d]& =\mathcal{L}[h(0)\cdot \exp (-t/{\tau }_r)]\\ sF(s)-r(0)+\frac{1}{{\tau }_d}F(s)& =\frac{h(0)}{s+1/{\tau }_r}\\ F(s)& =\frac{h(0)}{(s+1/{\tau }_r)(s+1/{\tau }_d)}\end{array}\end{array}$$

ただし$r(t)$のLaplace変換を$F(s)$とした. ここで逆Laplace変換を行うと次のようになる．

$$\begin{array}{r}\begin{array}{rl}r(t)& ={\mathcal{L}}^{-1}(F(s))\\ & ={\mathcal{L}}^{-1}\left[\frac{h(0)}{(s+1/{\tau }_r)(s+1/{\tau }_d)}\right]\\ & ={\mathcal{L}}^{-1}\left[\frac{h(0)}{1/{\tau }_r-1/{\tau }_d}(\frac{1}{s+1/{\tau }_d}-\frac{1}{s+1/{\tau }_r})\right]\\ & =\frac{1}{{\tau }_d-{\tau }_r}[\exp (-t/{\tau }_d)-\exp (-t/{\tau }_r)]\end{array}\end{array}$$

この式の最大値$r_{max}$を求めておこう． $r(t)$を微分して0と置いた式の解$t_{max}$を代入すれば求められる．計算すると,

$$t_{max}={\displaystyle \frac{\ln ({\tau }_d/{\tau }_r)}{1/{\tau }_r-1/{\tau }_d}},r_{max}={\displaystyle \frac{1}{{\tau }_d}}\cdot {\left({\displaystyle \frac{{\tau }_r}{{\tau }_d}}\right)}^{\frac{{\tau }_r}{{\tau }_d-{\tau }_r}}$$

となる．

なお, $\alpha$関数の導出は逆Laplace変換をする前に$\tau ={\tau }_d={\tau }_r$とすればよく,

$$\begin{array}{r}\begin{array}{rl}{F}_{\alpha }(s)& =\frac{h(0)}{(s+1/\tau {)}^2}\\ \alpha (t)& =\frac{t}{{\tau }^2}\exp (-\frac{t}{\tau })\end{array}\end{array}$$

となる．若干の係数の違いはあるが, 同じ形の関数が導出された．

---

1

シナプス強度というのは便宜上の呼称で, 実際には神経伝達物質の種類や, その受容体の数など複数の要因によって決定されている. また, このシナプス強度はシナプス重みということもある．これはどちらかと言えば機械学習の表現に引っ張られたものである．このため, このサイトでは重みという語も使う．

2

薬学動態の静注1コンパートメントモデルと同じ式である．

3

薬学動態の内服1コンパートメントモデルと同じ式である．