戦略課題 3：予測医療に向けた階層統合シミュレーション［統括：高木 周（東京大学）］

Ⅲ

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

Ⅲ－１

高木 周（東京大学）

心筋梗塞シミュレーションに向けた血栓シミュレータと冠循環モデルの連成及び筋繊維の集合体としての骨格筋のシミュレーション
実施計画
本研究では、「予測医療に向けた階層統合シミュレーション」研究の一環として、「心筋梗塞シミュ

Ⅲ－１－１

レーションに向けた血栓シミュレータと心臓シミュレータの連成及び筋繊維の集合体としての骨格筋のシミュレーション」の研究開発を実施し、研究開発を統括する。
また、プロジェクトを推進する上で、関連する研究者と必要な協議等を行い、予測医療に向けた階
層統合シミュレーション」の研究を統括する。
平成 26 年度は、心筋梗塞シミュレーションに向けた血栓シミュレータと心臓シミュレータの連成として、本課題で開発を進めている血栓症シミュレータに関して、抗血小板薬の薬効を評価する
ためのモデルの精緻化を行い、 P2Y12 や GPIIb/IIIa の阻害薬に関するシミュレーションを実施する。
筋繊維の集合体としての骨格筋シミュレーションの研究課題に関しては、これまで開発を行ってきた運動ニューロンからのシグナルにより筋繊維が収縮し、その集合体として骨格筋全体が収縮するモデルについて、これまでは下腿三頭筋を中心にモデルの構築を行ってきたが、同モデルの全身筋骨格系との連成に向けて、中村仁彦（東大）が開発中の全身モデルへの導入を行い、大規模並列化のための調整を行なう。

1

Ⅲ－１－２

実施内容（成果）

Ⅲ－１－２－１

抗血小板薬の薬効再現に向けたシミュレーション（１）研究開発の実施状況
１）はじめに血小板血栓が形成されていく過程の中で、血管壁に粘着した血小板が活性化反応を伴い血小板同士が凝集していく段階では、血小板の活性化を惹起する物質によってもたらされる刺激が重要な役割を果たすことが知られている。惹起物質が血小板表面に存在するそれぞれに対応した受容体と結合することにより刺激を受け、それに応じて血小板内の活性化反応を生じさせ、カルシウムイオンの増加等の活性化状態に遷移させていく。活性化した血小板は、惹起物質を含んだ顆粒を放出し、放出された惹起物質が自らや近傍にある血小板の活性化を誘発する、ポジティブ・フィードバック機
構によって血小板の活性化の連鎖が促進されていくことになる。 本研究では、このような血小板を活性化する惹起物質と血小板表面にある受容体との結合、それによってもたらされる血小板の活性化、活性化によって生じる顆粒(惹起物質)の放出と活性化の増幅の動作を模擬する、血小板活性化の計算モデルについて検討を進めている。
２）血小板活性化の計算モデル
昨年度までは、血小板の活性化を惹起させる物質として ADP(アデノシン二リン酸)のみを対象に計算モデルの構築を行ってきたが、今年度は、惹起物質としてトロンビン、セロトニン、トロンボキ
サン A2 を新たに追加した。特にトロンビンは、血小板表面にある受容体との結合により血小板の活性化を促進することに加え、血液凝固過程において重要な働きをするフィブリン生成の過程において必要な物質である。血小板表面では、血液凝固因子の Va 因子と Xa 因子が結合しプロトロンビ
ナーゼ(プロトロンビン活性化複合体)が形成される。これに II 因子(プロトロンビン)が結合して活性型のトロンビンが大量に生産され、これが血液凝固過程におけるフィブリノゲンからフィブリンを生成する反応に使用され、フィブリン網の形成を促進させることになる。本年度は、血小板表面でのトロンビンの生成についても追加を行った。
i）惹起物質と受容体
血管壁に粘着した血小板表面において、血小板活性化を惹起する物質(ADP、トロンビン、セロト
ニン、トロンボキサン A2)が、それぞれの物質に対応した血小板表面に存在する受容体と結合する
ことにより刺激を受け、血小板内での活性化反応が進行していく。
計算モデルで対象とした、レセプタとリガンド、およびアンタゴニストを表 1.1 に示す。
表 1.1 計算モデルで対象とするレセプタとリガンド
レセプタ
P2Y12

個数
ADP 受容体

リガンド
ADP

600

(Adenosine
2

アンタゴニ
スト
Clopidogreldiphosphate)
ADP
P2Y1

150

ADP 受容体(Adenosine
diphosphate)

PAR1

トロンビン受容体
トロンボキサン受TP容体

5-HT2A

セロトニン受容体

2,000

Thrombin
TXA2

1,000

(Thromboxane A2)
5-HT

1,000

(Serotonin)vWF

GP IIb/IIIa(von Willebrand factor)

20,000粘着蛋白質

Fibrinogen

GP IIb/IIIa阻害薬ii）活性化モデル
惹起物質の刺激による血小板活性化の反応モデルを図 1.1 に示す。
ADP

5-HT

TBN

TXA2

VWF

FIB

CLO

IGP
(Secretion)

P2Y1

P2Y12

5-HT2A

PAR1

TP

SIG3

SIG4

SIG5

GP IIb/IIIa

DG
SIG1

SIG2

SIG6

血小板

Ca2+

血管壁
リガンド

レセプタ

ADP : アデノシン二リン酸
5-HT : セロトニン
TBN : トロンビン
TXA2 : トロンボキサンA2
VWF : フォンビレブラント因子
FIB : フィブリノゲン

P2Y1: ADP受容体
P2Y12: ADP受容体
5-HT2A : セロトニン受容体
PAR1: トロンビン受容体
TP: トロンボキサン受容体
GP IIb/IIIa : 粘着系蛋白質

アンタゴニスト
CLO : クロピドグレル
IGP : GP IIb/IIIa 阻害薬

図 2.1 血小板の活性化モデル
血小板表面の受容体(レセプタ)と惹起物質(リガンド)とが結合する確率モデルは次とする。第 m 番
目のレセプタがリガンドと結合可能な状態とする場合、レセプタ m の近傍のリガンド濃度を lm としたとき、このレセプタがリガンドと結合する確率分布関数を、
 l 
 l 
rm l m   1  1  m   exp  m  ,
 B 
 B 

0  rm  1
(1.1)

として考える。
血小板活性化モデルの動作は次のようになる。
①血小板表面の受容体(レセプタ)が活性化を惹起物質(リガンド)の濃度に応じて結合する確率モデ
ルを適用する。
3

②リガンドと結合することによって刺激を受け、血小板内に活性化シグナルが生じる。
③この活性化シグナルの累積に応じて、血小板内のカルシウムイオン(Ca2+)が増加し、活性化状態となり、濃染顆粒の放出が生じる。
④濃染顆粒の放出により、活性化の惹起物質である ADP とセロトニンの濃度が上昇する。
⑤活性化シグナルの累積によって増加した Ca2+濃度に応じて、GPIIb/IIIa を活性型にして、フィブ
リノゲン／フォンビレブラント因子との結合が可能な状態に移行する。
⑥活性化した GPIIb/IIIa とフォンビレブラント因子とフィブリノゲンとの結合には、リガンド－レ
セプタの確率モデルを適用する。
⑦酵素反応速度論(Michaelis-Menten 式)に基づいて、血小板表面でのトロンビンの生成濃度を計算する。
iii）濃染顆粒に含まれる惹起物質
血小板活性化の計算モデルに活性化物質の種類を追加したので、濃染顆粒の放出時に印加される活性化物質に、セロトニンを追加した。血小板が活性化することに伴い濃染顆粒が放出される。通常、一つの血小板に含まれる濃染顆粒の個数は 5 個程度である。昨年度までは、惹起物質として ADPのみとしていたが、今年度はセロトニンを追加した。顆粒の放出は 5 回とし、放出量は文献値を参考に、ADP は 653mM、セロトニンは 65mM とした。放出の場所は、血小板表面の 1 計算格子分をランダムに選択し、放出量の ADP 濃度およびセロトニン濃度を印加する。放出の契機は、血小板内の Ca2+貯蔵量が閾値に達した時点で顆粒放出が起こるものとし、放出後は Ca2+貯蔵量を初期状
態にリセットする。
iv）血小板表面でのトロンビン生成
血小板表面では、血液凝固因子である Va 因子と Xa 因子が結合しプロトロンビン活性化複合体(プロトロンビナーゼ)が形成される。この複合体にトロンビンの前駆体である II 因子(プロトロンビン)が結合して活性型のトロンビンが大量に生産される。計算モデルでは、血小板表面にプロトロンビ
ナーゼ複合体を 2,000 個配置することとした。ある一つの複合体について、その複合体が存在している血小板表面上の計算格子での II 因子の濃度に応じて、確率モデルに基づいて II 因子との結合を決定する。II 因子と結合した複合体は活性化状態に遷移する。活性化した複合体では、酵素反応
速度論で用いられている Michaelis-Menten 式に基づいてトロンビンの濃度を計算する。計算で使用する Km(Michaelis 定数)および Vmax(最大反応速度)は文献値に基づいて行うこととした。トロンビンの濃度を TBN(i、j、k)、II 因子の濃度を FII(i、j、k)とすると、Michaelis-Menten 式から、それぞれの反応速度は次で得られる。
V
 FII i , j , k dTBN i , j , k   max
