平成２３年度における研究成果
戦略課題 3：予測医療に向けた階層統合シミュレーション[統括：高木 周（東京大学）]

Ⅲ

これまで個別に開発が進められてきた各種生体シミュレータ(血栓症、心臓、筋骨格、脳神経系等)を統合し、心筋梗塞やパーキンソン病等、様々な疾患に対してより複雑なプロセスを再現する。そのために、基盤ツールを整備するとともに、特定高速電子計算機施設を中核とする HPCIを活用することで病態予測や治療支援を目指す。

Ⅲ－１

高木 周（東京大学）
心筋梗塞シミュレーションに向けた血栓シミュレータと冠循環モデルの連成及び筋繊維の集合体としての骨格筋のシミュレーション

実施計画
本研究では、「予測医療に向けた階層統合シミュレーション」研究の一環として、「心筋梗塞シ

Ⅲ－１－１

ミュレーションに向けた血栓シミュレータと冠循環モデルの連成及び筋繊維の集合体としての骨格筋のシミュレーション」の研究開発を実施し、研究開発を統括する。
また、プロジェクトを推進する上で、関連する研究者と必要な協議等を行い、予測医療に向けた階層統合シミュレーション」の研究を統括する。
平成 23 年度は、心筋梗塞シミュレーションに向けた血栓シミュレータと冠循環モデルの連成として、現在進行中のグランドチャレンジプログラム「次世代生命体統合シミュレーションソフ
トウェアの開発」で開発中の血栓シミュレータと、心臓を取り巻く血管である冠動脈のシミュ
レーションを連成させて、フィブリン網の発達により血栓が凝固して血管を閉塞させるプロセスを再現する手法の開発を開始する。また、筋繊維の集合体としての骨格筋のシミュレーションとして、ヒトのふくらはぎ部分にあたる下腿三頭筋を対象として、個々の筋繊維の収縮により筋肉全体が収縮・変形する計算モデルを構築する。

1

Ⅲ－１－２

実施内容（成果）Ａ）心筋梗塞シミュレーションのための血栓症のモデリング
１）はじめに血栓が形成される過程は、血小板が粘着・凝集する一次血栓と、血液凝固機能の発動による二次血栓がある。
① 血小板の粘着・凝集(一次血栓/血小板血栓)
血管内壁の損傷部位のコラーゲンに vWF(von Willebrand Factor)が結合して、血小板が付着す
る土台をつくる。血小板は粘着系の糖蛋白である GP(glycoprotein) Ib と vWF との結合を介して、損傷部位に粘着をはじめる。この段階の結合は比較的弱く、可逆的で解離可能な状態であり一次凝集と呼ばれている。粘着した血小板は結合部位での刺激を受けて内部活性化し、別の粘着系の糖蛋白である GP IIb/IIIa を発現させる。GPII b/IIIa はフィブリノゲンと呼ばれる物質を架僑として血小板同士の接着を促進し、これにより、血管損傷部位への血小板の凝集が進む。この段
階では、非可逆的な強い結合状態であり二次凝集と呼ばれている。その際、血小板の活性化を促進する TX-A2(トロンボキサン A2)、セロトニン、ADP(アデノシン二リン酸)などの物質を放出して活性化を更に惹起する。
②血液凝固(二次血栓/フィブリン血栓)
次の段階で血液凝固の機能が働きはじめる。血管内皮組織や血漿、あるいは血小板から放出される、12 の血液凝固因子が多段の化学反応を経て、フィブリノゲンが重合してフィブリンが生成される。これが凝集した血小板を網状に覆い強固な血栓となっていく。
血栓形成においては、血小板の解離可能な一次凝集から非可逆的な二次凝集に進む過程が重要なポイントになるものと考えられる。光散乱法などによる血小板凝集能を測定した実験では、ADPが低濃度では一次凝集から二次凝集に至らず、ADP が高濃度になると二次凝集が起こるというよ
うな実験の知見があり、ADP は血小板凝集の過程において重要な役割を果たすものと考えられている。そのため、本年度は、ADP による血小板凝集の過程に着目し、ADP によって血小板が活性化される過程をモデル化することから着手した。
２）ADP 刺激による血小板活性化モデル
ADP は血漿などにも含まれているが、血小板が活性化した際に放出する濃染顆粒にも含まれており、従って血小板が活性化することに伴い、自らも ADP を放出して活性化を増幅することになる。
血小板の表面には 2 つの ADP 受容体があるが、その中の P2Y12 受容体に ADP が結合すると、その刺激を契機に血小板内では G 蛋白をはじめとする種々の化学反応が連鎖して、GP IIb/IIIa の活性化や濃染顆粒の放出などが誘発される。このような血小板内部での化学反応は活性化シグナ
リングと呼ばれている。P2Y12 受容体は ADP との結合刺激により、Gi(抑制性 GTP 蛋白質)によりAC(アデニル酸シクラーゼ)の活性化を抑制、cAMP(環状アデノシン 1 リン酸)濃度が抑制されることで Ca2+濃度が上昇することになり、結果として GP IIb/IIIa の活性化が誘導される。P2Y12 受
容体と ADP との結合は解離するため、P2Y12 受容体は断続的に ADP と結合する。このような血小板内部で生じる化学反応の連鎖について、反応速度論に基づき、常微分方程式の連立式としてのモデル化を試みた（図 1）。

2

図 1．反応速度論に基づいた ADP 刺激による血小板活性化モデル
ここで示した血小板内部での活性化シグナリングに関する化学反応については、関与する物質を含めて、まだ多くの部分が解明されていない状況であり、このまま評価を行うことは難しい。
ADP 刺激による血小板活性化モデルにおいて特に重要な点は、P2Y12 受容体が ADP と結合することにより血小板が活性化することで生じる GP IIb/IIIa 活性化と、ADP 放出(濃染顆粒)によるポ
ジティブ・フィードバックの機構であるため、この点に着目して次の簡易モデルで検討を行った。
・血小板表面に P2Y12 受容体と GP IIb/IIIa を配置し、ADP 濃度は各血小板単位で均一値とする。
・P2Y12 受容体が ADP との結合を契機に Ca2+の貯蔵量を一定量増加させる。
・P2Y12 受容体と ADP との結合は、当該血小板の ADP 濃度と所定の確率密度によって発生させる。
・P2Y12 受容体での ADP との結合は一定時間を経過すると解離する。
・血小板内部に貯蔵された Ca2+濃度が一定量を超えると、GP IIb/IIIa を活性化させるとともに、一定量の ADP を放出する。
３）ADP 刺激による血小板活性化の効果
GPII b/IIIa はフィブリノゲンを架僑として血小板同士が接着して二次凝集を促進する。つまり、GP IIb/IIIa が活性化する数は血小板の二次凝集の形成(大きさ)に直接的に寄与する指標の一つと考えることができる。ここでは、先に示した簡易モデルを用いたシミュレーションにより、
ADP 刺激で血小板が活性化されて、GP IIb/IIIa の活性化が増幅される効果について考察した。
血小板 100 個を正方格子状(x 方向 10 個×y 方向 10 個)に配置し、中央位置に一定値の ADP 初
期濃度を与えて拡散させる条件とした。実際の血小板表面に存在する P2Y12 受容体は 500～600
個程度、GP IIb/IIIa は 4 万～8 万個程度あり、これらの蛋白質の分布には局在性があることが示唆されている。ここでは各血小板について、500 個の P2Y12 受容体が 5 個単位で 100 か所に局在し、また 6 万個の GP IIb/IIIa が 200 個単位で 300 か所に局在しているものとして考え、一つ3

