Juliaを使ったニューロンの発火時系列の生成と推定

脳科学若手の会という団体のイベントでHonda Research Instituteの島崎秀昭先生に点過程によるニューロンの発火時系列の生成と推定に関するワークショップを依頼した。

ニューロンは発火活動という膜電位の急激な変化を通して互いに情報を伝達しあっている。その発火系列は各時刻ごとに発火or非発火のバイナリで表される時系列データのため、点過程という統計的手法を用いてモデリングを行うことができる。ワークショップではニューロンの発火時系列の解析の基礎となる点過程について島崎先生からご講演いただき、プログラミング言語Juliaを用いて参加者が各自コーディングに取り組んだ。

自分は運営スタッフとしてワークショップのコーディング部分のチューターを担当し、サンプルコードの作成とコーディングの補助を行った。ここではその内容を備忘録として簡単に記事にまとめる。内容としては発火頻度が時間変化する非定常ポアソン過程と、発火後の不応期やバーストを考慮するリニューアルガンマ過程の2つについて、発火系列の生成、背後の発火頻度の推定、推定モデルのバリデーションを実装した。

ポアソン過程

ニューロンの発火時刻Xがある時刻tより大きくなる確率をP(X > t)で表す。このとき、以下の仮定

P(X > s + t \mid X > t) = P(X > s),

を満たす過程をmemoryless processと呼ぶ。発火が起きる確率は1つ直前の発火時刻のみに依存するという意味になる。この仮定を満たす連続変数Xの確率分布は以下の密度関数

f(x) = \lambda e^{- \lambda x},

を持つ指数分布のみであることが知られている。なお $\lambda$ は指数分布のパラメーターで、ニューロンの単位時刻あたりの発火頻度に対応する。このようにmemorylessな特徴を持つ確率過程を(定常)ポアソン過程と呼ぶ。

ポアソン過程に従う乱数の生成には逆関数法が広く知られている。指数分布の累積密度関数は

F_X(X) = \int_{-\infty}^X f_X(s) ds = 1 - e^{-\lambda x},

で与えられ、このように変換された $U = F_X(X)$ は累積密度関数の定義より一様分布になる。一様分布の乱数はコンピューター上で簡単に生成できるため、累積密度関数の逆関数を用いて

\begin{aligned} X &= F^{-1}(U) \\ &= - \lambda^{-1} \log U, \end{aligned}

とすることでポアソン過程に従う乱数を生成することができる。 Juliaのコードで表現すると以下のようになる。以下のpprocess_poisson関数は発火頻度パラメーターlambdaを取り、発火時系列xsを出力として返す。

function pprocess_poisson(lambda)
    xs = [0.0];

    R = 0;
    
    X = 0.0;
    while x <= L*dt
        x = x + - log(rand()) / lambda;
        push!(xs, x);
    end
    
    return xs
end

Time-rescaling theorem

定常ポアソン過程では計測した全ての時刻で発火頻度が一定であることを仮定していたが、実際の実験データの解析を考えるとむしろ発火活動の時間変化を推定したい。ではそのような発火頻度が時々刻々変化する発火系列はどのように生成すればよいか？ナイーブなアプローチとしては時間変化のモデル化ごとに別々にアルゴリズムを構築することが挙げられるが、より怠惰な方向を探したい。そんな中で、任意の点過程を同一のアルゴリズムで取り扱うためにtime-rescaling theoremという理論がある。

Time-rescaling theoremは直感的には逆関数法の点過程バージョンである。本理論によれば、任意の点過程は単一頻度のポアソン過程に変換することができる。この関係を逆に利用することで、前述した単一頻度のポアソン過程に従う乱数を生成し逆変換することで任意の点過程に従う乱数が得られる。

任意の点過程について、ISIの密度関数 $f(\tau)$ と条件付き強度関数 $\lambda(\tau)$ の間には以下の1対1の関係が成り立つ。

\begin{aligned} \lambda (\tau) &= \frac{ f(\tau) }{ 1 - \int_0^{\tau} f(u) du }, \\ f(\tau) &= \lambda(\tau) \exp \left[ - \int_0^{\tau} \lambda(u) du \right]. \end{aligned}

なお、ここで $\tau = t - t^*$ であり、1つ前の発火からの経過時刻を表す。この関係式を用いると、以下の擬似コードで示す生成アルゴリズムを一般の点過程に用いることができる