dt
K m  FII i , j , k 

(1.2)

V
 FII i , j , k dFII i , j , k    max
dt
K m  FII i , j , k 

(1.3)

4

(a) t=0.0 msec(b) t=0.623 msec(c) t=0.625 msec(d) t=1.25 msec

図 2.2 複合体の活性化とトロンビン濃度
図 2.2 は、プロトロンビン複合体(紫色)に因子 II が結合して活性化(赤色)し、トロンビン(青色)を生成して濃度が上昇する計算結果例を示している。
４）まとめ血小板が凝集に至る過程で重要な役割を果たす惹起物質の刺激による血小板活性化の計算モデルについて検討を進めている。昨年度までは ADP のみを対象としていたが、惹起物質として新たにトロンビン、セロトニン、トロンボキサン A2 を加えた血小板活性化の計算モデルを構築した。更にトロンビン表面での産生計算を導入し、これにより血小板活性化の動作について確認した。
参考文献[1] 池田康夫，丸山征郎（編）， 血小板生物学， メディカルレビュー社 (2004)．
[2] 池田康夫（監修）， 血栓症， 南江堂 (2004)．
[3] 内山，堀（編）， 抗血小板薬の新しい使い方， 医薬ジャーナル社 (2006)．
[4] Alan D． Michelson， Platelets， Academic Press (2007)．
[5] Noriko Tamura et．al．，”Important Regulatory Role of Activated Platelet-Derived Procoagulant
Activity in the Propagation of Thrombi Formed Under Arterial Blood Flow Conditions，” Circ J 2009;73: 540–548．
5[6] 中川正雄， 真壁利明，確率過程，培風館 (2009)．
[7] Holm Holmsen et．al．，”Secretable storage pools in platelets”，Ann． Rev． Med．， Vol．30，
pp．119-134，1979[8] William H．Kane et．al．，”Blood Coagulation Factors V and VIII: Structural and Functional
Similarities and Their Relationship to Hemorrhagic and Thrombotic Disorders．”，BLOOD VOL 71，
NO3，MARCH，pp．539-555，1988．

6

Ⅲ－１－２－２

マルチスケール骨格筋シミュレータの開発（筋骨格系 FEM モデル）

本プロジェクトでは，コンピュータシミュレーションによる脳神経疾患の運動機能障害の病態予
測と治療支援を目指して，ヒト全身の神経-筋-骨格系の統合シミュレーションシステムの構築を行なっている．本システムは，Fig.1 に示すように，脳神経系シミュレータ，脊髄神経系シミュレー
タおよび筋骨格系有限要素法（FEM）シミュレータを統合したシステムである．
脳で生成された運動の司令は，脊髄内で感覚受容器からのフィードバック情報などと統合され，運動ニューロンを通して筋線維へ伝達される．脊髄から送られた運動指令により筋線維は収縮し，その結果として，関節運動が生じることになる．本統合シミュレータはこれら一連のダイナミクスの再現を目指している．

Fig. 1 Integrated simulation system of neuro-musculo-skeletal system for human whole body.
本研究では，CT や MR などの医用画像から再構成した 3 次元モデルを使用し，生体軟組織の変形と骨格筋の筋活動挙動および関節運動の解析を行うための筋骨格系有限要素法（FEM）シミュ
レータの開発を行っている 1-3)．本年度は筋骨格系 FEM シミュレータの高精度化に向けた動的解析
モデルの導入，計算速度向上に関する検討および単関節モデルから全身筋骨格モデルへの拡張に向けた並列化について検討を行った．
(1)

動的解析モデルの導入

骨格筋の FEM シミュレーションモデルの定式化において，骨格筋の慣性の影響を考慮をするた
めに加速度項を考慮すると，解くべき仮想仕事の原理式は以下のように与えられる．

   u  udV    WdV   S
V

0

V

V

act

:  ΕdV   t s   udS    0 g   udV
S

V

(1)

ここで，V および S は基準配置において連続体が占める領域および境界である．  は後に続く物理 は加速度ベク量の任意の変分である． 0 は初期配置における連続体の密度，u は変位ベクトル， uトルである． S act は筋収縮による第 2Piola-Kirchhoff 応力テンソル，E は Green-Lagrange ひずみテ7

ンソルである．t s および g はそれぞれ表面力ベクトルおよび重力加速度である．W は連続体のひずみエネルギーで次式で与えられる 4)．

W  W (C )  W H  Q(2)

ここで，右辺第 1 項は等容変形（ C  J  C ， C は右 Cauchy-Green 変形テンソル）によるひずみエ
2
3

ネルギーの寄与分，第 2 項は体積変化によるひずみエネルギーの寄与分で，WH=k/2(J-1)2 ある（ J  det F ，F は変形こう配テンソル）．体積弾性率 k は非圧縮性を微圧縮性として導入したことによる体積変化を拘束するためのパラメータである．第 3 項は混合法の導入によって付加される項で，
Q=-2/k( p -p)2，要素レベルで導入される圧力変数 p と変位場から計算される静水圧 p の不つり合いを補正する項である．
平衡方程式は次式の等価節点力ベクトル f および f p を考慮して，

W S act  dV   [ Nu ]T t s dS   [ Nu ]T 0 gdV
S
V
 E


 f   V [ Nu ]T 0 udV  V [ B]T : 

 f p   V [ N p ] Qp dV(3)

T

次式で与えられる 5)．

 u  f    p  f p   0
T

T

(4)

ここで，[ N u ], [ N p ] および [ B] は形状関数およびその微分で，有限要素節点上で定義される変位{u}，
圧力{p}を用いて，それぞれ次式で与えられる．

u  [ N u ]{u}, p  [ N p ]{ p},  E  [ B ]{ u}(5)さらに，微小増分の関係は次式で表される．

M

0

 G 0   u   KUU
0   u
0   p   0 0   p   K PU

KUP   u   f K PP   p   f p (6)

左辺第 1 項および第 2 項が動的解析モデルの導入により，新たに追加された項である．その係数はそれぞれ次式で与えられる．

M   [ N u ]T  0 [ N u ]dV(7) S pas S act G   [ B]T :  :[ B]dV

V
E 
 E

(8)

V

ここで， S pas  W E である．

(6)式には未知量として変位だけでなく，速度や加速度が含まれているため，これらを何らかの関係で結びつけることにより，一変数に減じて解析が行われる 6)．ここでは，位置，速度，加速度の関係として，次式で与えられる Newmark β 法を用いた．



t t



t t



  (1   )  t uu   t u  t   t t u





1
       t u u   t u  t  t u  t 2    t t u
2



8

(9)