の局在箇所内の蛋白質(P2Y12 受容体は 5 個、GP IIb/IIIa は 200 個)は同じ動作をするものとして粗子化した。活性化した GP IIb/IIIa の数が多いことは、血小板凝集の大きさに相関するものと考えられる。血小板 100 個の GP IIb/IIIa の全体総数 30,000(=300 か所/血小板×100 血小板)の中で活性化した数の比率を GP IIb/IIIa 活性化率として、その経時変化をみることで、ADP による血小板凝集過程への影響について調べた。
①ADP 濃度による GP IIb/IIIa 活性化への影響
血小板が一次凝集から二次凝集に至る過程で、ADP 濃度が関係することが知られている。図 2
は、ADP の初期値を 5 つのパタン(0.1,0.3,0.4,0.5,1.0)で与えた場合について、GP IIb/IIIa の活性化率の経時変化を比較した結果例を示したものである。ADP 初期値が 0.4,0.5,1.0 の場合には、GP IIb/IIIa 活性化率は時間経過とともに増加傾向を示し、これは凝集が進展していく挙動として考えることができる。一方、ADP 初期値が 0.1,0.3 の場合には、GP IIb/IIIa 活性化率は、一旦は増加するが、時間経過とともに減少しており、これは凝集が進むことなく解離していく挙動として捉えることができる。今回の条件では、ADP 初期値が 0.3 と 0.4 との間を境に、凝集過程の違いが生じていることがわかる。

Rate of activated GP IIb/IIIa

0.8
0.7
0.6
0.5

ADP初期値=0.1

0.4

ADP初期値=0.3

0.3

ADP初期値=0.4
ADP初期値=0.50.2

ADP初期値=1.0

0.1
0
0

50

100

150

200

Time step
図 2．ADP 濃度による GP IIb/IIIa 活性化への影響②ADP 放出によるポジティブ・フィードバックの効果
図 2 で示した ADP 初期値が 0.4,0.5 の場合は、GP IIb/IIIa 活性化率は時間経過ともに増加しており、これは二次凝集に至る場合として捉えることができる。この 2 つのパタンについて、P2Y12受容体と ADP との結合により活性化した際に ADP 放出を行わない条件として行った計算とを比較することで、ADP 放出による GP IIb/IIIa 活性化へのポジティブ・フィードバックの効果について調べた。図 3 は、ADP 放出の有無による違いを比較した結果例である。ADP 放出を行わない場
合の結果を見ると、2 つのパタンの双方とも、時間経過の初期の段階では ADP 濃度が高いため
P2Y12 と ADP との結合が生じて GP IIb/IIIa の活性化率が増加するが、時間が経過していくとともに GP IIb/IIIa の活性化率が減少傾向を示している。これは凝集の成長が抑制されている挙動と考えることができる。血小板からの ADP 放出は、活性化のための ADP 濃度を維持して ADP 結合
による活性化と GP IIb/IIIa 活性化率のサイクルを増幅しているものと考えられる。このことか
4

ら、ADP 放出によるポジティブ・フィードバックは、血小板の二次凝集を促進することに大きく寄与する機構であることが示唆される。

Rate of activated GP IIb/IIIa

0.8
0.7
0.6
0.5

ADP初期値=0.4
ADP放出有0.4

ADP初期値=0.4
ADP放出無0.3

ADP初期値=0.5
ADP放出有0.2

ADP初期値=0.5
ADP放出無

0.1
0
0

50

100

150

200

Time step
図 3．ADP 放出による GP IIb/IIIa 活性化への影響
４）ADP 結合阻害による血小板活性化抑止の効果
抗血栓療法として、血栓形成の過程に対応して、血小板の粘着・凝集の過程を阻害する抗血小板療法と凝固系機能を阻害する抗凝固療法がある。動脈系の抗血栓療法においては、血小板血栓の抑止が有効であると考えられており、本研究課題では、抗血小板薬が血小板凝集を抑止する機能に着目する。抗血小板薬は、血小板の活性化を惹起する物質との反応を阻害することにより、血小板の活性化を抑え、結果として血小板凝集を抑止することを狙ったものである。代表的なものの一つにクロピドグレルがあり、臨床でその効果が確認されている。クロピドグレルは、P2Y12受容体と結合することで ADP との結合を阻害する。クロピドグレルと P2Y12 受容体との結合は永続的に保持されるため、クロピドグレルが結合した P2Y12 受容体からの活性化は一律抑制される
ので、GP IIb/IIIa の活性化や濃染顆粒の放出が抑えられることになる。ここでは、各血小板のP2Y12 受容体をクロピドグレルで一定の割合で塞いだ初期条件でシミュレーションを行うことにより、ADP との結合を阻害することによる GP IIb/IIIa 活性化抑止の効果について考察した。図
4 は、ADP 初期値 0.4 の場合(二次凝集に至ると考えられる場合)について、P2Y12 受容体をクロピ
ドグレルで 0%,30%,40%,50%の割合で塞いだ初期条件(ここでは CLO 閉塞率と呼ぶ)で計算した結果例を比較したものである。CLO 閉塞率が大きくなるにつれて、GP IIb/IIIa の活性化率が上昇する速さが抑制されている。今回の条件では、CLO 閉塞率が 40%と 50%との間で凝集過程の違いが生じており、CLO 閉塞率が 50%の場合には、時間経過とともに GP IIb/IIIa 活性化率が減少している。これは P2Y12 受容体と ADP との結合がクロピドグレルによるブロック効果により、ADP 結合による血小板活性化が抑制されて GP IIb/IIIa 活性化を減少させ、結果として二次凝集への進展が抑止されている挙動であることを示唆している。

5

Rate of activated GP IIb/IIIa

0.8
0.7
0.6
0.5

ADP初期値=0.4
CLO閉塞率=0%0.4

ADP初期値=0.4
CLO閉塞率=30%0.3

ADP初期値=0.4
CLO閉塞率=40%0.2

ADP初期値=0.4
CLO閉塞率=50%

0.1
0
0

50

100

150

200

Time step
図 4．ADP との結合阻害による GP IIb/IIIa 活性化抑止の効果
５）まとめ
本年度は、血小板が一次凝集から二次凝集に至る過程で ADP が重要な役割を果たすことに着目
し、ADP 刺激により血小板が活性化される過程についてのモデル化に着手した。P2Y12 受容体に
ADP が結合することによって生じる血小板の活性化において、ADP との結合、血小板内部での Ca2+濃度の上昇、GP IIb/IIIa の活性化、ADP の放出の連鎖反応についての簡易モデルを用いて、ADP
放出によるポジティブ・フィードバック機構が GP IIb/IIIa 活性化を増幅させる挙動についての確認を行った。今後、有限差分法に基づく流体構造連成計算手法を用いた血小板周りの流れ場の計算と連動して、ADP 刺激による血小板活性化と血小板凝集が形成される過程についての計算手法の構築につなげていく。
参考文献[1] 池田康夫、丸山征郎(編)、「血小板生物学」、メディカルレビュー社 (2004)[2] Alan D. Michelson, Platelets, Academic Press (2007)

