大阪大学医学部 Python会

Now is better than never.

非負値行列因子分解(NMF)を用いたグリッドパターンの創発

2022-04-26(Tue) - Posted by 山本 in 技術ブログ    tag:Python tag:Neuroscience

Contents

    内側嗅内皮質(MEC)にあるグリッド細胞 (grid cells) は六角形格子状の発火パターンにより自己位置等を符号化するのに貢献している．この発火パターンを生み出すモデルは多数あるが，その中でもシンプルなのが(Dordek et al., eLife, 2016)である．場所細胞(place cells) の発火パターンを非負主成分分析(nonnegative principal component analysis) で次元削減するとグリッド細胞のパターンが生まれる，という内容である．このモデルはfeedforwardであり，場所細胞→グリッド細胞という逆向きの投射を考えているので完璧ではないが，グリッドパターンは最速で出すことができる．(Sorscher et al., NeurIPS2019, 2019)では非負値行列因子分解 (NMF: nonnegative matrix factorization) で問題ないと書いてあったのでお手軽に実行できるコードを書こうと思い今回の記事を書いた．

    ちなみに今回の内容はJuliaで学ぶ計算論的神経科学に書く予定だが，あちらではライブラリを必要最低限なものにある程度縛っているのでどういう構成にするか悩み中である．

    場所細胞の発火パターン

    訓練データとなる場所細胞の発火パターンを人工的に作成する．まずライブラリを読み込む．

    In [1]:
    import matplotlib.pyplot as plt
import numpy as np
from sklearn.decomposition import PCA, NMF
from scipy.signal import correlate2d

    場所細胞の発火パターンはDifference of Gaussians (DoG) で近似する．DoGは大きさの異なる2つのガウス関数の差分を取った関数であり，画像に適応すればband-passフィルタとして機能する．

    網膜神経節細胞等の受容野のON中心OFF周辺型受容野のモデルとしても用いられる．受容野中央では活動が大きく，その周辺では活動が抑制される，という特性を持つ．

    2次元のガウス関数とDoG関数を実装する．

    In [2]:
    def gaussian2d(center, width, height, step, sigma, scale=1):
    pos = np.array(np.meshgrid(np.linspace(-width/2, width/2, step),
                   np.linspace(-height/2, height/2, step)))
    g = np.exp(-((pos[0]-center[0])**2 + (pos[1]-center[1])**2) / (2*scale*(sigma**2)))
    return g / np.sum(g)
    In [3]:
    def DoG(center, width=2.2, height=2.2, step=55, sigma=0.12, surround_scale=2, rescale=True):
    g1 = gaussian2d(center, width, height, step, sigma)
    g2 = gaussian2d(center, width, height, step, sigma, surround_scale)
    dog = g1 - g2
    if rescale:
        dog = (dog - dog.min()) / (dog.max() - dog.min()) # recale in [0,1]
    return dog

    モデルのパラメータを設定する．場所細胞の数は増やせば増やすほど後の結果は綺麗になる．

    In [4]:
    Np = 5000          # 場所細胞の数
sigma = 0.12       # 場所細胞のtuning curveの幅 [m]
surround_scale = 2 # DoGのsigma^2の比率
box_width = 2.2    # 箱の横幅 [m]
box_height = 2.2   # 箱の縦幅 [m]
step = 40          # 空間位置の離散化数

    先にガウス関数とDoG関数がどのような見た目になるか確認しよう．

    In [5]:
    c_eg = [0,0]
gau_eg = gaussian2d(c_eg, box_width, box_height, step, sigma)
dog_eg = DoG(c_eg, box_width, box_height, step, sigma, surround_scale, rescale=False)
    In [6]:
    plt.figure(figsize=(5,5))