ここで，β および γ は定数である．
時刻 t における位置，速度，加速度が既知であるとして，時刻 t+Δt の解を Newton-Raphson 法により求める．k 回目の反復において，次式の関係が満たされるので，

( k ) }{u( k ) }  t 2  {u(10)( k ) }{u ( k ) }  t {u
最終的に，解くべき方程式は次式で与えられる．

 M  t  M
ここで， M0t tG ( k 1)  t 2  t t K ( k 1)u  f 

0   G 0    KUU,G
, K  K0 
 0 0
 PU( k 1)

(k )

(11)


KUP 
 f 
u}    , {f }  , {u である．

K PP  p 
f p 動的解析モデル導入の検証のため，片端固定の棒の縦振動問題において，理論解との比較を行った．自由端の変位の時間履歴について，理論解と FEM による計算結果を Fig.2 に示す．FEM による計算結果は理論解と良い一致を示している．また，Fig.3 に(a)有限要素の分割数および(b)時間増
分量 t が棒の振動数へ及ぼす影響を示す．それぞれ，有限要素の分割数および時間増分量を十分小
さくとることで，理論解と FEM 解析結果の振動数の相対誤差は十分小さくなることが確認できる．

Fig. 2 Time course of displacement at free edge(a) Effect of number of elements(b) Effect of time increment

Fig.3 Relative error of frequency between theoretical and FEM results9(2) 有限要素法解析の計算速度向上に関する検討
有限要素法の計算は時間増分ごとに，(11)式で記述される連立一次方程式を解くことで進められる．有限要素解析時間の大部分はこの連立方程式の係数マトリックスの作成と連立一次方程式の求解に費やされる．そこで，これら計算にかかる時間の短縮について検討した．本研究では，一つの骨格筋（有限要素モデル）の計算に，「京」の 1 ノードを割り当て，ノード内の並列化のみ考慮することにした．
有限要素法では，有限要素ごとに要素の係数マトリックスが独立に計算され，それを足し合わせることで全体の係数マトリックスが作成される．そこで，OPENMP（共有メモリ複数コアマシンの並列化 API）を用いて，要素係数マトリックス作成部分の並列化を行った．
要素数 3880 の有限要素モデルにおける，1 回の全体係数マトリックス作成にかかる時間とその並列化の効果を Table 1 に示す．コア数の増加とともに，係数マトリックス作成にかかる時間はほぼ線形に短縮し，並列化前に比べて約 1/8 の時間となった．

Table 1 Time spent creating global stiffness matrix
Number of cores
Time [sec]

8

4

2

1

0.049

0.096

0.191

0.382

連立一次方程式の求解では，LU 分解法による直接法を用いることとした．LU 分解法による連立一次方程式の求解において，現在，様々な数値計算ライブラリが開発されている．そこで本研究では，いくつかの線形計算ライブラリを用いて，「京」上で有限要素解析の計算時間を比較することにした．ここでは，標準的な数値計算ライブラリである LAPACK 7)，FUJITSU SSLⅡの疎行列ソル
バー， SuperLU 8)および SuperLU_MT 8)について比較を行った．
各要素数 3880 の有限要素モデルにおける，各線形計算ライブラリの 1 回の計算ステップ（全体係数マトリックスの作成，LU 分解および求解）の計算時間の比較を Table 2 に示す．SuperLU_MT（共有メモリ版）の 8 コアを用い場合の計算結果が最も短時間であり，標準的な数値計算ライブラ
リ LAPACK を用いた場合と比べて，約 180 倍高速化した．

Table 2 Computation time of 1 iteration step
SSLⅡ

LAPACK

SuperLU

SuperLU_MT

Number of cores

1

8

4

2

1

8

4

2

Time [sec]

377.4

65.9

113.4

214.7

2.70

2.10

2.55

3.47

(3) 全身筋骨格モデルへの拡張に向けた並列化
ヒト全身の筋骨格系は 400 個以上の骨格筋および 200 個以上の骨と 250 以上の関節からなり，全身筋骨格モデルの解析ではこれらをすべて考慮する必要がある．さらに，骨格筋の変形をミクロス
ケールからマルチスケールで計算するためには，細かな有限要素分割が不可欠であり，計算規模の増大にともない，効率的な並列化処理が必要となる．
骨格筋の変形は脊髄の運動神経からの入力による筋線維の収縮により生じる．脊髄内の運動神経
10とそれが支配する筋肉は筋ごとに独立して存在するため，各骨格筋の変形は，筋同士の接触などが無い限り，独立に考えることができる．そこで本研究では， 1 つの筋肉に対し 1 ノードを割り当てることで，全身の骨格筋を MPI により並列化する．

Fig.4 は 15 個の筋肉について並列化を行った計算例である．各筋肉は筋収縮のタイミングの違い
で 5 つのグループに分けられている．Fig.4(b)はある時刻の筋の横方向の変位を示したもので，筋収縮のタイミングが早いほど筋の膨らみは大きくなっており，同一グループ内の筋の変形は同一となることが確認できた．

Fig.4 Multi muscles simulation by MPI parallelization
上の例では骨格筋の変形のみを考えたが，全身の筋骨格系を考える上では，骨や関節を考慮する必要がある．考慮すべき筋骨格系の計算モデルの概略を Fig.5(a)に示す．骨格筋は骨と付着しており，筋収縮による力が骨へ伝達されると関節運動が生じることになる．筋骨格系の計算では，骨格筋の変形に加えて，関節でつながる骨の運動を考える必要がある．ここでは，骨を剛体と仮定し，関節を拘束条件として与えることで，剛体リンク系の運動方程式を定式化した．筋骨格系の計算で解くべき方程式は，骨格筋の運動と剛体リンク系の運動，さらに，骨と骨格筋の有限要素節点の結合を表す拘束条件を Lagrange 未定乗数法を用いて導入することで導出される 9）．
筋骨格系の計算の並列化では，1 つの骨格筋の計算に 1 つのノードを割り当て，各ノードが 1 つの筋の変数と関連する拘束条件に関する Lagrange 未定定数を管理する．剛体リンク系の自由度は骨格筋に比べて小さいため，剛体リンク系の変数と関節の拘束条件に関する Lagrange 未定定数は 1つのノードで管理する．骨と骨格筋の結合により，剛体リンク系の情報を通信する必要があるが，通信するデータ量は各骨格筋ノードで O((骨と接続する節点数)2) である 9)．

Fig.5(b)は下腿三頭筋（ヒラメ筋・腓腹筋）と前頸骨筋からなる足関節の筋骨格モデルの計算例である．下腿三頭筋への筋収縮指令により，各筋および剛体リンク系の運動が並列計算され，最終的に足関節の底屈運動が再現されている．

11

(a)

(b)
Fig.5 Simulation of musculo-skeletal system

まとめ
中枢神経系および筋骨格系の統合シミュレーションに向けて，筋骨格系計算モデルの高精度化および計算速度の高速化を行った．さらに，全身筋骨格系の計算に向けた効率的な並列化について検討を行った．筋骨格系計算モデルの高精度化では動的解析モデルの導入を行い，棒の縦振動問題において，理論解と比較することで本手法の検証を行った．また，有限要素計算の連立一次方程式の係数マトリックスの作成部分を OPENMP を用いて並列化を行い，さらに，数値計算ライブラリの線形ソルバーを検討することで，従来に比べて 180 倍の高速化を実現した．さらに，全身筋骨格モデルへの拡張に向けた並列化では，1 つの筋肉を 1 つのノードに対応させることで，通信量を大幅に増加させること無く，大規模な筋骨格系の計算が可能となった．今後，筋骨格系シミュレータを脳神経系シミュレータや脊髄神経系シミュレータと統合し，脳神経からの運動指令による筋骨格系の運動を再現する．さらに，安静時振戦や固縮などのパーキンソン病の病態再現を目指して開発進める予定である．
参考文献

1)