for $i$ in $1...\infty$

一様乱数を用いて指数分布に従う乱数を生成する: $\zeta \sim - \log \xi$ .
条件付き強度関数が $\zeta$ に達するまで積分する: $\zeta = \int_{t_{i-1}}^{\eta} \lambda(u \mid H_u) du$ .
$i$ 番目の発火時刻を $t_i = \eta$ とする。
$\eta \ge T$ なら終了。

なお、ここで $H_u$ は時刻 $u$ までの発火履歴であり、 $T$ は総計測時間である。つまり $\lambda(u \mid H_u)$ という表記は $\lambda(\tau)$ のより一般的な形式と見ることができる。ここからは具体的に非定常ポアソン過程とリニューアルガンマ過程に対して本アルゴリズムを適用した結果を紹介する。

非定常ポアソン過程

非定常ポアソン過程は、まさにポアソン過程の確率密度関数を以下の式

f(t \mid t_{i-1}) = \lambda(t) \exp \left[ - \int_{t_{i-1}}^t \lambda(u) du \right],

のように発火頻度パラメーターλを時間依存のλ(t)に拡張した形になっている。この過程に前述のアルゴリズムを適用すると入力としてλ(t)を直接受け取った際の発火系列の生成となるので、以下のようなコードで実現できる。

function pprocess_inhomopoisson(lambda)
    xs = [0.0];
    rs = zeros(L);
    R = 0;
    rs[1] = lambda[1];

    eta = - log(rand());
    for i = 2:L
        tau = i*dt - x[end];

        r = lambda[i];
        R = R + r*dt;

        if eta <= R
            push!(xs,i*dt);

            R = 0.0;
            eta = - log(rand());
        end
        rs[i] = r;
    end

    return xs, rs
end

ここでrは各時刻の発火頻度λ(t)、Rはその積分値を意味する。以下にその生成結果を示す。下図に青い線で示すような時間変化する発火頻度に対して、赤い点で示すような発火系列が得られた。

リニューアルガンマ過程

リニューアル過程は、発火が起きることでそのニューロンの発火頻度が一時的に上下するような挙動をモデル化する過程である。パラメーターによって例えば発火不応期であったり、逆にバースト発火するようなニューロンの特徴を捉えることができる。計算の容易さと実データへのフィッティングの観点から、典型的にはガンマ過程が用いられる。

ガンマ過程はISIの密度関数が以下の式

f(\tau) = \frac{1}{\Gamma(\kappa)} (\kappa \lambda)^\kappa \tau^{\kappa-1} e^{-\kappa\lambda\tau},

に従う過程である。 $\lambda,\kappa$ はそれぞれベースラインの発火頻度を決めるパラメーターと形状パラメーターである。 $\kappa$ の値の取り方によってニューロンの性質としてregular spikingかfast spikingかをコントロールすることができる。この過程に前述のアルゴリズムを適用する場合は入力としてこれらのパラメーター $\lambda,\kappa$ を取りISIの密度関数を計算し、 $\lambda(\tau)=\frac{f(\tau)}{1-F(\tau)}$ の関係を用いて逐次的に条件付き強度関数を更新する。これは以下のようなコードで実現できる。

function pprocess_gamma(lambda,k)
    x = [0.0];
    rs = zeros(L);

    P = 0; R = 0;

    rs[1] = lambda;

    eta = - log(rand());
    for i = 2:L
        tau = i*dt - x[end]

        p = 1.0 / gamma(k) * (k * lambda).^k * tau.^(k-1) .* exp(-k*lambda*tau);
        P = P + p*dt;

        r = p / (1.0 - P);
        R = R + r*dt;

        if eta <= R
            push!(x,i*dt)

            P = 0.0; R = 0.0;
            eta = - log(rand());
        end
        rs[i] = r;
    end

    return x, rs
end

ここでpはISIの密度関数 $f(\tau)$ 、rは各時刻の発火頻度 $\lambda(t)$ 、P, Rはその積分値を意味する。非定常ポアソン過程のときと同様に発火系列の生成を行うと、以下のような図が得られた。

ここまでは点過程に従う乱数の生成について紹介してきた。次に、得られた発火系列から発火頻度の時間変化を推定することを考える。

ガウスカーネル回帰を用いた非定常点過程の推定

