Julia で眼優位性マップ形成のシミュレーションをしてみる

NOTE

本記事で使用しているコードは, 「はじめての神経回路シミュレーション」¹における C コードを Julia で書き換えたものです.

大脳皮質の構造と機能

大脳皮質は, 表面の深い溝により前頭葉, 頭頂葉, 後頭葉, 側頭葉に分けられ, 各脳葉はそれぞれ異なる機能を持つ. Korbinian Brodmann は, 大脳皮質を解剖学的構造の違いに基づいて $52$ 個の領域に分割した² (図1). これはブロードマンの脳地図とよばれ, 概ね機能の違いに対応した構造をしている.

図1: ブロードマンの脳地図

また, 大脳皮質は6層構造³ (図2) である.

図2: 大脳皮質の6層構造

表面に最も近い I 層は II/III 層のニューロンの樹状突起と, 大脳皮質の別の領域からやってくる軸索がある. II/III 層は主に小型の錐体ニューロンからなっており, 大脳皮質の別の領域にあるニューロンから信号を受け取り, 大脳皮質間の水平結合を構成している. IV 層は小型の球状ニューロンからなっており, 視床からの感覚入力が主要な入力信号である. V/VI 層は大型の錐体ニューロンからなっており, 他の大脳皮質領域や皮質下構造へと信号を送る.

大脳皮質の学習

大脳皮質では, 教師あり学習における教師信号や強化学習における報酬予測信号等の入力がはっきりとしていないため, 大脳皮質で行われる学習は教師なし学習であると考えられている. これは, 感覚情報の学習において顕著であり, 特に第1次視覚野において膨大な研究がなされてきた.

大脳皮質の IV 層は身体からの感覚情報が最初に入ってくる場所であり, 特に視覚情報を最初に受け取る第1次視覚野 (V1) では厚い層を形成している. 以降では, V1 における情報処理について考える.

動物には右目と左目があり, それぞれの目から入ってきた視覚情報は視床を介して V1 の IV 層に入力し, 大脳皮質での最初の情報処理が行われる. IV 層のニューロンは右目・左目のどちらか一方の入力に対して選択的に応答する性質を有しており, これを眼優位性とよぶ. また, 皮質上で近くにいるニューロン同士は同じ眼優位性を有しており, 右目優位・左目優位の領域が皮質の表面上で縞模様もしくは斑点模様を形成していることが知られている. この空間的な構造を, 眼優位性マップとよぶ.

眼優位性マップ形成のシミュレーション

モデル

左目・右目の網膜内の光に応答する細胞である視細胞に対応するスパイキングニューロンを1つずつ用意し, それぞれは独立にスパイクを発射するものとする. なお, 本シミュレーションでは $100\, \mathrm{spikes/s}$ のポアソンスパイクとし, $\Delta t\, (=1\, \mathrm{ms})$ 毎に乱数を振ってスパイクを生成し, 以下の変数に格納する：