Naoto YAMAMURA, Jose Luis ALVES, Toshiaki ODA, Ryuta KINUGASA and Shu TAKAGI, "Effectof tendon stiffness on the generated force at the Achilles tendon - 3D finite element simulation of a humantriceps surae muscle during isometric contraction", Journal of Biomechanical Science and Engineering, Vol.9,
No.3 (2014), p.13-00294.
2) 山村 直人，高木周，“運動単位を考慮した骨格筋の 3 次元 FEM シミュレーション”，日本機械学会 2012 年度年次大会講演論文集，（2012）．

3) 山村 直人，清水和弥，高木周，“パーキンソン病振戦の 3 次元 FEM シミュレーション”，日本機械学会 2013 年度年次大会講演論文集，（2013）．

4) Sussman, T. and Bathe, K.-J., “A FE Formulation for Nonlinear Incompressible Elastic and Inelastic
Analysis”, Computer & Structures, 26 (1987), pp. 357-409.
5) 鷲尾 巧, "マルチスケール•マルチフィジックス心臓シミュレーションに関する研究−サルコメア
力学から心筋細胞構造を経て心拍動に至る解析手法の開発と応用−", Ph.D. Thesis，The University ofTokyo，2013.
6) 久田 俊明，野口裕久，”非線形有限要素法の基礎と応用”，丸善．
7) http://www.netlib.org/lapack/8) http://crd-legacy.lbl.gov/~xiaoye/SuperLU/
129) 平澤 謙章，”スパイク脊髄神経モデルと全身筋有限要素モデルの統合とその超並列計算”，Master
Thesis, The University of Tokyo，2014.

13

Ⅲ－２

久田俊明（東京大学）

心疾患の合理的治療のためのマルチスケール・マルチフィジックス心臓シミュレーショ
ン
Ⅲ－２－１ 実施計画
本研究では、「予測医療に向けた階層統合シミュレーション」研究の一環として、「心疾患の治療法のためのマルチスケール・マルチフィジックス心臓シミュレーション」の研究開発を実施し、研究開発を統括する。
本年度は昨年度に引き続きマルチスケール・マルチフィジックス心臓シミュレータ UT-Heart
のサルコメア動力学モデルの高度化とその応用を図ると共に、冠循環マルチスケールシミュレー
ションを通じて新たな医学的・生理学的知見を得ることを目標とする。具体的な例としては、冠動脈カテーテル治療（PCI）において、適応の判定のための指標として広く用いられている一方でその生理学的妥当性に疑問が呈されている FFR（Fractional Flow Reserve）の真の意味を明らかにし、より合理的な指標について検討する。更に、高木グループで開発が進められている血栓症シミュレータにおいて血小板粘着から血栓形成のプロセスを再現するため、UT-Heart から冠血管径変化、流量、圧などのデータを提供する。これにより最終年度の両シミュレーションプロ
グラムの 連成へ向けての準備とする。
また、平成２５年度に行われた中間評価の指摘に従い「京」以外の計算資源を積極的に活用するため、研究室が所有するスパコンを４か月間の間、本課題実施のため専有して利用し、その電気使用料を本委託費より支出する。
Ⅲ－２－２

実施内容（成果）（１） サルコメア動力学モデルの高度化とその応用
１）理化学研究所生命システム研究センターの Lorenzo Marcucci 研究員との共同研究により、
化学エネルギーと力学的エネルギー間の変換を定めるミオシン分子パワーストロークモデルの改良を試みた。現状では図 1 上段(黒)に示すように、ミオシンヘッドの回転量を反応座標とする
狭い井戸型の自由エネルギーPMF(Potential of mean force)モデルを想定し、限定された狭い範囲にミオシンヘッドの状態が局在する状態間の遷移を統計力学的考察からモデル化している。そこで改良モデルにおいては、図１上端(赤)のように幅をもったポテンシャルを想定し、ミオシン
ヘッドの状態の拡散効果を考慮して広い領域内での分布を許容する状態遷移モデルを実装し、ミクロレベルでの分子の振る舞い、筋線維レベルでの長さ変化と収縮力の関係性、および臓器レベ
ルでの拍動性能に対して従来モデルと比較してどのような違いが現れるかを検討した。
今回の改良により、各レベルにおいて以下のような効果が現れた。
・分子レベル：異なる状態間で近接する状態同志を許容することになった結果、揺らぎによる遷
移効果がより覿面に現れ状態間遷移を許容しやすくなった。
・筋線維モデル：急激な長さ変化に対する張力応答などの実験結果がより観測結果に近い形で再現できるようになった。また、少ないミオシンヘッドの結合率でより大きな収縮力を発生できるようになった。
・臓器レベル：より少ない ATP 消費量で正常拍動が再現できるようになり、拍動性能の効率が上がった。
14

一例を挙げると、臓器レベルにおいては図２左に示すように拡散型のモデルにより大幅に左心室における ATP 消費が削減されていることがわかる。さらに、図２右のオレンジ、緑、赤の線は強結合状態にあるミオシン分子濃度の時間変化を示すものであるが、拡散型モデルでのその割合が非常に少なくなっていることがわかる。

図 1 ミオシンヘッドパワーストロークポテンシャルモデル

図 2 エネルギー消費(左)と状態濃度の変化(右)の比較
15

２）先天性心疾患外科手術については、従来の UT-Heart マクロ解析機能の充実と安定化作業を行うとともに３症例の術前術後シミュレーションを実施し、結果の妥当性を確認した。具体的には、以下の解析がより安定に実行できるように改良を行った。
・臓器外壁および弁尖間の接触解析
・弁尖の折れ曲がりによる不安定性の回避
・弁尖と血流の相互作用の取扱い
・流体部 ALE メッシュ変形計算の安定化
以上の改良の結果、上記の３症例に対して安定に拍動解析が実行できるようになった。さらに、流体部の流入出部に接続する体循環、肺循環モデルの回路パラメータを患者の血圧データから自動的に決定するシステムも実装した。図 3 にそのうちの１症例のシミュレーション結果を示す。図 3A はシミュレーションによる動脈血と静脈血の混合の様子を表している。これは算出された血流の流速データのもとで移流拡散方程式を解くことにより得たものである。心室中隔欠損を通して酸素を含んだ動脈血が右室起始の大動脈に送り込まれる様子を表している。B は
この症例に対して適用された手術のイメージ図であり中隔欠損より大動脈にかけて人工血管により再建流路を構築し、左心室と右心室を分離することを表している。C はシミュレーショ
ンによる術前、術後の血流の様子、D は ATP 消費の比較を示したものである。この拍動解析は、先に述べた分子モデルに基づくものであり、血行状態のみならず、心筋にかかる負荷を評価す
る上でも分子モデルは重要である。この例では、術後の右心室壁のエネルギー消費面での負荷が軽減されていることがわかる。

図 3 先天性心疾患外科手術シミュレーションの１症例

16心臓サポートネットに関しては、右室拘束軽減ネットの有効性を確認した。本研究には研究室のスパコンを使用した。心筋サポートネットは、拡張型心筋症に適用され、心室をネットで覆うことにより心室の拡張を抑えることを目的とするものである。本研究においてはネット装着
が拍動性能に及ぼす影響やネット着圧の評価などをシミュレータにより行った。
３）リモデリング機構については、マクロレベルでのモデル化を行い、その有効性を確認した。
心筋細胞の特徴である枝分れ構造(図 4A)に着目し、それをマクロレベルでモデル化した心室モ
デルで拍動解析を行い(図 4B)、各拍動で得られた局所的な情報を利用して線維構造の主方向を改変していくように線維構造リモデリング機構を実装し、ほとんど一様に水平な線維構造から出発した場においても、実測と一致する内壁から外壁にかけた捻じれ構造が現れることを確認した(図 4C)。さらに従来からその発生理由が謎であった線維に沿ったバンド構造が再現できることも確認した(図 4D,E)。

図 4 線維構造リモデリング解析（２）冠循環マルチスケールシミュレーションに関する研究
UT-Heart には冠循環がマルチスケールモデル、すなわちマクロとミクロの複合モデルとして実装されている。図 5 はマクロモデル、すなわち心筋壁に埋め込まれた動脈系（赤色）と静脈系（青色）の一部をズームアップして示したものである。これらのマクロモデルの血管系はある圧力差に対して血液が心筋の各所に均等に流れるように生成されている。