6Ｂ）筋繊維の集合体としての骨格筋のシミュレーション
Ｂ－１） マルチスケール筋骨格系シミュレータの開発
１） はじめに筋骨格系のシミュレーションに関して、筋線維から筋肉全体までマルチスケール筋骨格系シミュレータの開発を引き続き行なった。具体的な対象としてヒトのふくらはぎ部を選び、腱や膜組織も含めてシミュレーションを実施した。平成 23 年度は特に、膜要素を新たに導入することで筋表層の膜組織のモデリングを行い、膜組織が筋の力発揮挙動に重要な役割を果たすことを明らかにした。さらに、パーキンソン病振戦の再現を目的として、剛体骨モデルおよびピ
ンジョイントモデルの導入を行い、筋収縮による関節運動を再現することが可能となった。
２）筋骨格系シミュレータの定式
骨格筋や腱などの軟組織は超弾性、非圧縮性、大変形、非線形性などの特性を有する。そこで本研究では、非圧縮性超弾性体の 3 次元混合型 FEM プログラムをベースにシミュレータの開発を行った(1), (2)。本手法の基礎となるひずみエネルギー関数は次式で与えられる(3)。

k1W  W (C )  ( J  1) 2  ( p  p ) 2
2
2k

(1)

ここで、右辺第 1 項は等容変形（ C  J  C ， C は右 Cauchy-Green 変形テンソル）によるひずみ
23エネルギーの寄与分で、超弾性材料構成式により記述される項である。第 2 項は体積変化（ J  det F ， F は変形こう配テンソル）によるひずみエネルギーの寄与分である。体積弾性率 kは非圧縮性を微圧縮性として導入したことによる体積変化を拘束するためのパラメータである。第 3 項は混合法の導入によって付加される項で、要素レベルで導入される圧力変数 p と変位場から計算される静水圧 p の不つり合いを補正する項である。
腱組織の材料モデルには次式の Mooney-Rivlin 体を用いた。
2

W tendon   cij  I1  3  I 2  3
i

j

(2)

i  j 1

ここで、cij は材料定数， I1 および I 2 はそれぞれ、右 Cauchy-Green 変形テンソルの第 1 および第(4)2 低減不変量である。骨格筋の材料モデルには次式の異方性モデルを用いた 。
(3)

W muscle  W matrix  W fiber

ここで、W matrix は筋内の結合組織による寄与分に相当し、(2)式の Mooney-Rivlin 体を仮定した。
pas
W fiber は筋線維による寄与分で、筋線維方向の Cauchy 応力成分  fiberと筋線維のストレッチ  の関係として次式与えられる。

W fiber paspas  fiber  0pasf fiber

0

(4)

pas
ここで、  0pas は材料定数、 0 は筋の至適長におけるストレッチ、 f fiberは次式の正規化した静止張力である。

0,
*pas
  1[e 2 (  1)  1],
f fiber0.4 2*( *  1.4)   1[e0.4 2  1],
 1 2 e7

for  *  1
for 1   *  1.4
for   1.4*(5)

ここで,  i は材料定数、  *   / 0 である。
(5)

骨格筋の筋収縮モデルには Hill-type model を採用した 。筋収縮によって生じる筋線維方向の Cauchy 応力成分  act
fiber は次式で与えられる。
isom
 act

fiber   0

 ft  fl  f v
0(6)

ここで、  0isom は筋線維の至適長（ストレッチ 0 ）における最大等尺性応力で材料定数である。
f t は筋の活動レベルを表す関数、 f l ,f v はそれぞれ筋長‐力関係(6)、収縮速度‐力関係(7)を表す関数で、次式で与えられる。

exp  b1 3f fl ( f )  2exp  b2 f 

for  f  0
for  f  0(7)

for  f  12.5

 0

 (1  k1 )k2  k
1  f  k2
c1
f v ( f )  c  f  c3 2  f  c3  2 k5c1
 c2

k5


for  12.5   f  0
for 0   f  8.556(8)for 8.556   f

k2 (1  c2 ) 2c, c2  k4 and c3  1(1  k1 )k3
1  c2*
ここで、 b1 , b2 および k1 , k2 , k3 , k4 , k5 は材料定数、 f    1 および  f は筋線維のひずみお
with

c1 

よびひずみ速度である。
３）筋表層の膜組織が骨格筋の力発揮挙動に及ぼす影響解析に用いる 3 次元モデルは、CT や MR などの医用画像から再構成される。本研究で用いた解析モデルは医用画像データをもとに、ヒフク筋、ヒラメ筋およびアキレス腱の領域抽出を行い、六面体要素に分割して作成した（Fig.1a）。医用画像の利用は、非侵襲的に生体から直接データを取得できるという利点がある。しかし、その取得手法の特性や解像度の問題から、画像データで認識できない組織や抽出が困難な組織も存在する。骨格筋では特に、筋表面を覆う
筋膜あるいは腱膜と呼ばれる薄い膜組織を画像データで識別することが困難となっている。そのため、計算モデルの再現性に問題があるのが現状である。本研究では、解剖学的な情報をもとに膜組織を再構成し、下腿三頭筋の等尺性収縮を対象に、膜組織の有無が筋の力発揮挙動（アキレス腱付着部位に発生する力）に及ぼす影響を調査した。
筋表面の膜組織は筋表面の有限要素節点を用いて膜要素（四角形要素）として作成し、その表面に固定した。膜組織は解剖学的な知見からコラーゲン線維が豊富である白色の膜組織の領域のみとした（Fig.1b, c）。この膜組織は腱膜（Aponeurosis）に相当する。
(8)

筋内の筋線維（筋束）方向は文献 を参考に、各有限要素積分点に単位ベクトルとして設定した。筋および腱の材料定数は文献(4), (9)より算出した。膜組織の材料モデルおよび材料定数は腱組織と同じであると仮定し、膜厚は死体解剖の結果より 0.5mm とした。
骨と付着する部位の節点の自由度を固定し、ヒフク筋およびヒラメ筋に 100％MVC (Maximum
voluntary contraction)の筋活動を与えることで、等尺性収縮の解析を行った。
8

