時空間神経ダイナミクスの数理モデリングと脳型計算機応用に関する研究

算機応用に関する研究

著者

守谷 哲

学位授与機関

Tohoku University

学位授与番号

11301甲第19244号

時空間神経ダイナミクスの

数理モデリングと脳型計算機応用

に関する研究

Mathematical Modeling of Spatio-temporal Neuronal

Dynamics and its Application to Brain-morphic Hardware

東北大学 大学院 工学研究科 電子工学専攻

博士課程後期

₃

年の課程

守谷 哲

⟨ABSTRACT⟩

Mathematical Modeling of Spatio-temporal Neuronal

Dynamics and its Application to Brain-morphic Hardware

Satoshi Moriya

The human brain performs advanced information processing, such as cognition and decision making, with only 20 W of metabolic energy. Elucidating the neural basis of brain functions is thus not only important in biology and medicine but also in engineering. In fact, the brain serves as a role model for solving the power consumption issues in current computing technologies based on the von Neumann architecture. In designing such “brain-morphic”

hardware, it is critically important to mechanistically understand how the neural connectivity underpins various types of information processing that is realized within the brain. Artificial neural networks are the mathematical models of neuronal networks in the brain. In the field of machine learning based on artificial neural networks, deep learning has led to a major breakthrough. As hardware to execute deep learning, e.g., Tensor Processing Unit developed by Google, and Jetson developed by NVIDIA, have begun to spread. However, current neuromorphic system has not yet fully utilized the spatiotemporal dynamics of neural spikes for information processing. In addition, the network connectivity is different from the real brain.

Recent studies on the neural connectivity in the brain revealed the prevalence of the so-called “modular” structure, at both microscopic and macroscopic scales. A modular structure is a network structure in which a subset of nodes forms strongly interconnected subgroups (referred to as modules) which then are interconnected sparsely with others. Such organization is universal in real-world networks and is known to be observed in social, engi-neering, and biological networks. The present work is aimed to advance the understanding of the relationship between the structure and dynamics in biological neuronal networks in order to elucidate how connectivity defines the spatiotemporal patterns of neural activity for future engineering applications. In particular, I analyzed and discussed the influence of modular structure on the neural activity based on computational simulations and graph theoretical analysis. In parallel with this study, the dynamics of living neuronal networks have been experimentally studied in dissociated culture, in which the network structures are manipulated by growing the neuronal cells on microfabricated substrates. By fusing such ex-periments and computational simulations, I aim at constructive understanding of the effect of network structure on neuronal dynamics. This Thesis consists of six chapters. Chap-ter 1 explains the background and purpose of this research, as well as a concise review on preceding research on neuromorphic computing and properties of biological neural networks. Chapter 2 describes the generative models for random modular networks, and presents the results of a detailed graph theoretical analysis of network structures. In particular, assuming

a simple case in which the connection densities are uniform for all intra-modular and inter-modular couplings, a mean-field approximation of the connectivity matrix is derived, which enables the analytical formulation of the relationship between the structural parameters and the structural features of networks. The analytical values of modularity, clustering coefficients, and small-worldness agreed well with numerical ones of networks generated from the same configuration parameters. In particular, the analysis show that the modularity can be controlled linearly by a single parameter, i.e. the ratio of the number of intra-modular connections to total connections, while maintaining constant the other parameters such as to the connection degree or the number of neurons. This finding is important in discussing the relationship between modularity, dynamics, and information processing capability in neuronal networks.

In Chapter 3, the effect of network modularity on dynamics is computationally investigated using spiking neural network models. A particular focus has been devoted to the analysis of spontaneous activity, which is generated intrinsically in neuronal networks in the absence of any external stimuli, synchronous bursting activity, which characterizes the spontaneous activity of cultured neuronal networks, and its spatial complexity. The dynamics of network nodes are modeled using the leaky-integrate-and-fire (LIF) neurons which are extended to model the behavior of cultured neurons. The synchronous firing frequency of modular net-work increased approximately seven times compared to a random netnet-work. This is to due strong activities triggered by densely-connected neurons inside modules, which then propa-gate throughout the whole network. Furthermore, introducing inhibitory neurons increased the spatial complexity of the dynamics. This effect was particularly profound in networks with high modularity. These results suggest that modular structures play an important role in generating complex dynamics in biological networks with limited resources.

Chapter 4 investigates the relationship between spatiotemporal dynamics in a cultured neuronal network model and the performance of time-series information processing. For this, the framework of reservoir computing is applied to the dynamics of cultured neuronal networks to perform a speech recognition task. By using random and modular networks as the reservoir layer, speech recognition accuracy of approximately 80 to 90% was achieved in both networks. Contrarily to our expectation, however, difference in the learning and classi-fication performances was not observed for reservoirs with different network structures and functional complexity. Introduction of a learning rule that takes into account the asymmetry of input/output connection and the temporal spike information to utilize the dynamics of modular network is likely to be necessary to exploit the impact of internal organization and dynamics within a reservoir.

Chapter 5 reviews the results of several studies devoted towards a hardware realization of neural network functions, such as ultra-low power Izhikevich neuron circuits and synaptic elements. It was shown that brain-morphic calculations could be performed using a non-volatile analog memory as a synapse, and multiple neuronal spike-like voltage waveforms were obtained for the neuron circuit with the same physical implementation. These results

contribute to reducing the area and power consumption of brain-morphic hardware.

Chapter 6 is the conclusion. The results presented in this Thesis are summarized, along with their current limitations, potential solutions, and future perspectives. In summary, this Thesis reveals the effect of network connectivity on neural spatiotemporal dynamics by integrating cell culture experiments and computational simulations, and discusses the application of neural dynamics to time-series information processing to realize the efficient brain-morphic hardware. Elemental technologies required for hardware implementation of such functions is also described. These achievements is expected to advance the practical design of next-generation brain-morphic hardware.

目次

第1章 序論 1 1.1 本研究の背景 . . . 1 1.2 脳型計算機の開発動向と課題 . . . 2 1.3 脳神経回路における情報表現 . . . 3 1.4 脳神経回路の構造的特徴 . . . 4 1.5 神経情報処理の解明に向けた生理学実験 . . . 5 1.6 本研究の目的 . . . 6 1.7 本論文の構成 . . . 6 第2章 モジュール構造型神経ネットワークの構築とネットワーク構造の平均場 近似解析 7 2.1 本章の概要と目的 . . . 7 2.2 ネットワークモデル . . . 8 2.3 ネットワークの平均場近似 . . . 11 2.4 モジュール性 . . . 13 2.5 クラスタ係数 . . . 15 2.6 スモールワールド性 . . . 18 2.7 結言 . . . 24 第3章 モジュール構造型神経ネットワークにおける時空間神経ダイナミクスの 解析 25 3.1 本章の概要 . . . 25 3.2 本章の目的 . . . 25 3.3 ニューロンモデル . . . 26 3.4 シミュレーションプログラム . . . 29

3.5 ネットワークダイナミクスの解析手法 . . . 34 3.6 モジュール性とネットワークダイナミクス . . . 36 3.7 抑制性結合を有するモジュラーネットワークのダイナミクス . . . 43 3.8 結言 . . . 48 第4章 時空間神経ダイナミクスに基づく時系列信号分類 49 4.1 本章の概要と目的 . . . 49 4.2 Reservoir computing . . . 49 4.3 音声識別課題 . . . 50 4.4 ネットワークダイナミクスと音声識別性能の関係 . . . 55 4.5 結言 . . . 65 第5章 脳型情報処理のハードウェア実装に向けた要素技術の開発 66 5.1 本章の概要と目的 . . . 66 5.2 スピントロニクスメモリデバイスを用いた連想記憶システム . . . 66 5.3 アナログスパイキングニューロン回路 . . . 71 5.4 結言 . . . 78 第6章 結論 79 謝辞 81 受賞・表彰等 82 本研究に関する発表 83 参考文献 88