17

図 5 マルチスケール冠循環マクロモデル（動脈系を赤色、静脈系を青色で示す）

図 6 実測に基づき生成された冠循環ミクロモデル（黄緑色のネットワークは毛細血管）18

一方ミクロモデルは図 6 に示すようにマクロ動脈の抹消とマクロ静脈の末梢を繋ぐ形で接続されている。従って大動脈基部から冠動脈マクロモデルに流入した血液はミクロモデルの微小循
環系、毛細血管を経て最終的に静脈系から右心房へ流出する構造が再現されている。
1）本年度は、先ず UT-Heart のマルチスケール冠循環モデルにおける、マクロ部とミクロ部のパラメータ調整により、各所の流量や圧力分布がより適切なものになるよう改良した。図 7 は大動脈から流入する血液の流量（実線）と右心房へ流出する流量（破線）を示したものである。
4 つのケースはミクロモデル血管の条件を多少変えたものであるが、70～100cc/min/心筋 100g

図 7 冠動脈入口および冠静脈出口での流量の変化
程度の範囲にあり、また流量の波形も動脈は拡張期優位、静脈は収縮期優位と実測を再現するものとなっている。さらに例えば小静脈（直径 840μm）の流速を図 8 に示すが、このシミュレーション
結果も定性的・定量的に実測と良好な一致を示している。

図 8 小静脈における流速19２）冠動脈カテーテル治療（PCI）において、適応の判定のための指標として広く用いられている一方でその生理学的妥当性が議論されている FFR（Fractional Flow Reserve）の検討を行った。図 9 は UT-Heart において左回旋枝（LCX）の基部を狭窄させた場合の流量変化を表す。
横軸は狭窄率を示すが、狭窄率 0 の健常な状態から狭窄を増して行っても流量に殆ど変化はな
く、狭窄率が実に 0.85（血管径が 0.15 まで減少）付近に至った段階で急激に流量が減少し始める傾向にあることが分かる。

図9

UT-Heart の左回旋枝基部に狭窄を与えた場合の流量の変化

FFR は次式で定義される。
(Pd-Pv)/(Pa-Pv)
ここで Pa:狭窄開始点の時間平均血圧、Pd:狭窄終端部の時間平均血圧、Pv:右心房の時間平均血圧である。これに対し、従来から用いられてきたもう一つの指標：CFR(Coronary Flow
Reserve)は安静時における流量と冠血管拡張剤を投与した後の最大冠血流量の比として定義され、冠動脈の機能的狭窄度や微小循環障害の指標として用いられてきた。このように一般に冠血流は末梢の血管の障害の程度に左右され、しかも実際は拍動流であるので、上式で定義される狭窄部の圧力損失比だけから PCI による血流の改善を予測できるかどうかについては自明ではない。解析結果の掲載は控えるが、以上の観点から検討を進めた。
３）毛細血管床における圧や流量を出力し高木グループにデータを提供する準備を行った。

20

Ⅲ－３

後藤

信哉（東海大学）

抗血小板薬薬効予測を可能とする階層統合シミュレータの開発に向けたモデリングⅢ－３－１ 実施計画
本研究では、「予測医療に向けた階層統合シミュレーション」研究の一環として、「予測医療の実現に向けた抗血小板薬薬効予測を可能とする階層統合シミュレータの開発に向けたモデリング」の研究開発を実施する。
また、「予測医療の実現に向けた抗血小板薬薬効予測を可能とする階層統合シミュレーターの開発と大規模計算」の研究を行う上で、関連する研究者と必要な協議等を行うとともに、本格実施に必要な研究体制の整備を行う。
平成２６年度は以下の研究開発を実施する、（１）東京大学・高木グループの「心筋梗塞・脳梗塞のマルチスケールシミュレーション」モ
デルに対して、フローチェンバーを用いた血小板細胞接着動態の実証実験に裏打ちされた、かつ血小板細胞の GPIbαと von Willebrand 因子の分子結合から血小板細胞接着を連成する血小板細胞接着モデルを作成し精緻化する。
（２）実証実験に裏打ちされた血小板細胞代謝、活性化モデルの精緻化を継続する。
（３）実証実験に裏打ちされた血小板細胞活性化と活性化血小板上の凝固系活性化を連成するモデルを作成し精緻化する。

21

Ⅲ－３－２

実施内容（成果）（１）東京大学・高木グループの「心筋梗塞・脳梗塞のマルチスケールシミュレーション」モデ
ルに対して、フローチェンバーを用いた血小板細胞接着動態の実証実験に裏打ちされた、かつ血小板細胞の GPIbαと vonWillebrand 因子の分子結合から血小板細胞接着を連成する血小板細胞接着モデルを作成した。具体的には分子動力学計算の手法を用いて血小板細胞膜糖蛋白 GPIbαの N
末端の von Willebrnad 因子 (VWF)
結合部分と、その結合相手である
von Willebrand 因子の A1 ドメインの結合構造を図 1 に示すように予測した。構造予測において、既に他グループから発表された結晶構造から出発し、水分子を配置して水溶構造を予測した。VWF と GPIαの結合により惹起される力を高木グループのシミュレーターに提供する必要がある。両分子の重心を数Å毎に移動させ、両分子により形成される力場から potential of mean force を図 2 のように算出した。図 2 の上段には、重心間距離を移動させたときの分子構造の変化を示す。各分子構
造について揺らぎの計算を行っている（動画のため本報告書には掲載せず）。図 2 下段に各重心間距離における potentialof mean force の計算結果を示す。計算値である赤+を結んだ曲線を微分すると力を算出できる。両分子の結合により惹起される力が数十 pNであることを示し、高木グ
ループのシミュレーターの精緻化に寄与した。
本法の応用により次のステップとして活性化血小板における活性化構造を呈し
た GPIIb/IIIa と VWF、フィブ
リノーゲンの結合力を見積もるための基盤モデルの作成に成功した。

22

（２）実証実験に裏打ちされた血小板細胞代謝、活性化モデルの精緻化を継続した。しかし、図 3に示すように公表された反応値を用いて血小板細胞内の代謝活性化
モデル計算を行っても、細胞内の化学反応の反応定数が、細胞外と同一である保証がない。
研究としてモデル化は進んだが、現時点での問題点も明らかになった。
（3）実証実験に裏打ちされた血小板細胞活性化と活性化血小板上の凝固系活性化を連成するモデルの作成は定性
的には、現時点にて入れ込み得る因子のほとんどの入れ込みに成功した。モデル（図 4）を論文として公表した（Goto S,

et al. Circulation, 2015）。
従来、作成した定量モデルを、
図 4 に示した精緻モデルを用いてさらに精緻にすることを目指す。

23

Ⅲ－４

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