点列からその背後にある関数を推定する一番素直な方法はカーネル回帰だろう。発火時刻の系列

x(t) = \sum_{i=1}^{n} \delta (t - t_i),

が得られた時、背後の発火頻度の時間変化をカーネル関数 $k_w$ を用いて

\lambda(t) = \sum_{i=1}^n k_w (t - t_i) = \int_0^T k_w (t - s) x(s) ds,

でフィッティングすることを考える。カーネル関数には典型的にはガウス関数 $k_w(s) = \frac{1}{\sqrt{2\pi}w} e^{-\frac{s}{2w^2}}$ が用いられる。カーネル回帰は発火系列の生成過程について過程を置かないノンパラメトリックな手法になるので、生成モデルが未知の場合には有効な手法である。この時、 $w$ はカーネルのバンド幅でコントロールパラメーターとなる。真の発火頻度との二乗誤差

\mathrm{MISE} = \mathbb{E} \int \left( \lambda_t - \hat{\lambda_t} \right)^2 dt,

を最小化することを考えると、バンド幅の目的関数として

C(w) = \sum_{i,j} \int k_w(s-t_i) k_w(s-t_j) ds - 2 \sum_{i \ne j} k_w(t_i, t_j),

が得られることが知られている [Shimazaki & Shinomoto (2010)]。ここで紹介したバンド幅の最適化を含めた発火頻度の推定は以下のコードで実現できる。

unction compute_optimal_bw(C, ws)
    cs = C.(ws);
    cmin, idx = findmin(cs);
    wopt = ws[idx];
    wopt
end

function get_costfunc(data)
    C(w) = sum(
        [
            sum(gauss(x-data,sqrt(2)*w)) - 2.0*sum(gauss(x-data,w)) + 2.0*gauss(0,w) for x in data
        ]
    );
    C
end

function sskernel(data; ws=nothing)
    if typeof(ws) == Void
        ws = linspace(minimum(diff(data)) / 5.0, maximum(data), 300)
    end

    C = get_costfunc(data);
    w = compute_optimal_bw(C, ws);

    density(ts) = [mean(gauss(t - data, w)) for t in ts];
    rs = density(ts) * length(data);
    rs
end

先程非定常ポアソン過程で生成した発火系列から発火頻度を推定すると以下の図のようになる。真の発火頻度(青)に対して、発火頻度の時間変化を精度よく推定できていると言える(橙)。

リニューアルガンマ過程に対する最尤推定

発火活動の生成モデルにある程度あたりがついている場合はそれを用いた方がよい。今回のケースであるリニューアルガンマ過程を生成モデルとして仮定した場合の最尤推定について紹介する。今回は推定アルゴリズムとして、自然パラメーターの空間での勾配法を実装する。

ガンマ過程の密度関数

f(\tau) = \frac{1}{\Gamma(\kappa)} (\kappa \lambda)^\kappa \tau^{\kappa-1} e^{-\kappa\lambda\tau},

は自然パラメーター $\alpha,\beta$ を使って

\begin{aligned} \alpha &= \kappa - 1, \\ \beta &= \kappa \lambda, \end{aligned}

と変数変換すると

f(\tau) = \frac{1}{\Gamma(\alpha+1)} \exp \left[ \alpha \log \tau - \beta \tau + (\alpha+1) \log \beta - \log \Gamma(\alpha+1) \right],

の様に書き換えることができる。この時の対数尤度は

\begin{aligned} \mathcal{L}(\alpha, \beta) &= \log f(\tau_1, \tau_2, ..., \tau_n) \\ &= \alpha \sum_i \log \tau_i - \beta \sum_i \tau_i + n \left\{ (\alpha+1) \log \beta - \log \alpha - \log \Gamma(\alpha) \right\}, \end{aligned}

と書ける。それぞれのパラメーターについて勾配を計算すれば、

\begin{aligned} \frac{ \partial \mathcal{L}(\alpha, \beta) }{ \partial \alpha } &= \sum_i \log \tau_i + n \left\{ \log \beta - \frac{1}{\alpha} - \phi(\alpha) \right\}, \\ \frac{ \partial \mathcal{L}(\alpha, \beta) }{ \partial \beta } &= - \sum_i \tau_i + n \frac{\alpha+1}{\beta}, \end{aligned}