第

₁

章

序論

1.1

本研究の背景

近い将来，あらゆる人とモノがネットワークにつながるIoT (Internet of Things)の時 代が到来することが予期されている．これによって，利用可能な情報が爆発的に増加し， 分野を横断した知識や情報の連携が図られることによって，これまで想定されていなかっ たデータ利用が実現され，より豊かで暮らしやすい社会が実現されることが期待される． 一方で，膨大な情報から価値のある情報を取捨選択して分析する作業は未だ人の手によっ て行われていることが多く，作業量に限界があった．このような背景から，脳に学んだ情 報処理である脳型情報処理の重要性が益々高まっている． 脳型情報処理の重要な特徴の一つは，パターン認識や判断などの高度な情報処理を従来 のノイマン型コンピューターでは達成不可能な低消費電力で実行できる点である．これ は，情報を並列かつ分散的に処理していることに起因する．分散的な情報処理の利点とし て，高いクロック周波数を必要としない為に電力消費が抑えられることや，メモリが分散 してプロセッサの付近に配置されることによって，メモリとプロセッサの通信速度に計算 が律速されるノイマン・ボトルネックを緩和することが出来ることなどが挙げられる．ま た最近では，人工ニューラルネットに基づく機械学習の分野において，Deep learning（深 層学習）が大きなブレイクスルーをもたらした．Deep Learningは多層に接続したニュー ラルネットワークにより情報処理を行う機械学習手法である．多層に接続されたニューラ ルネットワークの先駆的な研究として福島らが提唱したNeocognitron [1]がある．しか し，計算量が膨大であったため，社会への普及はしなかった．その後，Hintonらの研究 により多層ニューラルネットワークを学習する際の勾配消失や局所最適化へのトラップな どの問題の解決が示された [2, 3]．計算機の演算性能やメモリの増大，web の普及に伴う

学習データの調達の易化なども手伝って，Deep LearningはIoT 社会において必須な画 像認識や情報の検索・分析などのアプリケーションに革新をもたらした．

1.2

脳型計算機の開発動向と課題

画像認識や異常検知などの，ある特定の情報処理を実行する際には，汎用性の高い CPUで計算するよりも，性能を制限してその情報処理に特化したハードウェアを利用し た方が消費エネルギーや高速化の観点から有利である．このような背景から，GPGPU

(General purpose computing on graphics processing units)の利用や脳型計算に特化し

た専用ハードウェア（：脳型ハードウェア）の開発が精力的に行われている．前述のDeep

Learningは，荷重値とノード出力の積和演算の計算負荷が大きいため，多数の積和演算

器を搭載したハードウェアで実行すると計算効率が良い．このようなハードウェアとして

Google社が開発したニューラルネットワーク演算専用のLSIである Tensor Processing Unitや，NVidia社が開発したDeep Learningの実行に特化したJetsonなどが市場に表 れ始めている． Deep Learningは画像認識や分類などで驚異的な性能を誇るが，その反面，時系列情報 などを処理する場合は必ずしも効率的とは言えない．これは，Deep Learningが主に画像 などの静的な情報を入力とし，多層ニューラルネットにより処理された情報を静的な空間 パターンとして表現する為である．特に構造に工夫を加えない限りは，ネットワークの出 力は現在の入力によってのみ決定され，過去の情報は反映されない．一方で，生体の脳は 過去に入力された時系列情報を反映した情報処理を行っている．これは，脳が図1.1に示 すように，再帰的な結合構造を持ち，神経細胞の時間的・空間的な発火パターン，すなわ ち時空間ダイナミクスにより情報表現・情報処理を行っているからに他ならない．時空間 ダイナミクスを積極的に活用することで，すなわち時間軸にも情報を展開することで，よ り少ないリソースで高い情報表現能力や情報処理能力を獲得しうる可能性がある． 神経細胞の発火ダイナミクスや再帰的な構造を取り入れたハードウェアとしては，Intel

社が開発したLoihi [4]やIBM社が開発したTrueNorth [5]などが挙げられる．これらの 脳型ハードウェアでは，神経細胞のスパイクダイナミクスを再現し，情報処理を行うこと が可能である．しかし，発火の時空間ダイナミクスを効率的に情報処理に活用する為のア ルゴリズムは未だ開発されておらず，社会への普及には至っていない．次節では，神経発 火による情報処理について，これまで生理学実験から得られている知見を整理する．

図1.1 階層型ニューラルネットワークにおける空間パターンによる情報表現・情報処 理（上図），再帰型ニューラルネットワークにおける時空間パターンによる情報表現・ 情報処理（下図）

1.3

脳神経回路における情報表現

脳は約1000 億の神経細胞とそれらを繋ぐ100兆のシナプスによって構成された複雑 ネットワークである．この多細胞ネットワークを物理的基盤として，認知・意思決定と いった高度な情報処理がわずか20 Wという消費電力で実現されている [6]．神経回路網 において，情報の大部分はニューロンが脱分極して活動電位を生じさせる発火と呼ばれる 現象によって表現されていると考えられている．しかし発火がどのように情報を表現して いるかは未だ明らかになっておらず，複数の情報表現方法が提唱されている [7]． 脳の情報表現として最も代表的なものは発火頻度表現と呼ばれるもので，細胞の発火率 が情報をコーディングしているという考え方である．発火頻度表現の主な情報表現方法と して，細胞1つ1 つが外界の情報に特化して対応する情報に対し強く発火するといった おばあさん細胞仮説や，多数の細胞が外界の情報に対してゆるやかな選択特性を持ってい て，その合算として情報表現が行われているというポピュレーションコーディングなどが ある． もう一つの情報表現の可能性として，各細胞の発火におけるミリ秒オーダーのタイミン グの差が情報表現において重要であるとするタイミング表現がある．タイミング表現には 各細胞の発火タイミングの時間差に情報が表現されているという時間差表現や，単一ある いは複数の細胞の発火タイミングの間隔が情報を表現しているという発火パターン表現な どがある．他にも，多数の細胞が数ミリ秒のオーダーで同期して発火することが情報を表

現しているとする同期表現がある．この発火の同期現象は大脳皮質で見られ，記憶や注意 などに関わっていると言われており [8, 9]，視覚野などにおいてもカラム間の同期的な発 火が視覚情報の統合に関連していることが示唆されている [10]．また，さまざまな空間的 に分離した領野間の情報の統合を行っているとの仮説もあり[11]，同期発火現象は脳機能 の実現において重要な役割を担っていると言える．

1.4

脳神経回路の構造的特徴