モデルにおける筋肉の体積効果の大規模並列化及び人の行動モデルに関する大規模計算」の研究開発を実施する。
また、「全身筋骨格モデルにおける筋肉の体積効果及び人の行動モデルに関する大規模並列計算」の研究を行う上で、関連する研究者と必要な協議等を行うとともに、本格実施に必要な研究体制の整備を行う。
平成２６年度は以下の研究開発を実施する、（１）脊髄の筋支配モデルを上肢から下肢へと展開し、全身モデルの研究開発を進める。
（２）沖縄科技大の銅谷らの研究グループで開発されたパーキンソン病脳モデルのソフトウェア(NEST)を、ソフトウェア MUSIC によって中村グループで開発した体積効果を考慮した四肢
筋骨格モデル(K-Body)と接続し統合計算の研究開発を行う。
（３）大阪大学の野村グループの患者データベースを利用することで、パーキンソン病による姿
勢制御破綻のモデルシミュレーションの研究開発を開始する。
Ⅲ－４－２ 実施内容（成果）（１） 脊髄神経系モデルの全身への展開
a. 全身α運動ニューロンの個数の推定
前年度には、全身ワイヤ筋骨格モデル[R3]から全身筋有限要素モデルを作成し、体積効果を考慮した筋骨格モデルを作成した。また脊髄の筋支配モデルの一例として、上腕二頭筋および上腕三頭筋のα運動ニューロンおよび感覚神経としての筋紡錘モデルを導入し、等尺性収縮の運動シミュレーションをおこなった。本年度は、脊髄の筋支配モデルに関し、
前年度までの上肢の一部モデルを全身へと展開することにより、全身での神経運動シミュ
レーションができるモデルを開発した。
筋は脊髄内に存在する複数のα運動ニューロンからスパイク信号を受けて収縮を行う。
過去の文献では限られた筋を支配するα運動ニューロンの個数は解明されているが[R1][R7][R8][R9]、全身の筋に対してそれらを支配するα運動ニューロンの個数を個別に設定できるほどの知見がない。そこで生理学文献を調査して得られる既知の運動ニューロンの個数から、全身の筋それぞれを支配するα運動ニューロンの個数を推定した。推定には筋の質量・筋紡錘の個数の情報を用いる。一つのα運動ニューロンが支配する筋繊維の数(神経支配比) は筋の制御の細かさにより異なる。同じ関節を駆動する筋は神経支配比が同程度となるという仮定のもと、筋を支配するα運動ニューロンの個数、筋の質量、筋紡錘の個数から相関を調査する。その結果、図１に示すようなグループ化を行うことで全身のそれぞれの筋を支配するα運動ニューロンの個数を推定する。推定の結果、α運動
ニューロンは全身で 255,546 個となる。運動ニューロンの数は脳幹及び脊髄に数十万個存在するといわれている[R10] ことから、妥当な数字であると考えられる。

24

図１．α運動ニューロン、筋紡錘および筋の質量に基づく筋のグループ分け

b. 脊髄内α運動ニューロンの配置
α運動ニューロンは脊髄内において、各筋に信号を伝達する役割のニューロンが集団（プール）を構成する。脊髄内での位置情報を基にニューロン間の結合を統計的に設定するため、解剖学的知見[R2] に基づきα運動ニューロンを脊髄内に配置する（図２）。ヒトの下半身の筋に関する運動ニューロンの脊髄内での分布に関して Sharrard[R2]を参考に、
また上肢および体幹の筋を支配するα運動ニューロンの分布は、Bossey[R11]らによる伸筋・屈筋配置に関する記述を参考にした。ある筋を支配する運動ニューロンに対しては、
脊髄長軸方向には運動ニューロンプールが存在する脊髄層の範囲において一様乱数で位置を決定する。脊髄断面方向は、ニューロンプールの中心位置から正規分布に従って周辺に点を分散させた。同じ色の点が同じ筋を支配するα運動ニューロンプールを表している。

図２．α運動ニューロンの脊髄内への配置

25

（２）全身神経筋骨格統合シミュレーション
a. スパイク入力による全身運動の生成
前年度では右上腕の神経信号を入力とする運動シミュレーションをおこなった。本年度
は全身へ展開した脊髄神経系のモデルの成果を利用し、α運動ニューロンのスパイク信号による全身有限要素筋骨格モデルの運動シミュレーションを実現する統合モデルを開発した。運動シミュレータとしては、前年度までの成果物であるニューロンシミュレータ
NEST[R4]、骨格筋有限要素解析シミュレータ V-Biomech[R5]、全身筋骨格シミュレータ
sDIMS [R3]、プロセス間の通信を可能にする MUSIC[R6]から構成される統合インターフェ
イスを利用した。
1Hz、10Hz、100Hz のスパイク信号を、上腕二頭筋、大腿二頭筋及び腓腹筋へ入力する。
また、筋の変形から筋感覚受容器の発火頻度を計算する。シミュレーションの実行にスーパーコンピュータ京を用い、筋 374 本各々1 プロセス、骨格 1 プロセス、神経系 166 プ
ロセスの計 541 プロセスで並列計算を行う。図３に 100Hz のスパイク信号を入力した際のスナップショットを示す。100Hz のスパイク信号を入力するとき 1 秒間の運動シミュ
レーションの計算時間は 1405 秒であった。

図３．α運動ニューロンからの 100Hz スパイク信号を上腕二頭筋および大腿二頭筋に入力した際の運動シミュレーション

a. 体性反射シミュレーション
脊髄内α運動ニューロンの位置情報をもとに、同名筋の感覚ニューロンとα運動ニュー
ロンを距離に応じて正規分布で重みづけし、結合することで、体性反射シミュレーションを行う（図４、図５）。シミュレーションでは 0.2 秒に下腿に外力を加え、筋が伸張することにより、筋紡錘発火頻度が増加し、感覚ニューロンがα運動ニューロンに信号伝達することで、α運動ニューロンの発火頻度が増加することが確認した。図４に膝蓋腱反射シ
ミュレーションの 0.2 秒毎のスナップショットを示す．外力が 0.2 秒に加わることにより，一旦屈曲し，その後膝関節伸展筋の筋活動度の増加により，膝が伸展する．図５に膝蓋腱反射シミュレーション結果を示す。図左上は大腿直筋の筋繊維方向歪速度、図右上は筋紡錘発火頻度、図左下にα運動ニューロン発火頻度、図右下は大腿直筋の運動単位筋活
動度をそれぞれ表す。感覚ニューロンと 運動ニューロンは脊髄で結合しているため，外力が加わり，筋紡錘発火頻度が増加する 0.2 秒から筋紡錘の発火が運動ニューロンに信号伝達されるまでの 0.035 秒の潜時の後， 運動ニューロンの発火が増加していることがわかる。

26

図４．膝蓋腱反射シミュレーション。20N の外力を脛骨に対し 0.2[s]のタイミングで印加する。

図５．膝蓋腱反射シミュレーション結果。
（左上）大腿直筋の筋繊維方向歪速度、（右上）筋紡錘発
火頻度、（左下）α運動ニューロン発火頻度、（右下）大腿直筋の運動単位筋活動度
27

（３）パーキンソン病による姿勢制御破綻のモデルシミュレーション
大阪大学の野村グループの患者データベースを利用することで、パーキンソン病による姿勢制御破綻のモデルシミュレーションの研究開発を開始した。具体的には、足底の皮膚の触覚に関わる感覚神経の空間分布モデルを開発した。また前提感覚に関わる三半規管の神経細胞
約 36000 個の空間分布モデルを開発した。図６に NEST で実装された足底の接触感覚神経細胞の空間分布モデルを示す。

15 SAI Neurons (Merkel corpuscle): Blue
16 SAII (Ruffini corpuscle): Green59 FAI (Meissner corpuscle): Red
14 FAII (Pacinian corpuscle): Yellow
図６

NEST で実装された足底の接触感覚神経細胞の空間分布モデル

参考文献[R1] R.M. Enoka. Morphological features and activation patterns of motor units.
Morphological features and activation patterns of motor units. 12, 6, 538-599, 1995.
[R2] W. J. W. Sharrard. The distribution of the permanent paralysis in the lower limb inpoliomyelitis. Journal of Bone & Joint Surgery, 37, 4, 540-558, 1955.
[R3] Y. Nakamura, K. Yamane, and I. Suzuki. Somatosensory computation for man-machineinterface from motion-capture data and musculoskeletal human model. IEEE Transactions on
Robotics, 21(1): pp. 58-66, 2005.
[R4] M-O. Gewaltig, M. Diesmann. NEST (Neural Simulation Tool), Scholarpedia, 2(4): 1430,
2007.
[R5] 山村直人, ALVES Luis, 小田俊明, TEODOSIU Cristian, ヒト骨格筋力学解析のための有限要
素シミュレータの開発. JSME annual meeting: 29-30, 2009.
[R6] 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.
[R7] F. Buchthal, F. Erminio, and P. Rosenfalck. Motor unit territory in different human28