となるので、これを用いて対数尤度を最大化する。実装は以下のようになる。

function fitgamma(isi)

    epsilon = 1e-6;

    n = length(isi);
    T = sum(isi);

    lambda = n / T;

    alpha = 0.2
    beta = 1.0

    eta = 0.001

    dll_dalpha = epsilon * 10
    dll_dbeta = epsilon * 10

    while max(abs(dll_dalpha), abs(dll_dbeta)) > epsilon

        dll_dalpha = - n / alpha + n * log(beta) + sum(log(isi)) - n * digamma(alpha);
        dll_dbeta = n * (1+alpha) / beta - T;

        alpha = alpha + eta*dll_dalpha;
        beta = beta + eta*dll_dbeta;

    end

    k = alpha + 1
    lambda = beta / k

    lambda, k

end

alpha, betaの初期値や学習率eta、収束判定あたりはだいぶ雑に作ってあるので改良の余地はあるかもしれない。非定常ポアソン過程のときと同様に、生成された発火系列から発火頻度を推定すると以下の図のようになる。こちらについても適切に真の発火頻度を推定できているように見える。

推定モデルのバリデーション

最後に、ガウスカーネル回帰やガンマ過程の最尤推定など、推定のモデルによらずにノンパラメトリックに推定モデルの評価を行う手法を紹介する。

発火系列推定の際に用いたtime-rescaling theoremは、そのままモデルの評価にも使用することができる。 Time-rescaling theoremでは任意の点過程を単一頻度のポアソン過程に変換することができる。この性質を利用して、任意の推定モデルで推定された発火頻度を一様分布の空間に写像し、単一頻度のポアソン過程を同様に写像したものと比較することで、推定結果の当てはまりのよさを定量化できる。これはQuantile-Quantile plot (Q-Q plot)として以下のようなコードで実現できる。

function QQplot(rs, data)
    R=cumsum(rs*dt);

    x_step = data/dt;
    x_round = Int.(ceil(transpose(x_step))+1);

    R_x = R[x_round][:];
    tau = diff(R_x);

    u_hat = 1-exp(-tau);
    u_hat_sorted = sort(u_hat);

    n = length(u_hat);
    u = ([1:n]-0.5)/n;

    cz =linspace(0,1,100);
    confb95 = 1.96*sqrt(cz.*(1-cz)/n);
    confb99 = 2.575*sqrt(cz.*(1-cz)/n);
    p_x = [0,1]; p_y = [0,1];
    plot(p_x, p_y, color="black", linewidth=2.0, linestyle=:dash, label="")
    plot!(cz, cz-confb95, color="black", linewidth=2.0, linestyle=:solid, label="")
    plot!(cz, cz+confb95, color="black", linewidth=2.0, linestyle=:solid, label="")

    plot!(u, u_hat_sorted, linewidth=0, marker="", label="")
end

実験として、非定常ポアソン過程で生成した発火系列に対して、ガウスカーネル回帰とリニューアルガンマ過程を仮定した最尤推定の2つの推定を行い、それぞれの結果をQ-Q plotで比較した。生成過程は非定常ポアソン過程なので、後者の最尤推定は悪い推定となっているはずである。

上の図は、ガウスカーネル回帰を行った際のQ-Q plotである。 x軸にデータ、y軸にモデルの推定した結果をプロットしており、ほぼ $y=x$ の軸に乗っていることから、よい推定であることがわかる。

一方で、上の図はデータの生成過程とは異なる、リニューアルガンマ過程の最尤推定を行った結果である。グラフの右部分で点列がy = xを大きく外れており、推定モデルで説明できない効果が見られたことが示唆される。

このように、Q-Q plotを用いることで、真のモデルがわからない場合にも推定モデルの当てはまりのよさを評価することができる。

Putting it all together

今回紹介したコードの詳細は以下のリポジトリで確認できる。

GitHub - ganow/julia-point-process: The generation/estimation/validation of spike sequences with point process.

The generation/estimation/validation of spike sequences with point process. All codes were written in Julia. These codes were used in a workshop on the point process of Student Association for Brain Science's training camp. We made them by referring to Dr. Shimazaki's implementation in Matlab. Detailed explanation is summarized in ganow's blog (written in Japanese only).

github.com

GitHub - ganow/julia-point-process: The generation/estimation/validation of spike sequences with point process.