発火の時空間ダイナミクスを規定する大きな要因として，ネットワークの接続構造が挙 げられる．顕微鏡やfMRI 等を用いた脳構造の研究によって，脳は様々なスケールにおい てモジュール構造（または階層的モジュール構造）と呼ばれる構造を有していることが明 らかになってきた[12]．モジュール構造とは，ネットワークの一部分がネットワーク全体 に比べ密に結合することでモジュールを形成し，それらのモジュールが疎に結合している 構造である [13]．階層的モジュールは，モジュールの中をさらにサブモジュールとして分 割できる入れ子のような構造である．モジュール構造は，社会や工学，生物などの全く異 なる分野における様々な複雑ネットワークに現れる普遍的な構造である． ヒトの脳に見られるモジュール構造として，全脳における階層的モジュール構造を示 す．図1.2は拡散MRIから得られたヒトの構造的結合行列から，モジュール構造を推定 したものである [14]．図中の樹形図は各ノードが含まれるモジュールを示しており，樹形 図の下に示された数値は脳をいくつのモジュールに分割したかを示している．数値は樹形 図を分割する赤線に対応しており，赤線が下に行くほど細かい分割となる．例えば図の左 上の最も荒い分割では，前頭葉や頭頂葉などの分類と良く対応している．このように，脳 はモジュールの中にいくつものサブモジュールが含まれるような階層的モジュール構造を 有している． 脳をより細かいスケール(10−3 m) で見たときにも，大脳皮質にはカラム構造と呼ば れるモジュール構造とよく似た特徴的な構造が存在する．1つのカラムは直径300− 500 textrmµm程度の面に対して垂直な円柱状の細胞集団で，カラム内の細胞集団は共通の入 力を受ける．カラムの構造は大脳皮質全体にわたって基本的に同じであるが，その活動は 部位ごとに大きく異なる．例えば視覚野のV1野の各カラムは視覚情報の線分の傾きに対 してそれぞれ特異的な反応を示す（反応選択性）[15]．このような反応選択性は，入力に 対する神経網の後天的な学習により獲得されるものである．興味深いことに，大脳皮質の サイズが大きく異なる他の動物においても，カラムの物理的サイズは大きく変わらない． よって，このカラムは脳の情報処理の最小単位であると考えられている．これらのカラム

図1.2 脳の階層的モジュール構造，および拡散MRIよりボクセル間のモジュール構 造を推定して得られた樹形図．文献 [14]より引用． （≒モジュール）が成す結合構造が，高度な情報処理を生み出していると考えられる．

1.5

神経情報処理の解明に向けた生理学実験

脳研究のアプローチの一つとして，生きた状態の神経細胞が構成する回路構造やその活 動を研究する領域がある．脳の活動を測定する技術として，代表的なものにfMRIがあ る．fMRIは核磁気共鳴現象を利用した測定方法であり，水素原子の磁場に対する共鳴現 象を利用する．脳の活性化した領野において酸素が多く消費されるが，この酸素が多く消 費された領野において，一時的に血流量が増加する．この関係を利用して，脳の血流量変 化をMRI で捉えることにより，脳の活動量を定量する方法である．大きな特徴として， 測定の方法上開頭する必要が無く非侵襲的であること，また脳の広い領野にわたって測定 することが出来るため，領野間の活動の相関などを調べることが出来ることなどが挙げら れる．一般的なfMRIの空間分解能は一辺が数ミリ程度のボクセルで，時間分解能は１フ レーム当たり数十_∼数百 msである [16]．侵襲的な方法としては，脳に電極を取り付け， 電気的な活動を記録する方法や，脳を直接Ca蛍光観察するなどの方法がある． 一方で，脳組織を体外に取り出し，単一の神経細胞になるまで分離してから培養皿上で 培養する分散培養という手法がある．分散培養の利点として，神経回路の形成過程の観察 やその活動の変化を観測できる点や，培養の温度や培養液の組成などの，細胞の生育環境 を自由に制御し規定できる点などが挙げられる．一方で，培養神経回路を用いた実験の課 題として，例えば実際の脳神経回路との活動の相違がある．1.4節で述べたように，脳で は三次元に配置された神経細胞がモジュラーなネットワークを形成しているのに対し，培

養神経回路では構造上の制約より二次元的な回路形成が行われ，また細胞間の結合は基本 的にはランダムである．このため，単純に細胞を分散させて形成した回路に見られる活動 は実際の脳と大きく異なる [17]．この課題を克服するために，培地の構造や物理的性質を 制御して細胞が形成する回路構造を制御する試みが行われている [18].これらの技術によ り，生体に近い構造をin vitroで作りだすことが可能になりつつある．

1.6

本研究の目的

本研究の目的は，次世代脳型ハードウェアの創出に向け，神経回路網の構造が時空間ダ イナミクスに与える影響を明らかにし，神経時空間ダイナミクスの時系列情報処理応用に ついて検討することである．生理学的に妥当なニューロンモデルを結合して構成した神経 回路のシミュレーションと，数理モデルから得られた結果を検証するための生理学的な実 験系を組み合わせ，神経回路の構造が時空間ダイナミクスをどのように制御しているかを 構成論的に明らかにする．特に生体神経回路を含む様々な複雑ネットワークにおいて普遍 的に存在するモジュール構造が，神経回路網の活動に与える影響を計算論的神経科学の 方法論に基づいて解析する．また神経時空間ダイナミクスを音声認識課題に適用し，モ ジュール構造と情報処理性能の連関について明らかにする．更に神経時空間ダイナミクス を効率的に実行するためのハードウェア実装に向けた要素技術の開発について示す．

1.7

本論文の構成

本論文の構成は次の通りである．第２章では培養神経回路のネットワーク構造を再現す るネットワークモデルについて概説し，ネットワーク構造と構造パラメータ間の数学的関 係を明らかにする．第３章ではネットワーク構造が神経回路ダイナミクスに及ぼす影響を 培養神経回路実験とシミュレーションの両面から明らかにする．第４章では神経回路ダイ ナミクスを用いた音声識別課題から，ダイナミクスと情報処理性能の関係について考察す る，第５章では時空間神経ダイナミクスのハードウェア実装に向けたニューロン回路やシ ナプスメモリ等の構成要素について検討する．第６章にて本論文の結論を述べる．

第

₂

章

モジュール構造型神経ネットワーク

の構築とネットワーク構造の平均場

近似解析

2.1

本章の概要と目的

本章では，ランダムネットワーク及びモジュール構造を有したネットワークモデルにつ いて示し，グラフ理論に基づいてネットワーク構造の詳細な解析を行った結果について示 す．特に，モジュール内およびモジュール間の接続密度がすべてのノードペアで同等であ る単純なケースを仮定し，ネットワークの構造パラメータとネットワークの構造的特徴と の関係を，結合行列の平均場近似から解析的に導出できることを示す [19]． 本章の構成は以下の通りである．2.2節では，ランダムな結合構造及びモジュール構造 を有するネットワークの構成方法について説明する．2.3節では，ネットワーク構造の解 析の為に，ネットワークの結合行列を平均場近似した方法について説明する．2.4, 2.5, 2.6 節では，グラフ理論や複雑性ネットワークの分野で定式化されているmodularity,クラス タ係数，スモールワールド性などのネットワークの構造的特徴の理論式を導出し，実際に 構成したネットワークから数値解析的に求めたこれらの指標と比較する．最後に，2.7節 で本章を総括する．

図2.1 ランダムネットワークの構造．左：N = 20, ¯k = 8，中央：N = 100, ¯k = 17.5， 右：N = 200, ¯k = 25

2.2

ネットワークモデル

本節では，ネットワークの詳細な構成方法を述べる．ネットワークの構成を行うにあ たって，基本的な要素であるニューロン（またはノード）の数をN，各ニューロンが平 均で形成するシナプス結合（またはエッジ）の数である平均結合次数を¯_{k とする．この} 時，結合は方向性を持つとした．また特に明示しない限り，多重結合および自己結合を許 容した．