Fig. 1 FE model of the human triceps surae muscle, reconstructed from MR images. (a) Finite element
mesh segmented into Achilles tendon, gastrocnemius and soleus muscles. Achilles tendon andaponeurosis (membrane element) on the surfaces of (b) gastrocnemius muscle (green) and (c) soleus
muscle (blue).
下腿三頭筋の等尺性収縮シミュレーションの結果を Fig.2 に示す。筋収縮により、アキレス腱はヒフク筋およびヒラメ筋に挟まれた腱膜（深部腱膜）を通して引っ張られ、アキレス腱付着部
位に力が発生する（Fig.2a）。アキレス腱付着部位に発生した長手方向の力は、膜組織なしの場
合には 645.6N、ありの場合には 1080.8N であった。また、同様な収縮形態での実験では 1354.0Nの力発揮であった。膜組織によって、下腿三頭筋が発揮する力は膜なしに比べて 1.67 倍に増加
し、実験値に近づいた。したがって、筋表層の膜組織が筋の力発揮挙動に大きな影響を及ぼすことが示された。
膜組織の存在による力の増加の要因の一つは、膜組織を通して、筋表面の筋線維端と骨が結合されることによるものと考えられる。膜組織がある場合のヒフク筋表面の下向きの変位量は膜組織なしの場合に比べて明らかに小さくなった（Fig.2b,c）。これは膜組織の存在により、ヒフク
筋表面側にある筋線維端が大腿骨と接続され、筋線維の自由な収縮が拘束されためと考えられる。
表面の筋線維端が膜組織によって拘束されることにより、筋収縮による力はその他端である深部

Fig.2 Simulation results of the isometric contraction of the human triceps surae muscle. (a) The nodal
force vectors at the end of the activation. The circle indicates the forces at the insertion with thecalcaneus. Distributions of the displacement in the longitudinal direction at the end of the activation(b) without and (c) with membrane.

9腱膜の変位の増加を引き起こし、その結果、アキレス腱付着部に発生する力が増加したと考えられる。
４）3 次元 FEM による関節運動のシミュレーション
本研究ではこれまで、3 次元 FEM による筋骨格系シミュレータの開発を行ってきた(1), (2)。本シ
ミュレータは生体軟組織の変形と骨格筋の筋活動挙動の解析を行うものである。現在、本シミュ
レータとその筋活動の入力を司る脳神経のモデルを組み合わせることで、全身筋骨格-神経系統
合シミュレーションによる、脳神経疾患の運動機能障害の病態予測と治療支援を目指している。
脳神経疾患のひとつにパーキンソン病が挙げられる。パーキンソン病は大脳基底核内の神経細胞の変性によるドーパミン不足がひとつの原因と考えられており、手足のふるえ（振戦）、筋固
縮、動作緩慢、歩行障害などの症状を示す。本研究ではその第一ステップとして、パーキンソン病振戦の再現を目指している。そこで、本研究では本筋骨格系シミュレータを拡張し、剛体骨モ
デルとピンジョイントモデルによる関節運動のモデリングを行った。
骨は骨格筋に比べて非常に大きな剛性を持つため、その変形は無視できると仮定し、剛体として取り扱った。ここでは剛体と変形体を結合させる方法として、Lagrange 未定定数法を用いた。
変形後の剛体の剛体回転を R とするとき、剛体と変形体の接合は次式で与えられる。

x  xg  R  ( X  X g )(9)

ここで、 X g および X はそれぞれ、初期配置における剛体の重心および接合する変形体節点の座
標、 x g および x はそれぞれ、変形後の剛体の重心および接合する変形体節点の座標である。
Lagrange 未定定数 λ を導入し、次式をエネルギー関数(1)式に付加することで、新たな付帯条
件付き停留問題に変換される。

  λ  { x  x g  R  ( X  X g )}(10)(10)式の変分を求めることで、次式の内力ベクトルを得る。

{Qu }  {}{Q }  {}
 U

T
T
{Q m }  { X  X g } {R} {}{Q m }  {x}  {x }  [ R]{ X  X }gg 
さらに、 [ K RE ] 

(11)

{Q}
より、次式の剛性行列を得る。
{x}

0 0
0 0

[ K RE ]  
0 0 I  I



0
I

T
T
T
T
{ X  X g } R {} { X  X g } [ R ] [ R ]{ X  X g }
0

I

0



(12)

(12)式は剛体と変形体の接合により、新たに付加される剛性行列である。同様に、剛体同士の結合も Lagrange 未定定数法を用いて定式を行った。
Fig.3 にヒト足関節の関節運動のシミュレーション例を示す。足関節は剛体 1 と剛体 2 がピンジョイントにより結合されている(Fig.3a)。下腿三頭筋のアキレス腱は踵骨で骨と結合され、下腿三頭筋の筋収縮により、足関節の運動を与えた。筋活動にしたがって、ヒフク筋およびヒラメ
10Fig. 3 Simulation results of the ankle joint movement generated by the activation of the triceps suraemuscle. (a) Ankle joint model consist of rigid bodies and pin joint. (b) Ankle joint movement during
the activation of the triceps surae muscle.
筋が筋収縮し、その結果、アキレス腱が上方へ引っ張られ、足関節まわりの運動を生じた(Fig.3b)。
５）まとめ筋骨格系のシミュレーションに関して、筋線維から筋肉全体までマルチスケール筋骨格系シミュレータの開発を引き続き行なった。具体的な対象としてヒトのふくらはぎ部を選び、腱や膜組織も含めてシミュレーションを実施した。3 次元 FEM によるヒト下腿三頭筋の等尺性収縮の解析において、筋表層の膜組織がアキレス腱付着部位に発生する力へ及ぼす影響を検討した。そして、膜組織がある場合に筋が発生する力は、膜組織がない場合に比べて 1.67 倍増加するという結果を得た。その結果、筋表層の膜組織が筋の力発揮挙動に重要な役割を果たすことを明らかとなった。
さらに、パーキンソン病振戦の再現を目的として、剛体骨モデルおよびピンジョイントモデルの導入を行い、筋収縮による関節運動の再現が可能となった。
今後、全身筋骨格-神経系統合シミュレーションに向けた脳神経系との連携、また、筋活動の不均一性やカルシウムイオン感度の導入など、よりミクロスコピックな筋収縮モデルの導入を行う予定である。
参考文献
1) Luis ALVES, Naoto YAMAMURA, Toshiaki ODA and Cristian TEODOSIU, Proc. of CMBBE 2010,
Valencia, Spain.
2) 山村 直人, Luis ALVES, 小田俊明, Cristian TEODOSIU，機械学会第 23 回バイオエンジニリ
ング講演会講演論文集，pp. 325-326(2011.1)．

3) Sussman, T., Bathe, K.-J., 1987, Computer & Structures 26, pp. 357-409.
4) Chi, S.W., Hodgson, J., Chen, J.S., Edgerton, V.R., shin, D.D., Roiz, R.A., Sinha, S.,
2010, Journal of Biomechanics 43, pp. 1243-1250.
5) Johansson, T., Meier, P., Blickhan, R., 2000, Journal of Theoretical Biology 206, pp.
131-149.
6) Jenkyn, T.R., Koopman, B., Huijing, P., Lieber, R.L., Kaufman, K.R., 2002, Physics in
Medicine and Biology 47, pp. 4043-4061.
7) Nagano, A., Gerritsen, K.G.M., 2001, Journal of Applied Biomechanics 17, pp. 113–128.
118) Kawakami, Y., Ichinose, Y., Fukunaga, T., 1998, Journal of Applied Physiology 85, pp.
398-404.
9) Louis-Ugbo, J., Leeson, B., Hutton, W.C., 2004, Clinical Anatomy 17, pp. 30-35.