muscles. Acta Physiologica Scandinavica, Vol. 45, No. 1, pp. 72-87, 1959.
[R8] A. J. McComas. Invited review: motor unit estimation: methods, results, and present
status. Muscle & Nerve, Vol. 14, No. 7, pp. 585-597, 1991.
[R9] B. Feinstein, B. Lindegard, E. Nyman and G. Wohlfart. Morphologic studies of motor
units in normal human muscles. Cells Tissues Organs, Vol. 23, No. 2, pp. 127-142, 1955.
[R10] C. T. Leonard. The Neuroscience of Human Movement. Mosby, 1997.
[R11] J. Bossy. Atlas of Neuroanatomy and Special Sense Organs. Saunders, 1970.

29

Ⅲ－５

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

Ⅲ－５－１ 実施計画
「予測医療に向けた階層統合シミュレーション」の研究では、循環器系、筋骨格系、脳神経系の数値モデルを構築し、それらをつなぐ大規模シミュレーションにより、心筋梗塞、パーキンソン病などの発症プロセスを解明し、治療手法の探索評価を可能にすることをめざしている。その一環として、脳神経系の感覚入力から運動出力までをつなぎ、脳全体の回路から細胞、分子レベルのダイ
ナミクスを含む多階層の脳神経系モデルの構築が不可欠である。
平成２６年度は、NEST により実装した大脳基底核モデルと視床–大脳皮質モデルを MUSIC により統合し、健常時の行動選択機能とパーキンソン病モデル動物の実験データを再現するため、シナ
プス、ニューロン、回路の各レベルでの検討を行う。さらにこれを東大高木チーム、中村チームが作成した脊髄神経回路モデルおよび筋骨格系モデルと統合し、単関節運動における振戦と固縮症状の発生機構の探索を行う。また立位姿勢制御モデルの実現に向け、小脳の神経回路モデルの NEST による構築を進める。
また、「行動制御の脳神経系の多階層モデルの構築」の研究を行う上で、関連する研究者と必要な協議等を行う。
Ⅲ－５－２ 実施内容（成果）(1) 大脳基底核モデルの構築（図１）
前年度に提案した視床下核－淡蒼球の局所神経回路網モデルを拡張し、健常状態およびパーキン
ソン病状態で観測される神経活動の特徴を、平均発火頻度、8-15 Hz での単一神経細胞の振動活動、およびバースト発火に関して再現するための、神経細胞およびシナプスのパラメタ探索を行った。
その結果、パーキンソン病患者の大脳基底核での異常な振動活動の発生には、視床下核神経細胞の発火特性および淡蒼球外節神経細胞の自発活動の低下が重要であることが明らかになった。
(2) 視床―大脳皮質モデルの開発（図２）
脊髄―筋骨格系と大脳基底核を結ぶ視床―大脳皮質モデル開発に取り組み、前年度は視床神経回路モデルではパーキンソン病で見られるシータ、アルファ周波数帯の振動的発火現象を再現し、運動皮質モデルでは側抑制により motor primitive に対応する特定細胞集団の活動選択を実現した。
今年度は、神経細胞密度、層の厚みなどについてより現実に即したパラメータの導入を行うとともに(Igarashi et al. 2014)、パーキンソン病の振戦の再現に必要となる約 5Hz の神経振動の逆相同期現象の発生について取り組んだ。これまでに、パーキンソン振戦での腕や指の約 5Hz の振動は、運動皮質で発生する 10Hz 程度のアルファ周波数帯の脳波と 1 対 2 同期をし、伸筋と屈筋の筋電信号は逆相同期を示し、また、視床では振戦に関連する発火活動がみられることが報告されている(Timmerman et al., 2003)。これらの報告から、パーキンソン病でみられる大脳基底核の異常な抑制性振動やドーパミン減少などの影響で、視床では low-threshold Ca2+ spike の 5Hz 前後のバースト発火による神経振動が発生し、視床や運動皮質の側抑制によって振動対が逆相同期して、さらにその信号が運動皮質から脊髄、筋肉へと伝わることで、振戦が発生するという可能性が考えられ
る。

30

図1

STN から GPe への結合強度（横軸）、GPe から STN への結合強度（縦軸）、GPe/GPi の興奮性（横３列）を変化させた時の、STN および GPe の発火頻度（左パネル）、ベータ帯域振動、バー
スト強度（右パネル）の変化。サルの神経活動データ(Tachibana et al., 2011)との比較により、健常状態およびパーキンソン病状態に対応するパラメタ領域を明らかにした。
この可能性を検証するため、我々の構築した視床―運動皮質モデルを用いてシミュレーションを行った。シミュレーションでは、視床網様体と介在神経の側抑制によって、同時に活動する視床皮質細胞の空間的なクラスタが形成され、異なるクラスタ間で振動の逆相同期が起きた。その視床内の活動の時空間パターンは運動皮質へ伝わり、運動皮質の出力細胞である 5B 層の錐体路細胞群の活動は空間的にクラスタ化し、異なるクラスタ間で、振動的活動は逆相同期を示した。ここで、視床と皮質それぞれでの個々の神経細胞を見ると、振戦に関連した 5Hz 程度のバースト様の発火活動がみられるが、集団としての活動では逆相同期により周波数が倍になるため、10Hz 程度のα振動が現われ、各神経細胞は、アルファ振動に位相固定する様子が見られた。これらの結果は、前述の生理実験で見られる、逆相同期、振戦と脳波の 1 対 2 同期の現象とうまく合致する現象である。このシミュレーションで生成された運動皮質の出力に関して、脊髄、筋骨格系モデルにその出力を与え、実際に振戦が起きるかどうかの調査を、東京大学高木研と中村研で行っている。
(3) NEST と MUSIC によるモデル統合（図３）
京コンピュータ上で神経シミュレータ NEST と神経シミュレータ間の相互接続をおこなうライブラリ
MUSIC を用いて、庄野らの大脳基底核モデルと五十嵐らの皮質視床モデルの結合を行った。前年度のモ
デルは、全体で 120,000 神経細胞からなる小規模のもので、大脳基底核のモデルは神経結合の空間構造が導入されていなかった。
今年度、我々はより現実的でスケーラブルなモデルを構築することに重点をおき、大脳基底核モデルに 2 次元の空間構造を導入した。
大脳基底核モデルはラットの大脳基底核全体のスケールとし、線条体（投射ニューロン約 265 万個、介在ニューロン 14 万個）、淡蒼球外節（4 万 6 千ニューロン）、視床下核（1 万 4 千ニューロン）、淡蒼球内節（3 千ニューロン）、により構成される。神経細胞モデルとしては、Izhikevich モデルもしくは31

図 2 （上）視床の振動的発火活動と視床皮質細胞の逆相同期発火。下：視床ー皮質モデルにおけるアルファ周波数帯の振動活動
intergate-and-fire モデルを用いた。視床―大脳皮質モデルは視床の４種類のニューロン、大脳皮質の９種類のニューロンから構成され、運動野の 3 平方ミリの領域のスケールとした。
両者を MUSIC(Multi Simulation Coordinator)により双方向に結合し、「京」の 8192 コアを用いてシ
ミュレーションを行った。６秒間の神経回路動作のシミュレーションは 631 秒で終了した。
統合シミュレーションにおいてパーキンソン病様の振戦状態を示し、約 14Hz の振動が大脳基底核で発生して、視床皮質領域、錐体路細胞に投射する現象がみられている。
参考文献
Igarashi J, Moren J, Yoshimoto J, and Doya K: Selective activation of columnar neural
population by lateral inhibition in a realistic model of primary motor cortex, Societf for neuroscience
2014, 2014, C66

Tachibana Y, Iwamuro H, Kita H, Takada M, and Nambu A: Subthalamo-pallidal interactions
underlying parkinsonian neuronal oscillations in the primate basal ganglia. European Journalof Neuroscience, 2011, 34, 1470–1484
Timmermann L, Gross J, Dirks M, Volkmann J, Freund HJ, and Schnitzler A: The cerebral
oscillatory network of parkinsonian resting tremor. Brain. 2003 126, 199-212.
Shouno O, Takeuchi J, Tsujino H: A spiking neuron model of the basal ganglia circuitry thatcan generate behavioral variability. The Basal Ganglia IX 2009, 58, 191–200.
Shouno O, Doya K: Local circuit model of the subthalamo-pallidal network for the generation of
parkinsonian oscillations. BMC Neuroscience 2014, 15(Suppl 1):P16832