2.2.1

ランダムネットワークの構成

1から N までの範囲の整数をランダムに二つ決定する．その整数をそれぞれi, j とし た時，ニューロンiからニューロンj へ向かう結合を形成する．このプロセスを，ネット ワークの総結合本数である T = ¯kN, (2.1) になるまで繰り返す．図2.1に構成したランダムネットワークの構造の例を示す．

2.2.2

モジュラーネットワークの構成

モジュール構造は，密な結合を有するモジュール同士が疎に結合している構造である． 以下に，本シミュレーションで使用したモジュラーネットワークの構成方法を示す．初 めに，ネットワークを構成するモジュールの数M を決める．その後，各モジュールに同 数のニューロンを割り当てる．N がM で割り切れる場合は，各モジュールにN/M の ニューロンが存在することになる．次に，ネットワークのモジュール化の度合いを制御す

るためのパラメータであるモジュール内結合確率pin を決定する．結合が同一モジュール 内で形成されるか，または異なるモジュール間で形成されるかを，50%の確率でランダム に選択し，同一モジュール内に結合が形成される場合には，同じモジュールに属している ニューロンのペアが選ばれるまでランダムにニューロンを選ぶ操作を繰り返す．最終的に 選ばれた同一モジュールに属するニューロン間に，pin の確率で結合を形成する．異なる モジュール間に結合が形成される場合は，同様に異なるモジュールに属するニューロンの ペアが選ばれるまでランダムにニューロンを選ぶ操作を繰り返し，選ばれたニューロン間 に1− pin の確率で結合を形成する．この操作を，T の結合が形成されるまで繰り返す． なお，pin = 1 の時には完全に分離したM 個のネットワークが形成され，pin = 1/M の 時にはランダムネットワークと等しくなる． このアルゴリズムの擬似コードを Algorithm 1 に示す．このアルゴリズムの特徴は， ネットワーク全体の最終的な結合数T を保ちつつ，各ニューロンの平均的な結合次数_k¯ が保存されることと，結合が形成されるかは確率的である為，最終的に構成されるネット ワークの同一モジュール内結合と異種モジュール間結合の比にばらつきが生じることの2 点である．図2.2にこのアルゴリズムで構成したネットワークの例を示す．

Algorithm 1 モジュラーネットワークの構成アルゴリズムの擬似コード

1: procedure ModularNetworkConstruction(N , M , ¯k, pin)

2: Initialize elements of connectivity matrix wij to zero

3: Allocate a module label m to each node

4: counter ← 0

5: repeat

6: x← Random int value INTRA or INTER

7: if x == INTRA then

8: repeat

9: i, j← Random int value from 1 to N

10: until mi= mj

11: if Random double value from 0 to 1 =< pin then

12: wij ← wij + 1 13: counter ← counter + 1 14: end if 15: end if 16: if x == INTER then 17: repeat

18: i, j← Random int value from 1 to N

19: until mi̸= mj

20: if Random double value from 0 to 1 > pin then

21: wij ← wij + 1 22: counter ← counter + 1 23: end if 24: end if 25: until counter == ¯k× N 26: end procedure

図2.2 モジュラーネットワークの構造．N = 100, ¯k = 17.5, M = 4左：pin = 0.25， 右：pin= 0.98．

2.3

ネットワークの平均場近似

ネットワークの構造を評価するために，グラフ理論や複雑ネットワークの分野で定式化 されている方法を用いてネットワーク構造の解析を行った．これらの解析を行うことに よって，ネットワークの構造的特徴を反映したスカラ量やベクトル量を得ることが出来 る．本研究では，指標としてモジュール性，クラスタ係数，平均経路長，及びクラスタ係 数と平均経路長を組み合わせて計算されるスモールワールド性について評価した．これら の指標はネットワークの大域的な構造を評価する代表的な指標である．これらのネット ワーク構造の解析指標については後の節で詳述するが，その前に，ネットワークの平均場 近似の方法についてここで説明する．

モジュラーネットワークを構成するアルゴリズムとして，stochastic block model [20, 21]が広く知られている．このモデルでは，ネットワークのノード数N，モジュール 数M，各モジュール内・モジュール間における結合確率行列Ψ = [ψmi,mj]を規定するこ とによりモジュラーネットワークを構成する．結合確率行列はM × M の行列で，行と 列はモジュールのラベルに対応する．要素ψmi,mj はモジュールmi，モジュールmj に 属するニューロンが結合する確率に対応する．この行列に従い結合を形成することで，各 モジュール内及びモジュール間の結合密度を自由に調整することが出来る．すべてのモ ジュールのサイズが等しい場合，本研究で用いたモジュラーネットワークの構成アルゴリ ズムはすべてのモジュール内またはモジュール間でそれぞれの結合確率が等しくなる． 前節で説明したアルゴリズムで構成されるネットワークの結合行列をB ={bij}とし，