12

Ｂ－２）

パーキンソン病の運動神経ネットワークモデルの開発

１）はじめにパーキンソン病のような神経変性疾患は日本ではその多くが特定疾患に指定されており，根本的な治療法も確立されていない．本研究では，スーパーコンピュータ「京」上に疾患状態の人体を再現し，健常者モデルとの比較等を通じて病因究明や新たな治療法の開発を目指す．
人体における運動は，大脳基底核から発せられた運動指令（電気的なスパイク信号列から成る）
が，神経ネットワーク内を修飾・伝達され，筋肉の収縮を引き起こすことで実現されている．そこで，神経変性疾患に伴う運動機能障害を再現するモデルを構築するにあたり，人体モデルを，脳モ
デル・神経ネットワークモデル・筋骨格系モデルの 3 つに大別し，次年度以後相互接続による統合を図ることを前提とした上で個々のモデルの開発を進めることとしている．本報告では，特に神経ネットワークモデルの開発について述べる．
運動神経回路は，脳から発せられた運動指令信号を筋肉に伝える役割を担うが，信号を脳から筋肉へと伝達するだけでなく，筋肉の状態（長さ変化や発揮している力）をモニターし，適切な運動を起こす，あるいは，適切な力発揮を持続するための信号修飾の役割も持つ．そこには数多くの神経細胞が信号を相互に伝達し合う複雑なネットワーク回路が存在する．脳から発せられた信号が筋肉の収縮へと至る過程を正確に再現する上では，このネットワーク回路モデルの構築も重要な要素の一つであり，脳モデルから得られる信号を受ける，また，筋骨格系モデルへと信号を送出するな
どの接続性を考慮することが必要である．本報告では，今年度開発した簡易的な運動神経ネット
ワークモデルについてその概要と結果を述べるとともに，詳細神経ネットワークモデル開発に関する取り組みについて報告する．
２） 巨視的モデルによる運動神経ネットワークの開発とドーパミン枯渇の影響に関する検討
脳を含む上位中枢神経系で作られた運動指令は，下位中枢神経系を経由して筋肉の運動を引き起こ
す．下位中枢神経系の中でも筋線維の収縮を直接支配する神経細胞は -運動ニューロンである．運動の最小単位は筋線維であるが，一つの -運動ニューロンが支配する筋線維は通常数本から数千本程度（筋肉によって異なる）であり，一つの筋肉の運動には数百から数千程度の -運動ニューロンの活動が関係する．同時に，同程度の数のフィードバック系の神経細胞や介在ニューロンなどが関係し，さらに一つの関節運動には少なくとも 2 つの筋肉が関係することも考慮すると，一つの関節運動に関係する神経細胞はFig. 1 Neural network of spinal cord (I: Ia fiber, X: Ib fiber, R: Renshaw cell, M: -motoneuron, S: static

-motoneuron, D: dyanmic -motoneuron, blue dot line: inbitatory signal, red arrow line: excitatory signal,
black arrow line: excitatory signal from cortex)
13

おおよそ数万から数十万程度になると見積もられる．
本研究を進めるにあたり，第一ステップとして脳から関節運動までを巨視的にモデル化した数理モデルを構築した(1)．神経ネットワークモデルとして，個々の神経細胞の活動とそれらの間での信号伝達を解析するのではなく，神経細胞の集合体としての活動度と，集合体間での信号授受をモデル化し，その状態に基づき関節運動を記述するモデルとなっている．Fig.1 にここで考慮した下位中枢神経系における神経細胞とそのネットワークの概略図を示す．このネットワークモデルには，-運動ニューロン，静的・
動的 -運動ニューロン，Ia 線維，Ib 線維，Renshaw 細胞が考慮されており，筋活動に必要な主要な神経細胞がフィードバック系も含めてモデル化されている．主動筋・拮抗筋それぞれに対してこれらの神経細胞を考慮し，また解剖生理学から得られた結合を基にモデルネットワーク回路を構築している．本来一つの筋肉に対してこれらそれぞれの細胞が数千程度ずつ関係するが，ここではそれらを一つの集合体として考え，その活動度をモデル化して数値的に解析する．例えば，-運動ニューロンの活動度は次のような常微分方程式によって記述する．

  

dM i  Bi  M i   DA7   Ai  P  Ei  M i  2   DA8  1  Ri    X i   I j
dt
dM j
dt



    X   I    B j  M j  DA8  A j  P  E j  M j  2 DA7  1   R j





j



i

ここで，M，R，X，I は Fig.1 内と同じ各神経細胞群の活動度を表す．下付添字 i，j は主動筋・拮抗筋を各々意味する．B，A，P，E，， は脳を含む上位中枢神経系から来る信号やそれに関連したモデル定数を表している．ここでは，大脳基底核から発せられる信号に，運動の最終目的位置と現在の位置との差からなる情報を合成したものを数値的にモデル化しており，脳科学的には，大脳基底核と運動野の働きをモデル化したものとなっている．，， はモデル定数である．Fig. 1 と見比べるとわかるように，主動筋の-運動ニューロンは，上位中枢神経系から興奮性信号を受けており，また，主動筋側の Renshaw 細胞，
主動筋側の Ia 線維および拮抗筋側の Ib 線維から抑制信号を受けていることがわかる．これらの合成として筋収縮を引き起こす興奮性信号を筋線維に送出している．このような数値モデルを各細胞について連立させた連立常微分方程式を解析することで，各神経細胞群の活動度を得ることができ，それに応じて肘関節の運動を検討することができる．
ここで構築した運動神経ネットワークモデルを用い，肘をある目的位置まで曲げる運動におけるドーパミ

Fig. 2 Input signal from basal ganglia (left). Experimental GPi post-stimulus time histogram in MPTP
monkey(2) (right).

ン枯渇による信号伝達の影響を検討した．ドーパミン枯渇による信号伝達の影響は，伝達される信号強度を簡易的に弱くすることで模擬している．例えば，正常時には R → M へと 1 の信号が伝達されるとき，

14ドーパミン枯渇によりこれが 0.7 に減少するとしてモデル化している．ここでは，(1) ドーパミン枯
渇がない場合（正常な中枢神経系の場合），(2) 上位中枢神経系においてのみドーパミン枯渇がある場合，(3) 上位・下位中枢神経系においてドーパミン枯渇がある場合，の 3 通りを検討した．大脳基底核から発せられる信号として Fig.2 を用いた．これは，薬物投与によりパーキンソン病状態にされた猿の随意運動時における信号波形(2)に基づいている．
Fig.3 にこのときに得られた運動履歴を示す．この図からわかるように，(1)と(3) の条件では，スムーズな運動によりほぼ目的位置に達した後に運動が停止しており，運動としてはほぼ正常な結果が得られている．これに対し，(2) の場合には最初の運動に続いて振動が見られており，最終目的位置に正常に達することができていない．このときに主動筋・拮抗筋それぞれが発揮していた力をFig.4 に，また，-運動ニューロン群の活動度の径時変化を Fig.5 にそれぞれ示す．正常時の (1) の条件では，初期に見られる大きな主動筋の活動によって運動が実現されており，拮抗筋は力を発揮していない．対応するように主動筋側 -運動ニューロン群の活動が初期に大きく見られる．これに対し(2) の場合，主動筋の初期の力発揮に続き，連続的な力発揮が現れている．-運動ニューロン群の活動
15