図 3：京コンピュータで大規模実行試験を行った大脳基底核―視床―皮質回路モデルの全体図

33

Ⅲ－６

野村泰伸（大阪大学・大学院基礎工学研究科）
コンプライアントな生体筋・腱系に駆動されるヒト骨格系の運動制御に対して
大脳基底核が果たす役割の数理モデル化

実施計画
本研究では、「予測医療に向けた階層統合シミュレーション」研究の一環として、平成２３年度

Ⅲ－６－１

から開発が進められている「全身筋骨格モデルにおける筋肉の体積効果の大規模並列化及び人の行動モデルに関する大規模計算」に導入可能な脳神経系による運動制御機構のマクロ機能レベルの数
理モデルを構築する。特に、ヒト静止立位姿勢と二足歩行運動の神経制御メカニズムをモデル化の対象とする。また、運動障害を伴う中枢神経疾患のひとつであるパーキンソン病の病態を、本研究で構築する制御メカニズムの崩壊現象として定量的に説明・再現する。
また、「コンプライアントな生体筋・腱系に駆動されるヒト骨格系の運動制御に対して大脳基底核が果たす役割の数理モデル化」の研究を行う上で、関連する研究者と必要な協議等を行う。
平成２６年度は、（１）東京大学・高木グループの「筋繊維レベルからの骨格筋モデル」、沖縄科
学技術大学院大学・銅谷グループの「脳神経系の多階層モデル」、東京大学・中村グループの「脊髄と四肢筋骨格モデル」を統合したインターフェイスの開発と利用を引き続き実施する。
（２）運動に関わる関節の受動的粘弾性が低くかつ能動的神経フィードバック制御のゲインも小さい（関節がコンプライントである）という生理学的に妥当な条件の下で全身機械力学モデルの直立姿勢と歩行を柔軟かつ安定化に実現する制御機構の数理モデル化を平成 25 年度に引き続き実施し、得られた知見を（１）における大規模詳細モデリングに対してフィードバックする。
Ⅲ－６－２ 実施内容（成果）（１）研究開発の実施状況
平成 26 年度は、前年度に引き続き、①ヒト静止立位姿勢制御に関して、関節の受動的粘弾性が低くかつ能動的神経フィードバック制御のゲインも小さい（関節がコンプライアントである）という生理学的に妥当な条件の下で、静止立位姿勢を柔軟に安定化する間欠制御仮説の検証を行った。これにより、東大高木・中村グループ、OIST 銅谷グループが構築する大規模並列モデルが実現すべき神経系による身体運動制御のマクロレベルの特性として、制御の間欠性が重要な役割を果たしていることを明らかにした。特に①-1 として、ヒト静止立位中の姿勢動揺データを統計学的に解析し、間欠制御を特徴付けるベキ関数的長期相関を伴う低周波の姿勢変動が、個々の被験者に依らず普遍的に存在することを明らかにした。また、大脳基底核が重要な役割を果たすと考えられている強化学習によって、間欠制御を実行する制御器が獲得される可能性について検討した。特に、銅谷らが過去に行った強化学習を用いた人型ロボットの立ち上がり動作制御のシミュ
レーションを実装・解析し、類似のアルゴリズムによって静止立位の間欠制御器が獲得できる可能性を示唆することができた。②ヒト二足歩行運動に関しても、関節のスティフネスを小さく抑えた状況でも、静止立位と類似した（より動的な意味において）間欠制御によって、動的に安定な運動が実現できることを、ヒト二足歩行の動的数理モデルの数値シミュレーションに基づき明らかにした。以下では国際論文誌への掲載に至った①-1 と②の成果を具体的に述べる。

１） 成果①（発表論文１）
ヒト静止立位姿勢制御に関して、関節の受動的粘弾性が低くかつ能動的神経フィードバック
制御のゲインも小さい（関節がコンプライアントである）という生理学的に妥当な条件の下で、静止立位姿勢を柔軟に安定化する間欠制御仮説の検証を行った。これにより、東大高木・中村グループ、OIST 銅谷グループが構築する大規模並列モデルが実現すべき神経系による身体運動制御のマクロレベルの特性として、制御の間欠性が重要な役割を果たしていることを明らかにした。特に、ヒト静止立位中の姿勢動揺データを統計学的に解析し、間欠制御を特34徴付けるベキ関数的長期相関を伴う低周波（0.05-0.5 Hz）の姿勢変動が、個々の被験者に依らず普遍的に存在することを明らかにした。さらに、ヒト静止立位中の姿勢動揺に含まれる高周波（0.5 Hz-2.5 Hz）の変動は、筋の粘弾性に依存して決まる関節インピーダンスや多リンク機構としての骨格系の機械力学的特性といった個々の被験者の身体特性を強く反映し、個人特異性が高いことを明らかにした。これらの結果は、東大高木・中村グループ、
OIST 銅谷グループが構築する大規模並列モデルが実現すべき神経系による身体運動制御のマクロレベルの特性として、制御の間欠性が重要な役割を果たしていることを示唆している。
また、神経疾患に起因する立位姿勢の不安定化を個々の患者毎に定量化するテーラーメイド
シミュレーションの実施と臨床応用において、脳・神経系のみならず、身体機械力学系（筋骨格系）の個別的詳細モデルを構築することも、制御系としての脳・神経系のモデル化と同様に重要であることを意味している。
２） 成果②（発表論文２）
平成 25 年度において、本研究課題の研究協力者である国立病院機構刀根山病院・佐古田院長とともに、すくみ足症状の無いパーキンソン病患者およびすくみ足症状はあるがパーキンソン病ではない患者の歩行運動計測データの解析と計測された運動を再現する単純な数
理モデル化とその数値シミュレーションを実施した。その結果、すくみ足症状を呈する患者の歩行では、平均的な振舞いとしては左右脚運動の相対位相を１ステップで（ある意味強引に）平均値に修正するタイプの左右脚協調が見られる反面、協調制御の神経指令はステップ毎に大きく揺らいでいる（位相修正量の確率的変動が大きい）ことが明らかになった。この結果は、パーキンソン病における肢間協調の崩壊において、運動指令の確率的なゆらぎの存在が重要な役割を果たしていることを示唆している。このことは、歩行運動リズムのゆらぎの生成メカニズムを明らかにすることの重要性を示唆すると同時に、健常者における歩行運
動リズムの揺らぎの生成に、静止立位姿勢の制御と類似した間欠制御が関与する可能性が高いと考えられた。そこで、この発表論文２の研究では、間欠制御による歩行運動の動的数理モデルの構築とその数値シミュレーションを実施した。その結果、ヒト二足歩行運動においても、関節のスティフネスを小さく抑えた状況下で（すなわち、運動に関わる関節がコンプ
ライアントである状況下で）あっても、間欠制御によって、二足歩行がロバストに安定化できることを明らかにした。
これら 2 つの成果は、マクロレベルにおけるヒト立位姿勢制御の特性と制御様式を理論的・定量的に明らかにするものであり、本プロジェクトで開発中のマルチスケールなヒト身体運動制御系が再現すべき現象とその背後にある制御メカニズムを提供すると考えられる。
発表論文リスト
1. Yamamoto T, Smith CE, Suzuki Y, Kiyono K, Tanahashi T, Sakoda S, Morasso P, Nomura T. Universal andindividual characteristics of postural sway during quiet standing in healthy young adults. Physiol Rep. 2015
Mar;3(3). pii: e12329. doi: 10.14814/phy2.12329.
2. Fu C, Suzuki Y, Kiyono K, Morasso P, Nomura T. An intermittent control model of flexible human gait using astable manifold of saddle-type unstable limit cycle dynamics. J R Soc Interface. 2014 Dec 6;11(101):20140958.
doi: 10.1098/rsif.2014.0958.

35