図2.3 ネットワークの結合行列と平均場近似した結合行列の概略図 結合行列の平均場近似をA = {aij}とする．ネットワークの結合関係に焦点を当てて解 析を進めるため，aij はbij の要素が1以上となる確率に対応させる．すなわち， aij = { 0 (bij = 0) 1 (bij ̸= 0) (2.2) で示される二値化の操作を行う． モジュール内，モジュール間の結合密度がすべてのモジュールで同一の場合，平均場近 似した結合行列Aの要素aij はstochastic block modelにおける結合確率行列Ψの要素

ψmi,mj と以下のように対応する． aij = ψmi,mj = { 1− (1 − pinM/N2) ¯ kN_{) ≡ P} I (mi = mj) 1− (1 − M(1 − pin)/N2(M − 1)) ¯ kN _{≡ P} D (mi ̸= mj) (2.3) ここで，PI 及びPD はそれぞれ同一モジュール(mi = mj)と異なるモジュール(mi ̸= mj)に属するニューロン間のaij の値である．図2.3に，平均場近似した結合行列の概略 図を示す．この平均場近似した結合行列を元に，構成したネットワークの各構造評価指標 の期待値を導出する．

2.4

モジュール性

2.4.1

Modularity

の定義

Modularityはネットワーク全体の構造がどの程度のモジュール性を有しているかを定 量する指標である[13]．バイナリ有向グラフのmodularity Qは以下の式で計算される． Q = 1 T ∑ i,j [ ai,j − k_iink_jout T ] δmi,mj (2.4) ここで，k_iin =∑_jaji はニューロンiの入力次数, kiout = ∑ jaij はニューロンiの出力 次数,δはクロネッカーのデルタ，mi はニューロンiが属するモジュールのラベルである. ネットワークの持つモジュール性が高いほどQは1に近づき，ランダムネットワークに 近い場合，Qは0に近づく． このQを最大化するようにmi を決定したとき，Qの最大値がネットワークの mod-ularity に対応する．mi の組み合わせは指数関数的に増加するため，Q を最大化する ようなmi の組み合わせを求める最適化問題を解く必要がある．この解法として広く用 いられているのがLouvain法 [22] であり，本研究で用いた MATLABのツールボック スであるBrain Connectivity Toolbox [23](URL:https://sites.google.com/site/bctnet)

のmodularityを求めるプログラムもLouvain法により計算している．

2.4.2

Modularity

の理論値

クロネッカーのデルタδmi,mj の値は mi = mj の時に1，mi ̸= mj の時に0となる． 即ち，式2.4の右辺の総和を計算する際には同一モジュール内の結合のみを考慮すれば良 い．式2.4におけるネットワーク全体の結合数T の期待値は E[T ] = N (NIPI + NDPD), (2.5) で計算される．ここで，NI = N/M は同一モジュールに属するニューロンの数，ND = N − NI = N (M − 1)/M は異なるモジュールに属するニューロンの数である．E[·]は期 待値であることを示している．上記の式は，個々のモジュールに含まれるニューロンの数 が同一であるという仮定の下で成り立つ． Algorithm 1 に示されるアルゴリズムでネットワークを構成した場合，式2.4右辺の

表2.1 ネットワークパラメータ Parameter M 2, 4, 10 N 100, 500, 1000 ¯ k 10, 50, 100 pin [1/M, 1] kin_とkout_{はどちらも平均がk¯の二項分布に従う．その積の期待値は} E[kinkout] = (NIPI + NDPD)2+ PI(1− PI), (2.6) と求められる．上式の右辺第二項は，相関を持つ項同士の積で生じる誤差を補正する為の 項である．式2.5, 2.6を2.4に代入すると，Qの期待値は， E[Q] = 1 E[T ] ∑ mi=mj [ ai,j − E[k_iink_jout] E[T ] ] δmi,mj = 1 M (NIPI + NDPD) ( N PI − E[PI(1− PI)] NIPI + NDPD ) −1/M, (2.7) となる.なおこの時，同一モジュールに含まれるニューロンのペアの数はM N_I2 = N2/M である．

2.4.3

Modularity

の数値解との比較

ここでは，式2.7から求めたmodularityの解析解と，Algorithm 1で示したアルゴリ ズムから構成したモジュラーネットワークの持つmodularityを式2.4に基づいて数値計 算した値を比較した．数値解は，100 回のシミュレーションの平均値を用いた．シミュ レーションで用いたネットワーク生成パラメータを表2.1に示す．なお，代表的な値は太 字で示した． 図2.4(a)に，結合次数が異なるネットワークのQとpinの関係を示す．この時，他の パラメータは太字で示したデフォルトパラメータ（N = 100，M = 4）とした．数値計算 の結果は，pinが高い領域での理論曲線とよく一致した．0.95 ≥ pin ≥ 0.5の範囲の平均 Err[Q]は，k = 10, 50, 100¯ でそれぞれ3.4%，3.3%，4.4%となった． 一方で，低いpin での誤差は大きくなった．通常，ネットワークのQの評価では，各

ルゴリズム [22]によりQを最大化する最適なモジュール割り当てが新たになされる．こ の各モジュールへのノードの割り当ての乖離が，低pin での誤差の原因である．一般にモ ジュール検出はNP困難な問題であり，モジュール構造があいまいな場合，こちらの意図 しないモジュールの割り当てが生じる [24]．事前に想定したモジュールの割り当てから Qを計算する場合，すなわち，モジュール割り当てを固定し，生成したネットワークの 結合行列を用いて式2.4からQを計算した場合，数値結果は理論曲線と完全に一致した． （図2.4，塗りつぶしマーカー）．従って，式2.7の理論式は事前に想定したモジュール割 り当ての下でのQを表している．数値計算の結果は，N と¯k を一定に保ちながらM を 変化させた場合の理論曲線ともよく一致した（図2.4(b)）．pin ≈ 1/M でのErr[Q]の増 加は，モジュールM の数に関係なくすべてのネットワークで観察された．一方，高いpin でのErr[Q]はM に比例して増加した．さらに，ネットワークのサイズが変わってもモ ジュール内およびモジュール間の接続密度は変わらないため，QはN に関係なく同様の 値を持つことがわかった（図2.4(c)）．

2.5

クラスタ係数

2.5.1

クラスタ係数の定義

クラスタ係数はネットワーク中の三つのニューロンの結合である三体結合が，どの程度 存在するかを定量する指標である．Modularityがネットワークの大域的な集団を定量す る指標であるのに対し，クラスタ係数は局所的な密結合の程度を定量する．クラスタ係数 はそのネットワーク実現されている三体結合の数を，実現可能な最大の三体結合数で割る ことで得られる．この指標が1に近いほど，ネットワーク中の三体結合の割合が高いこと を示す．

2.5.2

クラスタ係数の理論値

ニューロンiにおけるクラスタ係数Ciは前述の定義から，以下の式で計算される [25]． Ci = (A + AT)3_ii 2[ktot i (ktoti − 1) − 2k↔i ] . (2.8) ここで，AT _{はAの転置行列，(A + A}T₎3 ii は(A + AT)3 のi 番目の主対角成分である． ネットワーク全体のクラスタ係数C は各ニューロンのクラスタ係数の平均から求められ

図2.4 Modularity Qとpinの関係．実線は解析値，白抜きのマーカーは数値解の平 均値．マーカーは式2.4において，事前に割り当てたモジュールのラベルにおいて計 算した数値解の平均値．(a) 結合次数が異なるネットワーク (¯k = [10, 50, 100], N = 100, M = 4). (b) モジュール数が異なるネットワーク(M = [2, 4, 10], N = 100, ¯k = 10). (c) ニューロン数が異なるネットワーク(N = [100, 500, 1000], M = 4, ¯k = 10)． 文献 [19]より引用．

る．k_itot は入力次数と出力次数の和であり，以下の式で計算される． k_itot =∑ j̸=i aji+ ∑ j̸=i aij (2.9) k_i↔はニューロンiにおける双方向結合の本数であり，次の式で計算される． k↔_i =∑ j̸=i aijaji (2.10) 仮定した平均場近似の結合行列は対称行列なので，(A + AT)ij は2aij と等しくなる． よって，式2.8右辺の分子は (A + AT)3_ii = 8∑ j ∑ h aijaihahj, (2.11) となる．上記の式は，3つのニューロンの属しているモジュールを考慮することでより簡 単化できる．3つのニューロンの属しているモジュールの関係は，3つのケースが考えられ る．1つはすべてのニューロンが同じモジュールに属している場合，即ちmi = mj = mk で，この場合は aijaihahj = P_I3 となる．この状態にあるニューロンの組み合わせは (NI − 1)(NI− 2)通り存在する．次に，2つのニューロンが同じモジュールに属している 場合，即ちmi = mj ̸= mh で，この場合はaijaihahj = PIPD2 となる．このような組み 合わせは3ND(NI − 1)通り存在する．最後はすべてのニューロンが異なるモジュールに 属している場合，即ちmi ̸= mj ̸= mh で，この場合はaijaihahj = PD3，組み合わせの数 はND(ND− NI)である． 次に，クラスタ係数の分母を考える．クラスタ係数の分母の期待値は以下のようになる．

E[2{ktot_i (ktot_i − 1) − 2k↔_i }] = 2{E[(ktot_i )2]− E[k_itot]− 2E[k_i↔]}, (2.12)

ここで， E[k_itot] =∑ j̸=i aji+ ∑ j̸=i aij = 2{(NI − 1)PI + NDPD}, (2.13) E[k↔_i ] =∑ j̸=i aijaji= aij = (NI − 1)PI2+ NDPD2, (2.14) E[(k_itot)2]は自己相関を持つ為，分散を考慮する必要がある．式で表すと,

E[(k_itot)2] = E[(ktot_i )]2+ V [(k_itot)], (2.15)

となる．モジュールのサイズがすべて同じ場合を仮定しているため，すべてのノードのク ラスタ係数は同一となる．結果として，最終的なクラスタ係数Cは E[C] = (NI − 1)(NI − 2)P 3 I + 3ND(NI − 1)PIP 2 D+ ND(ND− NI)P 3 D {(NI − 1)PI + NDPD}2 − (NI − 1)PI2− NDPD2 (2.17) となる．モジュールのサイズが異なる場合は，同様にニューロンの属しているモジュール の関係とその数を計算すればよい．例えば，モジュールmに属しているニューロンiと同 じモジュールに属しているニューロンの組み合わせの数は(NIm− 1)(NIm− 2)となる．

2.5.3

クラスタ係数の計算値との比較

Modularityの比較と同様に，式 2.17から求めたクラスタ係数の解析解と，図1で示 したアルゴリズムから構成したモジュラーネットワークの持つクラスタ係数の値を比 較した．実験で用いたネットワークパラメータは，表2.1 と同じとした．図2.5(a) に， ¯ k = [10, 50, 100], N = 100, M = 4のネットワークにおけるpinとCの関係を示す．数値 解と理論曲線の間の相対誤差は，すべてのpinで2%未満だった．これは，クラスタ係数の 理論式2.17が，結合行列の確率表現から近似なしに導出されたためである．図2.5(b),(c) に示した，異なるモジュール数(M = 2, 4, 10)およびニューロン数(N = 100, 500, 1000) のネットワークにおいても，数値計算と理論曲線の一致が確認された．特に，ランダム グラフに対応するpin = 1/M の条件で，ネットワークのC はランダムネットワークに おけるクラスタ係数の理論値であるk/N¯ （ただし，k¯ ≪ N) に収束することを確認し た[26, 27]．

2.6

スモールワールド性

2.6.1

スモールワールド性の定義

ネットワークのスモールワールド性は，高いクラスタ係数と短い平均経路長によって特 徴づけられる性質で，WattsとStrogatzにより提唱された [26]．実世界の様々なネット ワークはスモールワールド性を有すると言われており，モジュラーネットワークも一般的 にスモールワールド性を有する．

図2.5 クラスタ係数Cとpinの関係．実線は解析値，白抜きのマーカーは数値解の平

均値．(a) 結合次数が異なるネットワーク(¯k = [10, 50, 100], N = 100, M = 4). (b)

モジュール数が異なるネットワーク (M = [2, 4, 10], N = 100, ¯k = 10). (c) ニュー ロン数が異なるネットワーク (N = [100, 500, 1000], M = 4, ¯k = 10).文献 [19]より 引用．

2.6.2

スモールワールド性の理論値

ネットワークの持つスモールワールド性を評価する指標であるsmall-worldness S は以 下の指標で求められる． S = C/Crand L/Lrand (2.18) ここで，Crand とLrandはそれぞれ同数のN , ¯k のランダムネットワークにおけるクラス タ係数と平均経路長に対応する．2.5節で既にクラスタ係数の理論値は導出したため，こ こでは平均経路長の理論値を導出する．平均経路長はすべてのニューロンペアの最短経路 の平均であり，以下の式で求められる． L = 1 N (N − 1) ∑ i,j̸=i Lij (2.19) ここで，Lij はニューロンij間の最短経路長である．バイナリネットワークにおける最短 経路とは，あるニューロンから目的のニューロンに到達するまでに経由した経路の数の中 で最小のものであり，直接の経路がある場合は最短経路長は1となる． 平均経路長の理論値は，ノードとエッジの増加に伴い通過可能な経路数が爆発的に増加 するため，ランダムネットワークでも導出が複雑で難しい [28]．そのため，最長の最短経 路を3とし，3を超える場合は3に丸め込む操作を行うことで簡単化した．この仮定の下 での，Lの理論値は以下の式で計算される． E[L] = 1 (N − 1){(NI − 1)(αI1+ 2αI2 + 3αI3) + ND(αD1+ 2αD2+ 3αD3)} (2.20) ここで，αIx とαDx はそれぞれニューロンが同じモジュール及び異なるモジュールに属 している場合の，最短経路長がLij = xになる確率である．Lij = 1の時，αI1 とαD1 は それぞれPI とPD に等しくなる． Lij = 2の場合，それらのニューロンがモジュールにどのように属しているかを考慮す る必要がある．出発点のニューロンi と目的のニューロンj が同じモジュールに属して いると仮定する．この場合，iと同じモジュールに属するニューロンhを経由したのちに j に到達するルート，即ちmi = mj = mh の場合と，iとは異なるモジュールに属する ニューロンgを経由したのちにj に到達するルート，即ち mi = mj ̸= mg の2 通りの ルートが存在する． mi = mj のニューロンペアi, j で経路長がxのパスが存在する確率をβIx とする．こ の時， βI2 = 1− {1 − PI2} NI−2{1 − P2 D} ND_{, (2.21)}

となる．αI2は経路長2のパスが存在し，かつ経路長が1のパスが存在しない確率に等し いので，以下のように求められる． αI2 = (1− αI1)βI2. (2.22) 次に，出発点のニューロン iと目的のニューロンj が異なるモジュールに属している と仮定する．この場合は，iと同じモジュールに属するニューロンhを経由したのちにj に到達するルート(mi = mh ̸= mj)，iとは異なるモジュールに属するがj とは同じモ ジュールに属するニューロンgを経由したのちにj に到達するルート(mi ̸= mj = mj)， iともj とも異なるモジュールに属するニューロンf を経由したのちにj に到達するルー ト(mi ̸= mj ̸= mj)の3通りのルートが存在する．先ほどと同様にmi ̸= mj のニュー ロンペアi, jで経路長がxのパスが存在する確率をβDxとすると，βD2は βD2 = 1− {1 − PIPD}2(NI−1){1 − PD2} ND−NI_{, (2.23)} この時，αD2は αD2 = (1− αD1)βD2， (2.24) となる．最短経路長が3である確率αI3 及びαD3は，仮定から最短経路長が1でも2で もない確率なので，αI3 = 1− αI1− αI2 とαD3 = 1− αD1− αD2となる．最終的に，平 均経路長の理論値は E[L] = 1 N − 1{(NI − 1)[3 − 2PI − (1 − PI)(1− {1 − P 2 I} NI−2{1 − P2 D} ND_)] + ND[3− 2PD− (1 − PD)(1− {1 − PIPD}2(NI−1){1 − PD2} ND−NI_)]} (2.25) 最終的に，式 2.18に式 2.17, 2.25を代入することによってスモールワールド性の理論 値が求められる．

2.6.3

スモールワールド性の計算値との比較

これまでの解析と同様に，スモールワールド性の解析解と，Algorithm 1で示したア ルゴリズムから構成したモジュラーネットワークの持つスモールワールド性の値を比 較した．実験で用いたネットワークパラメータは，表2.1 と同じとした．図2.6(a) に， ¯ k = [10, 50, 100], N = 100, M = 4のネットワークにおけるpin とS の関係を示す. 密 に接続されたネットワークの数値解（¯_{k = 100）は，理論曲線（Err[S] < 0.8%）と精度} よく一致した．ただし，疎に接続されたネットワークの理論値（k = 10¯ ）は，高いpin

において数値解から逸脱することが分かった．（pin＝0.95におけるErr[S] = 13%）．こ れは，スモールワールド性の理論式を導出するときにすべての Lij > 3 をLij = 3 と 近似したが，ネットワークの結合が疎になると，大きな Lij を持つニューロンペア ij の割合が増加する．従って，相対誤差は比例して増加する．この誤差は，理論式で考慮 されるLij の次数を増やすことで減少させることができるが，現在の近似でも十分に S − pin 曲線の全体的な傾向を捉えている．これは，S の計算において L が Lrand に 対する比率で表現されているため，絶対値が重要な意味を持たないためである．M が 増加すると，Lij > 3 のニューロンペアの数も増加し，これに伴って誤差も増大した （pin = 0.95, M = 10, Err[S] = 33%）．対照的に，M を固定したより大きなネットワー ク（N = 500, 1000）のErr[S]の増加は，N = 100のネットワークの相対誤差と同程度 のものであった(pin= 0.95, N = 1000, Err[S] = 13%）．これは，pin の増加によるLの 変動が比較的小さいためである．

図2.6 スモールワールド性Sとpinの関係．実線は解析値，白抜きのマーカーは数値解

の平均値．(a) 結合次数が異なるネットワーク (¯k = [10, 50, 100], N = 100, M = 4).

(b) モジュール数が異なるネットワーク (M = [2, 4, 10], N = 100, ¯k = 10). (c)

ニューロン数が異なるネットワーク (N = [100, 500, 1000], M = 4, ¯k = 10).文献[19]

2.7

結言

本章ではモジュラーネットワークの構成方法について示し，またグラフ理論に基づいて ネットワーク構造の解析を行った．特に，モジュラーネットワークの構造特性は，ネット ワークのパラメーターセット（N, M, ¯k, pin）から生成された平均場近似結合行列に基づ いて解析的に予測できることを示した．モジュール性Qは主にpinとM によって決定さ れ，N 依存性はわずかであった．疎に接続したネットワークの場合，Qはpin の線形関 数として近似される．クラスタ係数C は，pin に非線形に依存し，N, ¯k, M の値にも強く 依存する．これまでに述べたように，モジュラーネットワークはスモールワールド性を有 している [26]．しかし，Qはpinにほぼ線形に依存するが，C および平均経路長 Lの非 線形性のために，small-worldness Sはpinに非線形に依存することが示された．従って， QとS の関係は実際はより複雑である．例えば，非常に密に接続されたネットワークで は，Qが増加するとSが減少することさえある（図2.6(a)）． また，実際のシステムに存在するモジュラーネットワークは，しばしばランダムではな いモジュール内接続を持つ場合がある．例えば，人間の皮質のさまざまな脳領域間の接続 性は，モジュール内に入れ子のモジュール構造が存在する階層的なモジュール性によって 特徴付けられる [12]．細胞ネットワークのモジュール内結合におけるべき乗則分布の重要 性も報告されている [29]．本研究では，ネットワークの次数分布は二項分布であると仮定 した．したがって，このような実際のネットワークの分析に直接使用することは出来な い．原則として，ここで示した平均場近似によるアプローチは，複雑なサブモジュール構 造を持つモジュラーネットワークの分析に拡張出来うるが，式が複雑になり，簡素な解析 により近似が可能であるといるメリットが縮小する． 他の考察すべき点として，平均場近似した結合行列で示されているネットワークは物理 空間に埋め込まれておらず，実際のネットワークが受ける物理的制約は考慮されていな い．距離依存性などの物理的制約を適用した場合，ダイナミクスから生じる偏差は平均場 近似 [30]から予測され，外部ノイズ [31]によってさらに影響を受けると報告した．この ような考慮事項は，得られた結果を物理ネットワークに適用する場合に注意が必要である ことを示唆している．よって，この研究で調査した平均場ネットワークに似た神経回路を 実現する場合には，距離の影響を最小限に抑える必要が生じる．この問題を回避する方法 として，細胞間の距離のばらつきを減らすために，3次元的にニューロンを成長させるこ となどが考えられる．

第

₃

章

モジュール構造型神経ネットワーク

における時空間神経ダイナミクスの

解析

3.1

本章の概要

本章では，培養神経回路の挙動を再現するような神経細胞モデルについて示し，第2章 で示したネットワークモデルを用いてネットワークを構成した際の，ネットワークダイナ ミクスについて解析する．3.3節では，ニューロンの動作とシミュレーションで用いたパ ラメータ群について示す．3.4節では，シミュレーションプログラムのコードの詳細や動 作フローについて説明する．3.5節では，ネットワークのダイナミクスを解析する為の手 法について説明する．3.6節では，モジュール構造を有する興奮性神経細胞からなるネッ トワークにおいて，モジュール化の度合いによりどのように活動が変調されるかを述べ る．3.7節では，興奮性神経細胞と抑制性神経細胞からなるネットワークのモジュール性 が，活動の複雑性にどのような影響を及ぼすかを定量・考察する．最後に，3.8節で本章 を総括する．

3.2

本章の目的

これまでに述べた通り，脳神経回路を含めた実世界の複雑ネットワークにおいてモ ジュール構造は普遍的な構造であり，その構造が生体において実現されることの優位性を 議論することは重要である．培養神経細胞を用いた実験において，培養神経細胞を用いて

図3.1 モジュール性を制御した培養神経回路の顕微鏡写真． モジュール（同一区画に配置された培養神経細胞が作るネットワーク）を構成し，そのモ ジュール間を疎に結合（区画間を細い連絡線で結ぶ）することによってモジュール構造を 実現し，その活動を調べる研究が行われている（図3.1）．この結線構造をシミュレーショ ン上で再現し，神経回路の持つモジュール構造がネットワークダイナミクスに及ぼす影響 を明らかにすることが狙いである．

3.3

ニューロンモデル

3.3.1

ニューロンモデルの比較

神経回路を構成する最小単位である神経細胞（ニューロン）は，基本的にその内外に 電位差が生じている．細胞外を基準としたときの細胞内との電位差を膜電位と呼び，平 衡状態では特に静止膜電位と呼ぶ．基本的に膜電位は負の値を持ち，他の細胞から信号 を受け取ることにより上昇して電位差が無くなる方向，すなわち細胞外の電位に近づく． 逆に信号が無い場合には静止膜電位に漸近する．膜電位がある値（閾値電圧）を超える と，多数のCaイオンチャネルが開口して膜電位は短時間の間に急激に上昇し，その後， Kイオンチャネルが開口することにより再度分極が起こる．この一連のインパルス状の 電位変化を発火と呼び，細胞の情報表現の大枠はこの発火により担われていると考えら れている．ニューロンの構造を機能毎に大別した場合には，主に三つの部分から構成さ れる．一つ目は他の細胞からの発火活動に由来する信号を受け取る樹状突起である．こ れは，ニューロンの入力に相当する．二つ目は細胞体であり，樹状突起で受け取った信 号を時空間的に加算する．三つ目は細胞で発生した発火を他の細胞に伝播する軸索であ る．これは，ニューロンの出力に相当する．樹状突起と軸索の間の結合はシナプスと呼ば れ，主に化学物質により信号の伝達が行われている．このシナプスの伝達強度が変調さ

れることで，多様な活動パターンが生み出される．ニューロンのダイナミクスを表現す る数理モデルは数多く提案されており，最も有名なものは1952年にHodgkinとHuxley が定式化したHodgkin-Huxleyモデルである [32]．Hodgkin-Huxleyモデルは細胞のイ オンチャネルの開閉の度合いによるコンダクタンスの時間変化を抵抗とキャパシタから なる等価回路で表現したモデルで，実細胞の物理的なメカニズムを忠実に考慮したモデ ルとなっている．しかし，このモデルは4本の連立微分方程式を解かねばならず，計算 量が大きいといった問題が存在する．一方で，正確な描像の一部を捨て，活動電位など の重要な挙動のみを再現することで計算量を削減したモデル [33]や，パラメータの物理 的対応を抽象化し，Hodgkin-Huxleyモデルの振る舞いをより簡単な関数形で再現した モデル[E. M. Izhikevich, 2003] など，以降多くのモデルが提唱されている．一般的に 計算の精密さは計算時間とトレードオフの関係にあるため，これらのモデルは調べよう としている現象の複雑さや求めている結果から適切に選択する必要がある．本研究では， 細胞の発火現象を再現するspiking-neuron の中でも最も簡単で，計算コストの少ない leaky-integrate-and-fireモデルを用いた．同モデルで用いられるパラメータは，French らが論じた培養神経細胞の挙動を再現したモデルに用いられているパラメータを基にして 決定した [18, 34].

3.3.2

ニューロンダイナミクスの計算

本研究で用いたニューロンモデルのパラメータは，培養神経細胞の生理学的なパラメー タと対応が可能なように調整・決定した．神経細胞のタイプは，発火により後細胞の膜電 位を増加させる興奮性ニューロンと，後細胞の膜電位を減少させる抑制性ニューロンの二 種類を用いた．興奮性ニューロンと抑制性ニューロンを示すインデックスをここではそ れぞれE,Iとする．興奮性ニューロンの膜電位VE(t) および抑制性ニューロンの膜電位 VI(t)の入力電流に対する時間変化はそれぞれ3.1,3.2式に基づいて計算する． τEdV E_(t) dt = V E L − VE(t) + REin(IE(t) + II(t) + IK(Ca)(t) + Iref(t) + ξ(t)), (3.1) τIdV I_(t) dt = V I L− V I (t) + RI_in(IE(t) + II(t) + 0.5ξ(t)), (3.2) ここで，τEとτIは膜電位の時定数であり，膜電位が静止膜電位まで減衰する速さを規定 している．V_LE と V_LI は静止膜電位でおり，入力が無い場合V (t)はこの値に漸近する． RE in と RinI は細胞の入力抵抗である. 膜電位が閾値電圧Vthを超えた場合，細胞が発火し

たとみなし，発火の後に膜電位はVresetにリセットする．シミュレーションではこれらの 微分方程式は差分方程式として計算しており，計算の時間ステップδtは0.1 msとした． IE(t) 興奮性のシナプス電流で，ニューロンiに流れる電流I_iE(t) は以下の式で計算さ れる． I_iE(t) = [E_RE− Vi(t)]gE ∑ j ∑ k wji [ exp (_{t− (t} j,k+ ∆ij) τE sd ) − exp(t− (tj,k+ ∆ij) τE sr )] , (3.3) ここで，E_REは興奮性シナプス電流の平衡電位を表しており，gEは興奮性シナプスコンダ クタンスである．τE sr とτsdE はそれぞれシナプス電流の増加と減衰に関わる時定数，wjiは ニューロンj, i間の荷重値，∆ij は伝達遅延，tj,k はニューロンiと結合しているニュー ロンj がk回目に発火した時刻である．指数項はt ≥ tj,k のみにおいて評価し，t < tj,k では0とする． II(t) は抑制性のシナプス電流で，ニューロンiに流れる電流I_iI(t)は以下の式で計算さ れる． I_iI(t) = [E_RI − Vi(t)]gI ∑ j ∑ k wji [ exp (_{t− (t} j,k+ ∆ij) τI )] , (3.4) 興奮性シナプス電流との主な違いは，興奮性シナプス電流では指数関数の差分である double exponentialとして計算したが，抑制性シナプス電流は単一の指数関数から計算し ている． IK(ca)(t) は 細胞内のCa2+ イオン濃度に依存した抑制性のK電流を表している．細胞 が発火するとCaイオンチャネルを通じて細胞のCa2+ イオン濃度が上昇する．その濃度 に依存してKチャネルの開閉度が変化し，抑制性のK+ イオン電流が流れることで膜電 位は低下する．これは，ニューロンにおいて時間的に連続した発火の周波数が徐々に低下 するadaptationの効果に対応する．ニューロンiに流れる Ca依存性K電流I_K(Ca)i (t) は以下の式で計算される．

I_K(Ca)i (t) = gK(Ca)ci(t)[EK− Vi(t)], (3.5)

dci(t) dt = cstep ∑ k δ(t− tk)− ci(t) τK(Ca) , (3.6) ここで，gK(Ca) はCa依存性K電流のコンダクタンスを表している．EK はKの平衡電 位を表しており，cstep は1回の発火により細胞内に流入するCaイオンの量，τK(Ca) は Caの流出に関する時定数である．tkはk 回目の発火が起こった時刻である．

Iref(t) はニューロンが発火した後，一定期間発火が抑制される相対不応期に対応する電 流であり，ニューロンiの相対不応期電流Ii ref(t)は以下の式で計算される． I_refi (t) =−gref ( 1 + t− tk τref )−1 Pref(t− tk)[Vi(t)− Vreset], (3.7) Pref(tk) = { 1 for tk < t < tk+1 0 otherwise (3.8) gref は不応期電流のコンダクタンス，τref は不応期電流の時定数である．なお相対不応期 とは別に，ニューロンが発火した後にいかなる電位変化も生じない期間である絶対不応期 tabsを設けている． ξ(t) は細胞の熱雑音などに起因するノイズ電流を表しており，電流の大きさは以下のα 関数で表される． ξi(t) = ANα(t− tNk; rN; τN), (3.9) α(s; r; τ ) = e −s/τ _{− e}−s/r e−¯s/τ − e−¯s/r, ¯s = rτ ln(r/τ ) r− τ , (3.10) ここで，ANはノイズ電流の最大値，tNk はk回目のノイズが発生した時刻，rNはα関数 の立ち上がり時間，τNはα関数の減衰の時定数に対応する．なお，ノイズの発生頻度は 平均をλ とする定常ポアソン分布に従うとした．これまでの式で定義した各パラメータ の代表的な値を表3.1に示した．

3.4

シミュレーションプログラム

3.4.1

シミュレーションプログラムの構造

シミュレーションプログラムは C++を用いて記述した．そこから出力されたデータ は，pythonを用いて可視化・解析を行った． ネットワークの構造解析に関して，複雑

ネットワーク解析のためにインディアナ大学のOlaf Spornsらが開発したMATLAB用 のツールボックスであるBrain Connectivity Toolbox [35]を導入し，Matlab 上で実行 した．

3.4.2

全体の動作フロー

表3.1 ニューロンモデルにおける各パラメータの代表的な値 記号 値 単位 τE 20 mV τI 10 mV V_LE −74 mV V_LI −70 mV RE in 40 MΩ RI_in 50 MΩ V_thE −54 mV V_thI −50 mV Vreset −60 mV EE R 0 mV E_RI −80 mV gE _{5 nS} gI 5 nS τ_srE 0.2 ms τ_sdE 5.3 ms τ_sI 6 ms ∆ij 2.8 ms gK(Ca) 10 nS/µM EK −75 mV cstep 0.1 µM τK(Ca) 2700 ms gref 159 nS τref 12 ms tabs 1 ms AN 1000 pA rN 30 ms τN 50 ms λ 0.5 Hz