も最初の大きな活動に続き，複数回の活動が見られる．さらに特徴的な結果として，主動筋が力を発揮しているときに，同時に拮抗筋も力を発揮してしまっており，これは，-運動ニューロンの活動にも表れている．
これが，Fig.3 で見られた振動を伴う運動履歴の原因となっている．本来正常時の運動では，主動筋が収縮するときには拮抗筋は弛緩し，これによって関節運動が実現されるが，(2) の条件では両者が同時に収縮しており正常な運動が実現されていない．このような状態はパーキンソン病時の固縮に類似したものと考えられる．(3) の場合，運動履歴からは正常時の(1) と顕著な差は見られないものの，発揮している力は (1) に比べて弱く，また，-運動ニューロンの活動も複数回の連続的な活動が見られており，この点は (2) の状況に類似した結果となっている．
本モデルは，大脳基底核からの運動指令がパーキンソン病状態に見られるような非正常な信号であっても，ドーパミン枯渇部位の違いにより運動状態が異なることを示している．大脳基底核から発せられる信号が異常でも，上位・下位中枢神経系が正常であれば運動はほぼ正常に行われるが，大脳基底核の異常に加え，上位中枢神経系の異常が伴うとパーキンソン病時の固縮に類似した筋活動が見られることを示している．
３） 詳細神経ネットワークモデルの開発と大規模化への取り組み
先に述べた巨視的モデルは，簡易的なモデルでありながら運動時の神経細胞の活動状態やそのときの相互作用を検討することができ，第一ステップとしては優れたモデルである．しかしながら，モデル定数を多く含むため将来的な拡張性に乏しく，また，各神経細胞の活動度が厳密な意味での物理量と必ずしも対応しない等，脳モデル・筋骨格モデルとの統合という観点からは必ずしも優れているとは言えない面も持つ．そこで次のステップとして，Hodgkin-Huxley モデルのようなイオ
ンチャネルの開閉による膜電位変化を記述する数理モデルに基づき，運動神経細胞ネットワーク間での信号伝達を考慮可能な数値解析モデルの開発に現在取り組んでいる(3)(4)(5)．このようなモデル
では，運動神経ネットワーク内で伝達される信号は，物理的に意味のある膜電位であるために，脳モデルから得られる運動指令との接続に関しても都合がよい．また，巨視的モデルとは異なり個々の神経細胞の状態（膜電位）を直接解析するため，巨視的モデルに比べより詳細な筋線維の分布等と接続することができ，今後の拡張性も高いことが期待される．巨視的モデルに比べて不利な点としては，上で触れたように数千から数十万程度の細胞を個別に解くために計算負荷が大きくなることが挙げられるが，Fig. 6 Schematic view of detailed mathematical model of neural network.

16

「京」の計算能力を考慮すればそれほど大きな問題とはならないと考えられる．Fig. 6 に詳細モデルの概念図を示す．考慮する神経細胞種としては，巨視的モデルとほぼ同じでありフィードバッ
ク系も含む．大きな違いは，各運動神経細胞種についてそれぞれ数千個程度を数値的に扱い，それぞれの電位変化を個別に解くことである．そのとき，個々の運動ニューロンの発火閾値を連続的に変えることで，速筋・遅筋などの筋肉種に応じた連続した力発揮（recruitment）を考慮することも可能であり，三次元筋骨格系モデルとの統合という観点からも都合がよい．
詳細神経細胞ネットワークモデルの開発と並行して，筋骨格系モデルの大規模・並列化に向けた
コードの改修も行っており，今後はこれらを接続する数理モデルを開発することが当面の目的となる．
４） まとめ神経変性疾患に伴う運動機能障害の原因究明や新たな治療法の開発などを目指し，スーパーコン
ピュータ「京」において疾患状態の運動を解析可能なシミュレーションモデルの開発に取り組む中で，本稿では特に神経ネットワークモデルの開発に関して報告した．
第一ステップとして，中枢神経系の活動状態から関節運動までを記述する巨視的モデルを開発し，肘を曲げる運動を行うときの神経細胞の活動状態とドーパミン枯渇による信号伝達の異常に関して検討した．
大脳基底核から発せられる運動指令に異常がある場合においても，脳内の信号伝達が正常であれば，関節運動に大きな異常は見られないこと，大脳基底核から発せられる信号が異常であり，なおかつ脳内における信号伝達にも異常があるとパーキンソン病の固縮に類似した症状が見られることを示した．また，脳モ
デル・筋骨格モデルとの統合や「京」の計算能力を最大限発揮することを考慮して，詳細な運動神経ネット
ワークモデルの開発に取り組んでいる．今後は詳細モデルの開発をさらに進め，筋骨格モデルとの統合を図る数理モデルの構築に取り組む予定である．
参考文献
1)

Vassili Cutsuridis, Origins of a repetitive and co-contractive biphasic pattern of muscle
activation in Parkinson's disease, Neural Networks, 2011, pp.1-10.

2)LEon Tremblay, Michel Filion and Paul J. B6dard, Responses of pallidal neurons to striatal
stimulation in monkeys with MPTP-induced parkinsonism, Brain Research, 498(1989),
pp.17-33.

3)

Rogerio R. L. Cisi & André F. Kohn, Simulation system of spinal cord motor nuclei andassociated nerves and muscles, in a Web-based architecture, J. Comput. Neurosci, 25(2008),
pp.520-542.

4)

Madeleine M. Lowery and Zeynep Erim, A Simulation Study to Examine the Effect of Common
Motoneuron Inputs on Correlated Patterns of Motor Unit Discharge, J. Comput. Neurosci.,
19(2005), pp.107-124.

5)

Arno H. A. Stienen, Alfred C. Schouten, Jasper Schuurmans, Frans C. T. van der Helm,
Analysis of reflex modulation with a biologically realistic neural network, J. Comput.
Neurosci., 23(2007), pp. 333-348.

17

Ⅲ－２

中村 仁彦（東京大学）
全身筋骨格モデルにおける筋肉の体積効果及び人の行動モデルに関する大規模並列計算実施計画
本研究では、「予測医療に向けた階層統合シミュレーション」研究の一環として、「全身筋骨格モ

Ⅲ－２－１

デルにおける筋肉の体積効果の大規模並列化及び人の行動モデルに関する大規模計算」の研究開発を実施する。
また、「全身筋骨格モデルにおける筋肉の体積効果及び人の行動モデルに関する大規模並列計算」の研究を行う上で、関連する研究者と必要な協議等を行うとともに、本格実施に必要な研究体制の整備を行う。
平成 23 年度は、
（１）東京大学・高木グループの「筋繊維レベルからの骨格筋モデル」、沖縄科学技術研究基盤整備機構・銅谷グループの「脳神経系の多階層モデル」に対して東京大学・中村グ
ループの「脊髄と四肢筋骨格モデル」を統合するためのインターフェイスを定義する。
（２）体積効果を考慮した四肢筋骨格モデルの開発を開始する。
（３）人の行動モデルの大規模並列処理の基
礎アルゴリズム開発を開始する。