plt.subplot(2,2,1)
plt.title("Gaussian")
plt.imshow(gau_eg, cmap='turbo', extent=(-box_width/2, box_width/2, -box_height/2, box_height/2))
plt.ylabel("y [m]")
plt.subplot(2,2,2)
plt.title("Difference of\n Gaussians (DoG)")
plt.imshow(dog_eg, cmap='turbo', extent=(-box_width/2, box_width/2, -box_height/2, box_height/2))
plt.subplot(2,2,3)
plt.plot(np.linspace(-box_width/2, box_width/2, step), gau_eg[step//2])
plt.xlabel("x [m]"); plt.ylabel("Activity (firing rate)")
plt.subplot(2,2,4)
plt.plot(np.linspace(-box_width/2, box_width/2, step), dog_eg[step//2])
plt.xlabel("x [m]")
plt.tight_layout()

    場所細胞の活動パターンを生み出す．それぞれの場所受容野の中心は一様分布で生み出し，scalingにより境界上にcenterがある場所細胞を抑制しておく．

    In [7]:
    np.random.seed(0)
centers = np.vstack([np.random.uniform(-box_width/2, box_width/2, Np),
                     np.random.uniform(-box_height/2, box_height/2, Np)]).T
centers *= 0.9
    In [8]:
    X = np.array([DoG(c, box_width, box_height, step, sigma, surround_scale) for c in centers])

    9個の場所受容野を確認する．

    In [9]:
    plt.figure(figsize=(4,4))
plt.suptitle("Activation of place cells", fontsize=14)
for i in range(9):
    plt.subplot(3,3,i+1)
    plt.imshow(X[i].reshape(step,step), cmap='turbo')
    plt.axis('off')
plt.tight_layout()

    主成分分析 (PCA: principal component analysis)

    PCAとNMFの違いは，圧縮表現および変換のための重み行列の符号である．PCAは正負の両者を含むが，NMFでは圧縮表現と重みの両者が非負である．この意味は後で触れるとして，比較用にPCAを先に実行する．sklearn.decompositionのPCAを用いる．空間方向にflattenしてPCAを実行する．ここではグリッド細胞の数を9個に設定した．

    In [10]:
    pca = PCA(n_components=9, random_state=0).fit(X.reshape(Np, -1))

    圧縮表現を取得し，reshapeする．

    In [11]:
    H_pca = pca.components_.reshape(-1, step, step)

    結果を確認する．下段はやや六角形格子のようだが微妙である．

    In [12]:
    plt.figure(figsize=(4,4))
plt.suptitle("Dimensionally reduced activities (PCA)", fontsize=14)
for i in range(9):
    plt.subplot(3,3,i+1)
    plt.imshow(H_pca[i], cmap='turbo')
    plt.axis('off')
plt.tight_layout()

    自己相関マップ(autocorrelation map)を確認する．

    In [13]:
    plt.figure(figsize=(4,4))
plt.suptitle("Autocorrelation map (PCA)", fontsize=14)
for i in range(9):
    plt.subplot(3,3,i+1)
    corr = correlate2d(H_pca[i], H_pca[i], mode='same', boundary='symm')
    plt.imshow(corr, cmap='turbo')
    plt.axis('off')
plt.tight_layout()

    非負値行列因子分解 (NMF: nonnegative matrix factorization)

    sklearn.decompositionのNMFを用いる．コード自体はPCAの場合とほぼ同じである．

    In [14]:
    nmf = NMF(n_components=9, init='nndsvdar', solver='cd', 
          random_state=0, max_iter=5000).fit(X.reshape(Np, -1))
    In [15]:
    H_nmf = nmf.components_.reshape(-1, step, step)

    結果を確認する．六角形格子が生まれていることが確認できる．

    In [16]:
    plt.figure(figsize=(4,4))
plt.suptitle("Dimensionally reduced activities (NMF)", fontsize=14)
for i in range(9):
    plt.subplot(3,3,i+1)
    plt.imshow(H_nmf[i], cmap='turbo')
    plt.axis('off')
plt.tight_layout()

    自己相関マップ(autocorrelation map)を確認する．

    In [17]:
    plt.figure(figsize=(4,4))
plt.suptitle("Autocorrelation map (NMF)", fontsize=14)
for i in range(9):
    plt.subplot(3,3,i+1)
    corr = correlate2d(H_nmf[i], H_nmf[i], mode='same', boundary='symm')
    plt.imshow(corr, cmap='turbo')
    plt.axis('off')
plt.tight_layout()

    上記のようにPCAとNMFでは結果が異なる．大きな違いは圧縮表現が正負(PCA)か非負(NMF)か，という点である．神経活動は非負であるという特徴があり，こうした点からNMFの方が妥当であると言える．じゃあなぜ六角形になるのか，というのは現時点で自分が簡単に説明できないのでもう少しsurveyしてから本には書こうと思う．

    Tags Python  Neuroscience