\mathrm{s\_eye}[X]=\left\{ \begin{array}{ll} 1 & \text{if a spike was fired at that time} \\ 0 & \text{otherwise} \end{array}\right.

ここで, $X$ は左目 $(L)$ もしくは右目 $(R)$ の視細胞のいずれかである.

図3: 眼優位性マップ形成モデル⁴

視覚野は2次元平面を考え, 本シミュレーションでは $32\times 32$ 個のニューロンを設定する. 視覚野のニューロン $i$ への入力は, 視細胞からのスパイク入力 $\mathrm{aff}(i)$ と, 皮質内のニューロン $i,j$ 間の水平結合を介した他のニューロンからのスパイク入力 $\mathrm{lat}(i)$ からなる. $\mathrm{aff}(i)$ は視細胞 $\mathrm{L,\, R}$ のスパイクと, 結合強度 $w_{i,\mathrm{L}},\, w_{i,\mathrm{R}}$ の積和であり, 次式で表される：

\begin{equation*} \mathrm{aff}(i)=w_{i,\mathrm{L}}\cdot \mathrm{s\_eye[L]}+w_{i,\mathrm{R}}\cdot\mathrm{s\_eye[R]} \end{equation*}

$\mathrm{lat}(i)$ は, 視覚野のニューロン $j$ のスパイクを変数 $\mathrm{s\_ctx}[j]$ でニューロン $i,j$ 間の水平結合の関数を $I(i,j)$ でそれぞれ定義し, 次式で計算する：

\begin{equation*} \mathrm{lat}(i)=\displaystyle\sum_{j}I(i,j)\cdot\mathrm{s\_ctx}[j] \end{equation*}

また, $I(i,j)$ は次式で表される：

\begin{equation*} I(i,j)=k_{1}\exp \left(\dfrac{-|i-j|^{2}}{2\sigma^{2}_{1}}\right)-k_{2}\exp\left(\dfrac{-|i-j|^{2}}{2\sigma^{2}_{2}}\right) \end{equation*}

ここで, $|i-j|$ は, 皮質の2次元平面上でのニューロン $i,j$ の距離とする.

この2つの値 $(\mathrm{aff}(i),\, \mathrm{lat}(i))$ から電流 $g_{\mathrm{aff}}$ , $g_{\mathrm{lat}}$ を計算し, 膜電位の更新, スパイクの発射を行う. シナプス結合の強度の初期値は, 左右から等しく結合しているとして, $0.5\pm 0.1$ の範囲でランダムに決める. ニューロンはカレントベースの LIF モデル, シナプスは指数減数モデルをそれぞれ採用する.

シナプス可塑性

ここでは, STDP を実装する. まず, $i,j$ をそれぞれポスト側, プレ側のニューロンの添え字とする. プレ側のトレース $x_{j}$ を以下のように定義する：

\begin{equation*} x_{j}(t)=\exp\left(-\dfrac{\Delta t}{\tau_{\mathrm{pre}}}\right)x_{j}(t-\Delta t)+S_{j}(t) \end{equation*}

ここで, $\tau_{\mathrm{pre}}$ は減衰の時定数である. また, $S_{j}(t)$ は以下の式で定義される：

S_{j}(t)=\left\{ \begin{array}{ll} 1 & \text{if a spike was fired at time }\, t \\ 0 & \text{otherwise} \end{array}\right.

ポスト側のトレース $y_{i}$ は以下のように定義する：

\begin{equation*} y_{i}(t)=\exp\left(-\dfrac{\Delta t}{\tau_{\mathrm{post}}}\right)y_{i}(t-\Delta t)+S_{i}(t) \end{equation*}

ここで, $\tau_{\mathrm{post}}$ は減衰の時定数である.

これらの値を用いて, 結合強度を次式で更新する：

\begin{align*} &w_{i,j}(t+\Delta t)=w_{i,j}(t)+\Delta w_{i,j}(t),\\ &\Delta w_{i,j}(t)=-A_{-}y_{i}(t)S_{j}(t)+A_{+}x_{j}(t)S_{i}(t) \end{align*}

ここで, $A_{-},A_{+}$ は定数である.

シミュレーション

実際にシミュレーションしてみる. まずは以下のように記述する.

using Plots
using Random
using LaTeXStrings
using Interpolations

plot_font = "Computer Modern"
default(fontfamily=plot_font,
    linewidth=1, label=nothing, grid=true)

T = 4000 # ms
dt = 1   # ms
num_t = Int(T / dt) + 1
t_list = [t * dt for t in 0:num_t-1]
tau_syn = 5.0 # ms
tau_STDP = 10.0 # ms

R = 1  # MΩ
V_eq = -49.0  # mV 
V_reset = -60 # mV
V_ex = 20.0   # mV
theta = -50.0 # mV
tau = 20 # ms

num_nx = 32
num_ny = 32
num_N = num_nx * num_ny
A_plus = 0.02
A_minus = 0.02

s_eye_L = zeros(num_t)
s_eye_R = zeros(num_t)
eye_L_STDP = zeros(num_t)
eye_R_STDP = zeros(num_t)
Random.seed!(123)
for t in 1:num_t
    if rand() <= 0.1
        s_eye_L[t] = 1
    end
    if rand() <= 0.01
        s_eye_R[t] = 1
    end
end
eye_L_STDP[1] = s_eye_L[1]
eye_R_STDP[1] = s_eye_R[1]
for t in 2:num_t
    eye_L_STDP[t] = eye_L_STDP[t-1] * exp(-dt / tau_STDP) + s_eye_L[t]
    eye_R_STDP[t] = eye_R_STDP[t-1] * exp(-dt / tau_STDP) + s_eye_R[t]
end

g_aff = 2.0
g_lat = 4.0

sig_1 = 4.0
sig_2 = 16.0
k1 = 1.0
k2 = sig_1 / sig_2
norm_coeff = 1.0 / (k1 * sqrt(2.0 * sig_1 * pi) - k2 * sqrt(2.0 * sig_2 * pi))

function Calc_I_ij(dif)
    return norm_coeff * (k1 * exp(-dif / (2.0 * sig_1)) + (-k2) * exp(-dif / (2.0 * sig_2)))
end

I_ij_list = zeros((num_nx, num_ny))

for i in 1:num_nx
    for j in 1:num_ny
        d2 = i * i + j * j
        I_ij_list[i, j] = Calc_I_ij(d2)
    end
end

S_list = zeros((num_nx, num_ny))

mutable struct Neuron
    i::Int
    j::Int
    V::Float64
    S::Int
    STDP::Float64 # x_i,y_i
    w_L::Float64 # ω_{i,L}
    w_R::Float64 # ω_{i,R}
    I_ext::Float64 # (aff+lat)
    I_aff::Float64 # aff
    I_lat::Float64 # lat
    dif_list::Array{Float64,2}
end

function Neuron(i, j)
    V = V_reset + 10.0 * rand()
    S = 0
    STDP = 0.0
    w_L = 0.49 + 0.01 * rand()
    w_R = 0.49 + 0.01 * rand()
    I_ext = 0.0
    I_aff = 0.0
    I_lat = 0.0
    dif_list = zeros((num_nx, num_ny))
    for k in 1:num_nx
        for l in 1:num_ny
            dif_list[k, l] = I_ij_list[abs(i - k)+1, abs(j - l)+1]
        end
    end
    return Neuron(i, j, V, S, STDP, w_L, w_R, I_ext, I_aff, I_lat, dif_list)
end

function Calc_i!(neuron, time)
    # aff
    neuron.I_aff *= exp(-dt / tau_syn)
    neuron.I_aff += neuron.w_L * s_eye_L[time] + neuron.w_R * s_eye_R[time]

    # lat
    neuron.I_lat *= exp(-dt / tau_syn)
    temp = S_list
    temp[neuron.i, neuron.j] = 0.0
    lat_list = neuron.dif_list .* temp
    neuron.I_lat += sum(lat_list)

    # ext
    neuron.I_ext = neuron.I_aff * g_aff + neuron.I_lat * g_lat
    return 0
end

function LIF!(neuron)
    neuron.S = 0
    temp = -(neuron.V - V_eq) + R * neuron.I_ext
    neuron.V += +dt * temp / tau
    if (neuron.V > theta)
        neuron.V = V_reset 
        neuron.S = 1
    end
    return 0
end

function Update_Trace!(neuron)
    neuron.STDP = neuron.STDP * exp(-dt / tau_STDP) + neuron.S
    return 0
end

function Update_Omega!(neuron, time)
    neuron.w_L += A_plus * eye_L_STDP[time] * neuron.S + (-A_minus) * neuron.STDP * s_eye_L[time]
    neuron.w_R += A_plus * eye_R_STDP[time] * neuron.S + (-A_minus) * neuron.STDP * s_eye_R[time]
    neuron.w_L = max(0.0, neuron.w_L) 
    neuron.w_R = max(0.0, neuron.w_R) 
    neuron.w_L = min(1.0, neuron.w_L) 
    neuron.w_R = min(1.0, neuron.w_R) 
    return 0
end

N_list = [Neuron(i, j) for i in 1:num_nx, j in 1:num_ny]
Data_list1 = zeros((num_t, num_nx, num_ny))

function interpolate_data(data, new_size)
    old_size = size(data)
    x = range(1, stop=old_size[1], length=old_size[1])
    y = range(1, stop=old_size[2], length=old_size[2])
    xi = range(1, stop=old_size[1], length=new_size[1])
    yi = range(1, stop=old_size[2], length=new_size[2])
    interp = interpolate((x, y), data, Gridded(Linear()))
    return interp(xi, yi)
end

以下を実行してシミュレーションを可視化してみる.

@gif for t in 1:4001
    for i in 1:num_nx
        for j in 1:num_ny
            S_list[i, j] = N_list[i, j].S
        end
    end
    for i in 1:num_nx
        for j in 1:num_ny
            Calc_i!(N_list[i, j], t)
        end
    end
    for i in 1:num_nx
        for j in 1:num_ny
            LIF!(N_list[i, j])
        end
    end
    for i in 1:num_nx
        for j in 1:num_ny
            Update_Trace!(N_list[i, j])
        end
    end
    for i in 1:num_nx
        for j in 1:num_ny
            Update_Omega!(N_list[i, j], t)
        end
    end

    if (t % 50 == 0)
        for i in 1:num_nx
            for j in 1:num_ny
                Data_list1[t, i, j] = N_list[i, j].w_L - N_list[i, j].w_R
            end
        end
        data = Data_list1[t, :, :]

        new_size = (num_nx * 10, num_ny * 10)
        interpolated_data = interpolate_data(data, new_size)
        
        heatmap(interpolated_data, c=:seismic, clim=(-0.5, 0.5), title="t=$t")
    end
end every 50

上記により, 以下が得られる.

NOTE

ファイル読み込みに少し時間がかかります. また, 可視化の都合上, 補完しています.

図4: 眼優位性マップ形成のシミュレーション

図4において, 負は右目優位, 正は左目優位を表している. 図4より, 右目優位の領域 (青, 負値) と左目優位の領域 (赤, 正値) が縞模様のように現れ, V1 で見られる眼優位性マップと似たパターンが形成されたことが確認できる.

Footnotes

山﨑匡, 五十嵐潤, "はじめての神経回路シミュレーション: 1ニューロンからヒト全脳モデルまで", 講談社, 2021, ISBN: 978-4-627-85631-8. ↩
Amunts, Katrin. "Brodmann Areas." In Encyclopedia of Evolutionary Psychological Science, edited by Todd K. Shackelford and Viviana A. Weekes-Shackelford, 1-4. Cham: Springer International Publishing, 2018. https://doi.org/10.1007/978-3-319-16999-6_3341-1. ↩
University of Minnesota, "Neocortical Layers," accessed August 23, 2024, https://vanat.ahc.umn.edu/brain18/pages/neocorticalLayers.html. ↩
Miller, K., Keller, J. B., & Stryker, M. (1989). Ocular dominance column development: Analysis and simulation. Science, 245(4918), 605-611. https://doi.org/10.1126/science.2762813 ↩