Ⅲ－２－２ 実施内容（成果）（１）モデル統合インターフェイス「脳神経系の多階層モデル」、「脊髄と四肢筋骨格モデル」、「筋繊維レベルからの骨格筋モデル」
を統合したシステムの目標として、当面はパーキンソン病の症状として現れる運動障害の発生機序を再現することを目指すこととした。このために大脳基底核のドーパミン受容体の変性から始まる運動障害のモデルについて Cutsuridis et al. [R1] などの文献を調査し、統合モデルのインターフェイスの具体化を行ってきた。
東京大学の高木チーム、沖縄科学技術研究基盤整備機構の銅谷グループと議論を行いなが、われ
われが開発してきた筋骨格シミュレータと中枢神経系シミュレータをつなぐインターフェイスを定義する作業を行ってきた。統合モデルのインターフェイス開発を行うソフトウェア環境として
NEST (Neural Simulation Technology)[R2]、および大規模なプロセス間の通信を可能にする環境として MUSIC [R3]を利用する方針とした。また、骨格筋モデルの有限要素法計算と全身骨格モデルの剛体力学シミュレータの開発環境は、前者が Fortran を後者が C++を基盤としているため、両環境において相互に参照可能である。これらをつなぐインターフェイスは、剛体のマクロな力学情報
および筋有限要素全体の境界条件のやり取りを行い、通信情報量の低減による並列計算効率の向上を目標とする。
筋骨格モデルの入力信号は連続値として表された筋活動度である。これに対して、一般的に神経活動は離散的なスパイク信号として表されることが多い。スパイク信号から連続波形への変換は比較的に容易であるが、筋の活動情報を中枢神経系へ導く際に必要になるその逆の変換は単純ではない。このような変換の方法について研究開発を行ってきた。

（２）体積効果を考慮した四肢筋骨格モデルの開発
筋肉の体積効果を考慮した全身の筋張力の並列計算手法を開発した（図 1）。本手法では、ロボティ
クスにおける剛体リンク系の力学計算法と、筋形状モデルに対する非線形有限要素解析を同時に解く。筋繊維の三次元的な配向や骨接触部の面圧力などの詳細な力分布を再現できる。筋間の接触力
やそれに付随する筋の変形も将来的には統合可能である。

18

具体的な計算方法は、有限要素解析から求まる各筋の有限要素の非線形等式拘束条件の下で、運動計測から求まる剛体リンク系の関節トルクを最小化するような各筋の筋活動度を求める非線形最適化問題となる。逐次二次計画法を応用することで、一回あたりの収束計算の最中では、筋の有限要素計算を剛体リンク系の力学計算と独立に分離することができるため、一つ一つの筋の有限要素計算を並列に計算することができる。また一つの筋の有限要素計算単体において、有限要素解析の並列計算を利用することも可能である。
計算論的神経モデルから人間の身体の物理モデルへの結び付けるための人間の解剖学モデルの有用性を検討し、Zygote 社 [R4]の人間の全身のポリゴンモデルを採用した。Zygote モデルは人間の幾何形状情報として、詳細な骨格、筋、などの運動器の他に、神経系、循環器系、消化器系などを含んでいるものの、機能情報を持たない。このため骨格上の筋腱複合体の端点、求心性・遠心性の運動神経と筋との接合点などの機能情報を解剖学の教科書を見ながら追加する作業を行ってきた。
Zygote 社製の三次元人体形状モデルを基盤に、下肢の骨格モデルおよび膝関節回りの五つの主要な筋の有限要素モデルを作成した（図 2）。簡易モデルでは筋数が 5 つのため、4 コア CPU 一つを利用して 5 スレッドの並列計算を行った。作成したモデルを用い、筋活動推定の際の二次計画問題の評価関数の導出一回を並列計算と逐次計算の二種類で実行し実行時間を比較した。逐次計算、並列計算での実行時間はそれぞれ 1.275×104 秒、5.058×103 秒であり、並列計算により逐次計算の約2.4 倍の実行速度が得られた。また得られた 5 つの筋活動度は、解析対象とする運動中の膝関節の関節トルクを再現できた。並列性が高い計算アルゴリズムであることが示唆され、より詳細な全身筋骨格モデルにおける高速計算に大きく貢献できると期待される。

図 1．体積効果を考慮した筋活動度推定法

19図 2．下肢の骨格および有限要素筋モデル（３）人の行動モデルの大規模並列処理の基礎アルゴリズム開発
人間の多様な身体運動パターンを統計モデルの一種である隠れマルコフモデル(Hidden Markov
Model : HMM)を用いて学習・分類することによって人間の行動モデルを構築する枠組みを構築してきた[R5]。
運動の時系列データに対して、時間方向にわたって各 HMM のノード間の遷移確率と各ノードからの運動データの出力確率の積の積算として、運動データが HMM から生成される尤度が計算される。この尤度に基づいて運動を分類することが可能となる。しかし、行動モデルの中で扱うべき運動パターンの数や運動の詳細化による運動データの次元数の増大によって、運動の分類に要する計算時間が増加する問題がある[R6]。新たに HMM から運動データが生成される確率の計算を①HMM の各ノードから各時刻の運動データが生成される確率の計算②各時刻の運動データが各ノードから生成される確率の積算の 2 つに分け、各確率計算を並列化するアルゴリズムを開発した。①では、各ノード毎の生成確率計算を 1 つのスレッドとして、②では各 HMM 毎の生成確率計算を 1 つのスレッドとした。開発した並列処理の試験として、従来の計算法による尤度計算との比較を GPU マシン(NVIDIA TESLA C2075、
Clock Rate 1.15GHz)を用いて検証した。ノード数 30 個の HMM に関して、従来の尤度計算より 47.9倍の計算速度を実現することを確認した。この並列化アルゴリズムは、行動の分類のみならず学習における収束計算においても利用される計算であり、大規模な運動データの学習・分類計算の高速化に寄与することができる。
参考文献[R1] V. Cutsuridis. Origins of a repetitive and co-contractive biphasic pattern of muscle
activation in Parkinson’s disease. Neural Networks, 24(6):592-601, 2011.
[R2] M-O. Gewaltig, M. Diesmann. NEST (Neural Simulation Tool), Scholarpedia, 2(4): 1430,
2007.
[R3] M. Djurfeldt, J. Hjorth, J.M. Eppler, N. Dudani, M. Helias, T.C. Potjans, U.S. Bhalla,
M. Diesmann, J.H. Kotaleski, O. Ekeberg. Run-time interoperability between neuronalsimulators based on the MUSIC framework. Neuroinformatics, 8(1): 43-60, 2010.
[R4] Human 3D Geometric Model, Zygote Media Group, www.3dscience.com[R5] T. Inamura, I. Toshima, H. Tanie, Y. Nakamura, Embodied symbol emergence based on mimesis
theory, International Journal of Robotics Research, Vol.23, No.4, pp.363-377, 2004.
[R6] W. Takano, H. Imagawa, Y. Nakamura, Prediction of human behaviors in the future throughsymbolic inference, IEEE International Conference on Robotics and Automation,
pp.1970-1975, 2011.

20

Ⅲ－３

銅谷賢治（沖縄科学技術大学院大学）
行動制御の脳神経系の多階層モデルフレームワークの構築

実施計画「予測医療に向けた階層統合シミュレーション」の研究では、循環器系、筋骨格系、脳神経系のⅢ－３－１

数値モデルを構築し、それらをつなぐ大規模シミュレーションにより、心筋梗塞、パーキンソン病などの発症プロセスを解明し、治療手法の探索評価を可能にすることをめざしている。その一環として、脳神経系の感覚入力から運動出力までをつなぎ、脳全体の回路から細胞、分子レベルのダイ
ナミクスを含む多階層の脳神経系モデルの構築が不可欠である。
平成 23 年度はそのようなモデル構築のために使える既存のデータ、モデル、計算手法、ツールのサーベイをもとに、モデルの全体の最適なデザインを策定し、その実現に必要な研究開発項目を明らかにし、モデル構築と計算手法開発に着手する。
Ⅲ－３－２

実施内容（成果）（１） 研究開発の実施状況
１）多階層の脳神経系モデル構築のための既存のデータ、モデルとツールのサーベイ
大脳基底核やパーキンソン病に関連する文献のサーベイを行った。既存の大脳基底核のニューロン数やコネクションのパターン、電気生理実験等の生理学的な知見、大脳基底核の機能とパーキンソン病の発生のモデルまでを幅広く調べた。
またシミュレータのサーベイでは、NEST や Neuron をはじめとした神経ネットワークシミュ
レーターから GENESIS や PhysioDesigner 等の細胞内分子カスケードシミュレーションに向いたシミュレータまでの広範なサーベイを行った。また、モデルを記述する上で重要となる XML 系言語（NeuroML、NineML、PHML）とシミュレータの関係性についてもサーベイを行った。
また、最新のモデリング手法の修得のため、以下の学会やサマースクールに参加した。
・ OCNC 2011(2011/6/13-2011/6/30)
・ Berkeley course in mining and modeling of neuroscience data (2011/07/11-2011/07/25)
・ Neuroinformatics training course & Neuroinformatics 2011 (2011/09/01-2011/09/06)
・ Allen Institute Course in Molecular Neuroanatomy (2011/10/10-2011/10/22)２）モデル全体の最適なデザインの策定
大脳基底核の神経回路モデルとしては、スパイキングニューロンモデルを用いた庄野らのモデル（文献１）を基盤として、生理学的に精緻な大規模モデルを開発することに決定した。小脳の神経回路モデルとしては、マルチコンパートメントタイプのニューロンモデルを用いた Erik De
Schutter らの小脳モデルを用いることに決定した。また大脳基底核と小脳の両方とループを形成する大脳皮質の神経回路モデルとしては、現在グランドチャレンジで開発が進むスパイキング
ニューロンモデルを用いた深井、Diesmann らの大脳皮質モデルを基盤にする方向で考えている。
これらの異なる脳部位のモデル間におけるスパイク情報の交換には、MUSIC（MUlti-SImulation
Coordinator）というインターフェイスを用いることに決定した。

21３）研究開発項目の明確化
NEST2.1 に書き換えたモデルのシミュレーションの結果、庄野らの結果を再現することができた（図 1）。その結果を解析したところ、シミュレーションでは淡蒼球内節（GPi）、淡蒼球外節（GPe）
と視床下核
（STN）
が外部入力の無い状態でも約 20Hz のバーストパターンを示していることが解った。これは生理実験（文献 2）の結果と異なる振る舞いであり、モデルの修正が必要なことが明らかになった（図 2）。
このモデルを修正するためには、1）Hodgkin-Huxley タイプのニューロンモデルのパラメータ
空間の探索、2）パラメータ化された GPi-STN のコネクティビティのパラメータ空間の探索が必要となる。しかしながら、1）に関しては Hodgkin-Huxley タイプのニューロンモデルはパラメー
タ空間が大きいという困難がある。また、2）に関しては Hodgkin-Huxley タイプのニューロンの計算量が大きいためにコネクティビティのパラメータ空間を広く探索するのが難しいという障害がある。そのため、現状のモデルで生理実験の再現が難しい場合には、パラメータ空間が小さい Izhikevich タイプのニューロンモデルを用いてニューロンモデルのパラメータ推定を行い、その上でコネクティビティのパラメータ推定を行うという方法も開発項目として考えている。
４）モデル構築と計算手法開発
モデル構築のツールとしては、大規模なネットワークのシミュレーションの実績があり、「京」への移植作業が進められている NEST をシミュレータとして用い、モデルの開発には PyNEST を使うことに決定した。庄野らのモデルでは 2 つのチャネルの競合で行動が選択されるが、これを複数のチャネルに拡張し、選択される行動が時間と共に切り替わるようなネットワークの構築を進めている。また、多チャンネルのスパイク活動記録データが得られた場合に、ニューロン間のシナプス結合パターンを自動推定するためのアルゴリズムの開発を行った。
（２） ソフトウェアの開発・高度化の状況
１）大脳基底核モデルのシミュレーション用コード開発
まず、NEST1.9 で書かれた庄野らの大脳基底核モデルのニューロンモジュールを、「京」へのポー
ティングが終わり現在チューニングが進んでいる NEST2.1 に書き換えた。また近い将来に予定しているモデルの大規模化に備え、SLI で書かれているコードを、大規模なモデルの開発やデバッ
グがしやすい PyNEST（SLI の python ラッパー）に書き換えた。

22

全 890 ノードの Striatum-GPe-STN-SNr モデルの Thread 並列シミュレーションの結果は、
1theadで 77.21 秒(標準偏差 1.82)、2thread で 43.13 秒（標準偏差 0.89）、4thread で 33.25 秒（標準
偏差 2.32）となっており良好なスケーラビリティが見られる。
今後、大脳基底核モデルを視床、大脳皮質モデル、小脳モデルと結合するための準備として、
INCF のイニシアチブのもとに開発中の XML 言語（NineML）で記述した神経細胞のモデルから、NESTや Neuron で使える神経細胞モデルのコードを直接書き出せる nemo というプログラムを開発した。
また NEST シミュレータ用に生成されたコードの振る舞いを Neuron シミュレータ用のコードの振る舞いと比較するプログラムも開発した。
文献１
Shouno, O., Takeuchi, J., and Tsujino, H. (2009). A spiking neuron model of the basal ganglia
circuitry that can generate behavioral variability. The Basal Ganglia IX, 58, 191–200
文献２
Tachibana, Y., Iwamuro, H., Kita, H., Takada M., and Nambu A. (2011). Subthalamo-pallidal
interactions underlying parkinsonian neuronal oscillations in the primate basal ganglia.
European Journal of Neuroscience, 34, 1470–1484

23