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

Ⅲ－１

高木 周（東京大学）

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

Ⅲ－１－１

レーションに向けた血栓シミュレータと心臓シミュレータの連成及び筋繊維の集合体としての骨格筋のシミュレーション」の研究開発を実施し、研究開発を統括する。
また、プロジェクトを推進する上で、関連する研究者と必要な協議等を行い、予測医療に向けた階
層統合シミュレーション」の研究を統括する。
平成 25 年度は、心筋梗塞シミュレーションに向けた血栓シミュレータと心臓シミュレータの連成として、本課題で開発を進めている血栓症シミュレータを、マルチスケールマルチフィジックス
心臓シミュレータ UT-Heart の冠循環モデルと連成させることを目指す。具体的には、心臓シミュ
レータの計算結果として得られる時間とともに径が変化する冠循環の血管に対して、血小板粘着か
ら血栓形成のプロセスを開発中の血栓シミュレータにより実施する。
、また、抗血小板薬として
P2Y12 の阻害薬と、GPIIb/IIIa の阻害薬に関するシミュレーションを実施し、それぞれの阻害薬の薬効の違いについて調べる。
筋繊維の集合体としての骨格筋シミュレーションの研究課題に関しては、前年度開発された、運動ニューロンからのシグナルにより筋繊維が収縮し、その集合体として骨格筋全体が収縮するモデ
ルを拡張し、UT-Heart に用いられているサルコメア力学モデル、すなわち状態遷移理論に基づきアクチン分子とミオシン分子のクロスブリッジを扱うモデルをモンテカルロ計算で行う手法を導入する。これにより、脳からのスパイクシグナルによりタンパク質のクロスブリッジレベルから筋繊維の収縮、さらに筋肉全体の収縮そしてパーキンソン病特有の振戦・固縮といった病態が再現可能となる。このモデルを用いて、平成 24 年度に並列化を進めた計算コードを改良して、「京」上で並列計算を実施する。

1

Ⅲ－１－２

実施内容（成果）

Ⅲ－１－２－１

P2Y12 の阻害薬と GPIIb/IIIa の阻害薬に関するシミュレーション

１）はじめに血小板血栓が形成されていく過程の中で、血管壁に粘着した血小板が活性化反応を伴い血小板同士が凝集していく段階では、ADP 刺激による血小板の活性化が重要な役割を果たすことが知られている。活性化した血小板からは濃染顆粒に含まれた ADP の放出が起こり、自らの活性化を増幅するとともに、近傍にある血小板に対しても活性化の連鎖を誘発すると考えられており、これはADP ポジティブ・フィードバック機構と呼ばれている。本研究では、血管壁に粘着した血小板について、血小板表面に配置した ADP 受容体である P2Y12 や P2Y1 の活性化を契機として、血小板内の活性化反応によって誘発されるインテグリン GPIIb/IIIa の活性化、および血小板から ADP
放出する動作を模擬する計算モデルについて検討する。

ADP放出(濃染顆粒)

ADP

ADP
ADP

ADP
ADP

ADP刺激(結合)活性化反応： Platelet
自らの活性化への寄与

他の血小板活性化への寄与： Receptor

Fig. 1.1: ADP ポジティブ・フィードバックと血小板活性化連鎖のイメージ
このモデルを用いて、活性化惹起物質である ADP や血小板架橋タンパク質であるフォンビレブラ
ント因子（以後、VWF と略記）やフィブリノゲン（以後、FIB と略記）が移流拡散して濃度変化する環境下で、血小板表面の P2Y12 や GPIIb/IIIa の活性化が経時変化していく過程について考察する。さらに、ADP 結合を阻害するクロピドグレル（以後、CLO と略記）、GPIIb/IIIa と結合して架橋タンパク質との結合を阻害する GPIIb/IIIa 阻害薬（以後、IGP と略記）を課した場合でのGPIIb/IIIa の活性化の違いについてもみていく。
２）血小板表面のレセプタとインテグリンの活性化モデル
対象とするレセプタは、ADP と結合する P2Y12と P2Y1、および血小板同士の架橋タンパク質

VWF・FIB と結合する活性型 GPIIb/IIIa である。レセプタ P2Y12・P2Y1 がリガンド ADP と結合することによって活性化し、血小板内部の活性化シグナリングの伝達により Ca2+ の増加を伴って、インテグリン GPIIb/IIIa が架橋タンパク質と結合が可能な活性化状態に推移するモデルを考える。
本研究では、レセプタ近傍のリガンド濃度に応じて結合判定する確率を考えることで、リガンドとの結合およびレセプタの活性化を決定するようなモデル化を行う。
(i) リガンド結合によるレセプタ活性化モデル
レセプタとリガンドが結合する確率を定式化し、それに基づいてレセプタの活性・非活性の確率モデルを示す。
①リガンド結合の確率
一つのレセプタがリガンドと結合することを、次のような確率モデルで考えることにする。単位2

体積 dV の中に N 個のリガンド分子があり、ある時刻の瞬間に dV に比べて十分小さい微小体積
δv (≪ dV) の中に n 個のリガンド分子が存在する確率を考える。n 個の分子が体積 δv に入り、残りの (N - n) 個が dV にある確率は、次の二項分布に従う。
nv  v  pN , n  N C n  1 
 dV   dV 

N n

(1.1)

δv がレセプタの結合部分の大きさに相当するものとし、そこに１個のリガンド分子が入ることをレセプタとの結合だと思うことにすると、分子の総数 N が十分大きい場合にはポアソン分布の近似として考えることができる。
 v  v 
pN   N   exp N 

dV


 dV 
(1.2)

N は単位体積 dV の中にある分子数であるから濃度と比例した変数である。確率の正規性に留意す
ると、レセプタとリガンドとが結合する確率密度関数は、
  2 x  exp x px    0

x  0,x  0  v dV(1.3)となり、この場合の確率分布関数は、F x   1  1  x  exp x 

(1.4)

として得られる。レセプタとリガンドとの結合は、この確率に従ってレセプタ近傍のリガンドの濃度に応じて生起するものと考える。
②リガンドとの結合による活性化の生起
一つの血小板表面に NR 個のレセプタが存在することとし、その中の第 m 番目のレセプタがリガンドと結合可能な状態である場合を考える。レセプタ m の近傍のリガンド濃度を lm としたと
き、このレセプタがリガンドと結合する確率分布関数を、
 l 
 l 
rm l m   1  1  m   exp  m  ,
 B 
 B 

0  rm  1

(1.5)

として考える。ここで λB は、リガンド濃度に対するレセプタの結合化応答の時定数である。一様
乱数 R[0,1] を発生させて、その乱数値 R と rm とを比べて、rm ≧ R ならば、そのレセプタ m はリガンドと結合することにする。P2Y12・P2Y1 は ADP と結合すると活性化状態に遷移するが、この結合は永続的にではなく、時間経過とともに結合が解放される。
(ii) シグナリングによるインテグリン活性化モデル
リガンド結合でのレセプタ活性化によるシグナリングの生起と、インテグリンの活性化についてのモデルを示す。
①レセプタ活性化で生じるシグナリング
血小板表面上のある単位面積 dSn (= dx×dy) において、Np 個のレセプタがあるとする。その中で、リガンドとの結合で活性化しているレセプタの活性化率（数密度）を xn として次のように定
義する。
1
xn 
Np

Np



idSn

1 act pi , p i  0 otherwise(1.6)

この単位面積 dSn の中で、レセプタ活性化率 xn によってシグナリング物質濃度 sn に増幅されるものと考える。

3 x 
s n  1  exp  n  ,
 S 

0  s n  1

(1.7)

②インテグリン活性化の生起
一つの血小板表面に NG 個のインテグリンが存在することとし、その中の第 m 番目のインテグ
リンの活性化について考える。このインテグリン m が含まれる単位面積 dSm におけるシグナリン
グ物質濃度を sm としたとき、当該インテグリンが活性化する確率分布関数を以下のように考える。
 s u m  1  exp  m  ,
 A 

0  u m  1

(1.8)

ここで λA は、物質濃度に対するインテグリンの活性化応答の時定数である。一様乱数 R[0,1]を発生させて、その乱数値 R と um とを比べて、um ≧ R ならば、そのインテグリン m を活性型に遷移することにする。GPIIb/IIIa は血小板の活性化に伴い活性型への状態に遷移するが、これは永続的にではなく、時間経過とともに非活性型へ遷移するものとする。
(iii) ADP 刺激による血小板内の反応モデル
ここで考えている ADP 刺激による血小板での反応モデルを Fig. 2.1 に示す。
ADP

①

③

P2Y1

②

CLO

⑤

P2Y12
SIGi

VWF

④

FIB

⑨

(Secretion)

ADP

⑩

⑥

⑪

GP IIb/IIIa

⑦

SIGj

IGP

⑧
血小板

Ca2+

血管壁

Fig. 2.1: 血小板内の反応モデル
３）血管壁に粘着した血小板活性化の計算
ADP や架橋タンパク質が移流拡散して濃度変化する環境下で、血小板表面に配置した P2Y12、
P2Y1 および GPIIb/IIIa の活性化が経時変化していく過程について、次の観点から見ていくことにする。
・ADP との結合による P2Y12・P2Y1 活性化と GP IIbIIIa 活性化の伝播過程
・活性化した GPIIb/IIIa が VWF・FIB と結合していく過程
・CLO によって P2Y12 がブロックされた状態での GPIIb/IIIa 活性化の比較・活性化した GPIIb/IIIa が IGP によってブロックされていく過程(i) 血小板①粘着血小板の形状
血管壁に粘着した血小板は、扁平した楕円体の半分の形状で変形しないものとする。大きさは、直径 2.0 μm、血管壁からの高さ 0.5 μm とする。
a=2.0µmb=0.5µm
血管壁

Fig. 3.1: 粘着血小板形状
4

②血小板表面での各タンパク質の分布
血小板表面上の P2Y12、P2Y1、GPIIb/IIIa の配置は、それぞれは血小板表面にランダムに分布しているものとする。設定する各タンパク質の数は、P2Y12 = 300 個、P2Y1 = 100 個、GPIIb/IIIa =10,000 個とする。
③ ADP の放出血小板が活性化することに伴い濃染顆粒に含まれた ADP が放出される。通常、一つの血小板に含まれる濃染顆粒の個数は 5 個程度であることから、ADP の放出は 5 回とし、放出する ADP 量
は 653 mM とする。放出の場所は、血小板の頭頂部 1 cell 分に放出量の ADP 濃度を印加する。
放出の契機は、血小板内の Ca2+ 貯蔵量が閾値に達した時点で ADP 放出が起こるものとし、ADP
放出後は、Ca2+ 貯蔵量を初期状態にリセットする。
(ii) 計算条件①計算サイズ計算の空間領域として、x 方向と y 方向に 10 μm (= Lx = Ly) をとり分割幅を 0.1 μm (= dx = dy)、z 方向に 2 μm (= Lz) をとり分割幅を 0.1 μm (= dz) とした。時間については、計算時間は 15 ms (= t)、
時間の分割幅は 15 μs (= dt) とした。粘着血小板の配置は、9 個の血小板を Fig. 3.2 のように配置
した。

ADP

[µm]

[µm]

Fig. 3.2: 粘着血小板の配置②濃度分布の計算
流れによって輸送される物質 ADP・VWF・FIB・IGP のそれぞれの濃度分布は、流速 ui （血小板表面近傍の流速）を一定値とした移流拡散として計算する。ADP 濃度分布については、血小板活性化の条件に基づき ADP の放出が発生する。
 
  2  
 ADPt , x    x   ADP  D  u iii x 2 xi  t i  
  2  
 FLDt , x    x   0
  D
FLD  VWF , FIB , IGP  u iii x 2 xi  t i 

(1.9)

ここで、xi は座標を示し、 は血小板の内外を示す指標である。
0 inside of platelet 1 otherwise xi   (1.10)

各物質 ADP・VWF・FIB・IGP の濃度について計算条件として、 x, y, z の各方向について一定値の移流速度として、ux = 0.1 μm/s、uy = uz = 0 μm/s とする。物質の拡散係数は 100 μm2/s とし、それ5

ぞれの初期分布は x = 0 での濃度を、[ADP] = 0.1 mM、[VWF] = 0.1 mM、[FIB] = 0.1mM、[IGP] = 0.01
mM として設定する。
(iii) 計算結果例① P2Y12 活性化
ADP が移流拡散していく時間経過の中で、血小板表面の P2Y12 が ADP と結合して非活性状態（紺）から活性状態（赤）に変化する過程を示している。
[mM]

[mM]

(a) t = 0 msec(b) t = 1.0 msec
[mM]

[mM]

(c) t = 2.0 msec

(d)t = 3.0 msec

Fig. 3.3: P2Y12 活性化の経時変化② GPIIb/IIIa 活性化
P2Y12 と P2Y1 の活性化に伴って、GPIIb/IIIa が非活性状態（青）から活性化状態（赤）に遷移していく過程を示している。

(a) t = 0 ms(b) t = 1.0 ms

6

(c) t = 2.0 ms(d) t = 3.0 ms
Fig. 3.4: GPIIb/IIIa 活性化の経時変化③CLO で P2Y12 をブロックした場合の P2Y12 と GPIIb/IIIa の活性化
初期状態として、CLO（緑）によって 10 % の数だけ P2Y12 をブロックした条件で、P2Y12 お
よび GPIIb/IIIa の活性化状態について、CLO の有無による比較を示している。

(a) P2Y12 (t = 5.0 ms)(b) GPIIb/IIIa (t = 5.0 ms)(c) P2Y12 with CLO (t = 5.0 ms)(d) GPIIb/IIIa with CLO (t = 5.0 ms)

Fig. 3.5: CLO の有無による P2Y12 と GPIIb/IIIa 活性化の比較④GPIIb/IIIa 活性化と VWF・FIB との結合活性化した GPIIb/IIIa （赤）に VWF （黄）と FIB （緑）とが結合していく時間変化の経過を示している。

(a) t = 0 ms(b) t = 1.0 ms
7

(c) t = 2.0 ms(d) t = 3.0 ms

Fig. 3.6: 活性化 GPIIb/IIIa に結合する VWF と FIB の経時変化⑤IGP を付加した場合の GPIIb/IIIa 活性化
IGP は活性化した GPIIb/IIIa と結合することで、VWF・FIB との結合を阻害する。ここでは、IGP（紫）を付加した環境下での、VWF（黄）・FIB（緑）の結合過程の時間変化を示している。

(a) t = 0 ms(b) t = 1.0 ms(c) t = 2.0 ms(d) t = 3.0 ms

Fig. 3.7: IGP を付加した場合での VWF・FIB 結合の経時変化
４）まとめ
血小板血栓が形成されていく過程の中で、ADP 刺激による血小板活性化に着目し、血管壁に粘着した血小板表面での、P2Y12・P2Y1 が ADP と結合することを契機とした GPIIb/IIIa 活性化についての計算モデルの検討を進めてきた。確率過程をベースにしたリガンド結合と活性化する血小板モデルを考え、このモデルを用いて、ADP や架橋タンパク質が移流拡散して濃度変化する環境
下で、血小板表面に配置した P2Y12 や GPIIb/IIIa の活性化が経時変化していく過程についての確認を行った。今回の流れ計算は簡易な移流拡散の計算で行ったが、赤血球や血小板が存在する血管内での血流環境が、血小板から放出された ADP と合わせて、血管壁近傍での ADP 濃度分布形成にどのように関係するかについての検討を行うため、血流計算との連成を進めていく。また検証実験との比較が行えるようにするため、レセプタ反応の応答時間など実際の血小板活性化に対応した計算条件のチューニングを行っていく。さらに、 ADP 阻害薬と GPIIb/IIIa 阻害薬における薬効の比較評価のための検討を進めていく。
参考文献(1) 池田康夫、丸山征郎（編）、 血小板生物学、 メディカルレビュー社 (2004)。
8

(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。
(6) 中川正雄、 真壁利明、確率過程、培風館 (2009)。

9Ⅲ－１－２－２ 冠循環末梢抵抗の推定に向けた微小循環系のモデリング１） はじめに本テーマでは，UT-Heart の冠循環系の末梢抵抗のモデリングに関連して，末梢血管網の三次元数値シミュレーションにより冠循環における末梢血管系の流れを詳細に解析，特に非定常に変化する
圧力損失の支配的因子を理解することにより，微小循環器系の末梢抵抗モデルの精度を向上させることを目的としている．本年度はそのモデリングを開始した．以下，その内容について説明する．

２） 冠循環毛細血管網に対するモデリング
冠循環の流れの特徴として，拍動に合わせて流れの様子そのものが大きく変化することが挙げられる．つまり，冠動脈においては心臓が収縮する周期（心臓から血液が拍出する周期）で流量が少なく，拡張期では流量が多くなる．逆に冠静脈では，収縮期で流量が多く，拡張期では少ないことが知られている．また，その構造は階層的な樹状構造を持ち，心外膜側から太い血管が細い血管へと分岐・細分化しながら心内膜側へと到達していることなどが知られており[1]，このような構造が，心筋への血液供給に対して生理学的・力学的に重要な機能を有していると考えられている．冠循環での圧力損失は，血管の湾曲，分岐，断面積の変化等の単一液体相流れでも起こり得るものだけでなく，赤血球（狭い流路断面積に対して相対的に大きな弾性体）などが存在することによる複雑さも加わり，推算が難しいものとなっている．
冠循環の構造は，心筋の一部を取り出した後に染色などの処理を施すことで，画像として得ている研究成果も存在する[2]．本来そのような画像から毛細血管を抽出することで流路を構築し，数値
シミュレーションを実行することが望ましいと考えられるが，抽出した画像データの精度・品質等の点から直接格子を作成することが現状では難しい．ここでは，久田らが UT-Heart の末梢血管径のモデルとして採用している１次元毛細血管モデルである．解剖学的知見から得られた血管径・長さ（分岐から分岐までの長さ）などの統計的な分布を再現するような仮想的な毛細血管網流路のデータを利用することを考えた．そのため，まずは，UT-Heart で用いられている１次元血管網モデルのデータではなく，１次元モデルの可視化用に試験的に作った血管網データを利用した．この１次元データを基に，３次元の血管網を再構築し，３次元流路を対象とした流体シミュレーションを実行
した．本研究で用いた毛細血管網流路を Fig. 1 に示す．流体は図中左端から流入し，右端から流出する．それ以外は全て閉じた流路となっており，図中上下方向や紙面垂直方向に流入流出境界は存在しない．この流路形状では，血管直径は 4 µm から 21 µm，直線形状の一断片の長さは 30 µm
から 180 µm となっており，二次元的に分岐・合流するような形状をしている．また，屈折部や分岐・合流部は計算の安定性等の観点から滑らかな曲線形状として構築しており（Fig. 2）．実際の冠循環で見られるような洞状の部分や三次元的な分岐構造[3]は持たせていない. 屈折部を接続する曲面は, 半径を変化させつつ移動する球の軌跡として与えられる(Fig. 3)． 球の中心は, 両端の血管の中心線に接する円弧上を移動し, 半径は両端での血管径に滑らかに接続するように, 弧長の 3 次関数として与えた．これにより，血管の中心を通る円弧の半径を調節することで, 両端の血管と干渉しない範囲で接続の滑らかさを調節した．本稿では屈折部を滑らかに構築した血管網を用いている. また分岐部及び合流部はさらにこれを接続する 3 本の血管から 2 本を結ぶ 3 通りに対してそれぞれ構築し, 内側半分を除去して上下を平面で接続することで構築した. 血管網には平面上では血管が交差する部分が存在するが, こうした部分は本来の 3 次元形状では接合していないと考えられるので交差を回避するように構築した.
本解析には，野田らが開発してきたボクセル型データに適した計算手法を利用した．この手法は，非圧縮性流体に対する解析手法として HSMAC 法をベースに，直交格子上で AOF (Area fraction Of
Fluid) 値および VOF (Volume fraction Of Fluid) 値を用いて流体・個体を判別することで複雑流路を表現する FAVOR (Fractional Area/Volume Obstacle Representation) 法を用いている[3]．本稿では，複雑な分岐構造をもつ血管網に対するシミュレーションとして，まずは，血液を非圧縮性の単一液相として扱い，赤血球の存在などは考慮しない計算を行った．．密度は 1.06×10-3 kg/m3，粘性係数は 4.70×10-3 Pa･s とした．流入境界においてはポアズイユ分布とし，拍動時の最大流量として
およそ 4×10-12 m3/s を想定することで得られる最大流入速度 4.12 cm/s を速度に対する入口境界条件として課した．なお，この条件におけるレイノルズ数は 0.146 である。
10

結果の一例として上記条件に対する速度分布を Fig. 3 に示す．Fig. 4 中上図は全流速を示しており，黒色が濃い領域ほど流れは速い．また，下図は代表的な領域の速度ベクトル分布を表示している．なお B 図は A 図に比べて 2 倍で表示されている（矢印長さが同じ場合 A 図の流速は B 図の流
速の 2 倍となっている）．解析結果から，流入部付近（細動脈に相当）ならびに流出部付近（細静脈に相当）の血管径が比較的太い領域において流速が高い．一方で，多くの分岐が存在する中央部付近（毛細血管に相当）ではあまり高い流速は得られていない．速度ベクトル分布からは分岐部における流れは流路の湾曲に沿って滑らかに流れている様子が確認できる．また，全体的に左から右へ流れる速度が支配的であり，図中で上下に伸びている毛細血管部の流速が特に小さいことがわかる．この傾向は分岐部のベクトル分布からも確認できる．
今後は赤血球の変形流動構造も考慮に入れたより現実的な微小血管網流路のモデリングに取り組み，計算結果の詳細な考察を通じて，圧力損失に及ぼす血管径の変化や分岐の影響，また，拍動によって流量が大きく変動した際の流れの様子や流量と圧力損失との関係などについてモデリン
グを進めて行く．

Fig.1 Constructed microvasculature for the numerical simulation

Fig.2 An example of branch structure.

Fig.3 An example of bending structure.

11

Fig.4 Velocity distribution in microvasculature

参考文献[1] Matsumoto, T., et al., Fluid Dynamics Research, 37(2005), 60-81.
[2] Kaneko, N., et al., Am. J. Physiol. Heart. Circ. Physiol, 300(2011), 754-761.
[3] Hirt, C.W., J. Wind Eng. Ind. Aerodyn, 46 & 47(1993), 327-338.

12

Ⅲ－１－２－３ マルチスケール骨格筋シミュレータの開発（脊髄神経ネットワークモデル）
１）はじめに本アプリケーションでは、主に単関節運動時における脊髄神経系の挙動と筋収縮の挙動を連成させた詳細解析を可能とすることにより、安静時振戦・筋固縮などのパーキンソン病症状の原因究明に貢献することを目的とする。そのために、脊髄神経系のニューロンモデルならびにそれらのネッ
トワークモデルの構築、筋骨格系の有限要素解析モデルの構築、ならびにそれらを連成させる解析モデルの構築に取り組んでいる。
２）ニューロン発火の数理モデル
本研究では、個々の神経細胞内における電気化学方程式とそれらを結合させたネットワークで形成される電気シグナル、さらにそのシグナルにより筋繊維が収縮し、結果として筋肉全体が 3 次元的に収縮するまでをシミュレーションの対象とする。この際、脊髄神経系の各ニューロンの状態は、神経細胞モデルとして Hodgkin-Huxley 方程式に代表されるコンダクタンスベースの数理モデルに基づき、膜電位 V に関する常微分方程式として記述する(1)(2)。本研究で用いたモデルでは、3 種類のイオンチャネルの挙動を考慮しており、その場合、1 つのニューロンの膜電位の経時変化について以下のような連立常微分方程式で記述される。

CdV
  I ion  I Leakdt
I L  g L VS  EL 

(2.1)
(2.2)

I ion  I Na  I Kf  I Ks

 g Na m 3 hVS  E Na   g Kf n 4 VS  E Kf   g Ks VS  E Ks ds
  s 1  s    s sdt

s  m, h, n, q 

(2.3)

(2.4)

ここで、gNa、gKf、gKs は ナトリウムチャネルと 2 種のカリウムチャネルに関する電気伝導度を表すモデル定数であり、ENa、EKf、EKs はそれぞれの平衡電位である。また、gL・EL は Leak 電流に対応する電気伝導度・平衡電位を意味し、生理学的には先に挙げた 2 種類以外のイオン流出入による効果を総じて表したものである。
２）脊髄ネットワークの概要と運動ニューロンプールモデル
脊髄内では、同種の神経細胞が数百～千個程度集合したニューロンプールを形成し、それらの間で相互に興奮性・抑制性のシナプス結合を形成することで、筋肉（筋肉群）の活動を支配するネッ

Fig. 2.1: 脊髄神経ネットワークの概要図．
13

トワークを構成している。筋肉を直接支配するα運動ニューロンとその活動に大きく影響する主な
介在ニューロンから成る脊髄神経系の概要図を Fig. 2.1 に示す。図は一対の主働筋・拮抗筋から成る仮想的な関節を想定し、それに関連する脊髄内のニューロンプールの結合を概略図として図示したものである。前年度までに、関節運動に最も重要なα運動ニューロンと主要介在ニューロンの一つであるレンショウ細胞に関するモデリングを中心に行った結果、Recruitment および Rate codingと呼ばれる連続的な筋力調節機構を正確に再現するに至った。また、α運動ニューロン・レンショ
ウ細胞間結合の重要な機能である反回抑制機構の再現、また、α運動ニューロンとレンショウ細胞の発火頻度がほぼ線形な関係を有していることなどを確認し、解剖等から得られている重要な知見を再現するモデルであることを確認した(2)。今年度は、静的・動的γ運動ニューロン、Ia 介在ニューロン、Ib 介在ニューロンを加え、筋肉 1 つに対する脊髄神経系モデルの構築を行った。さらに、筋肉から脊髄、脳へのフィードバック機能を担う、筋紡錘・ゴルジ腱器官の 2 つの固有受容器モデルを導入し、脊髄・骨格筋連成モデルの完結に必要な部位の構築を行った。静的・動的γ運動ニュー
ロンは筋紡錘の感度調節機能を担う。この働きにより筋紡錘は筋肉の受動的な伸長に対して常にフィードバックを脊髄・脳へ送ることが出来る。Ia 介在ニューロンは、筋紡錘からの信号を中継し、伸長反射において重要な役割を担う他、脳からの運動指令を受け円滑な関節運動を可能とする拮抗抑制機能も担う。Ib 介在ニューロンはゴルジ腱器官からの信号を中継し、過剰な筋活動を抑制する機能（逆伸長反射）を担う他、ネガティブフィードバックとして機能し筋張力を一定に保ち関節を固定する機能も担う。筋紡錘は紡錘形をした受容器であり筋組織の一部が分化したものである。筋
肉内に筋線維と平行に一筋肉当たり数十～百個程度存在し、筋肉の受動的な伸長に応答し、長さと長さの時間変化に応じたスパイク信号をα運動ニューロン、Ia 介在ニューロン、脳（視床）へ返す(3)。ゴルジ腱器官は、筋線維の腱への移行部に一筋肉当たり最大 100 個程度存在し、筋線維の能動的な収縮に応答し、Ib 介在ニューロンにスパイク信号を返す(4)。筋紡錘が主に受動的な長さ変化に関するフィードバック機能を担うのに対し、ゴルジ腱器官が主に能動的な長さ変化に対するフィードバック機能を担うという特徴がある。ここでは、これらの機能を非線形ばねモデルを用いて数理モデル化した Mileusnic らによるモデル(5)(6) を採用した。
３）シナプス結合の数理モデル
上記のニューロン単体あるいはニューロンプールについての数理モデルに加え、神経ネットワー
クモデルを構築するためにはシナプス結合に関する数理モデルが必要である。シナプスとはニューロン間で信号が伝達される箇所であり、シナプス結合においては、① シナプス前ニューロンにおいて活動電位が生成され、軸索末端まで伝達される。
② シナプス前ニューロンの軸索末端からシナプス間隙に神経伝達物質が放出される。
③ シナプス後ニューロンの受容体（樹状突起あるいは細胞体）に神経伝達物質が結合する。
④ 神経伝達物質の結合によってシナプス後ニューロン内においてイオンチャネルが開き、その結果活動電位が生成される。
という一連の流れによって活動電位が伝達される。ここでは、上記の一連の流れをシナプス前ニューロンの発火に対応してシナプス後ニューロンに刺激電流が入力されることと同等とみなしモデル化する(7)(8)。この電流をここでは、シナプス電流 Isyn と呼ぶ。つまり、シナプス結合がある
場合、シナプス後ニューロン j の膜電位に関する常微分方程式は、以下のように修正される。

Cj

dV j
dt

  I ion  I Leak  I syn, j(2.5)

Isyn, j は神経伝達物質の移動を電流が流れることと等価とみなすことで以下のようにモデル化する。
シナプス前ニューロンプールの i ニューロンから、シナプス後ニューロンプールの j ニューロンへの電流は、

I syn, i j  g i j t Esyn, i j  V j (2.6)とモデル化されるので、シナプス後ニューロン j に流れ込む総電流値は

I syn, j   g i j t Esyn, i j  V j 
Ni

(2.7)

i 1

14

となる。ここで、gi→j(t) はシナプス結合間の電気伝導度を意味し、後述するように、シナプス前ニューロンの活動電位に依存して時間的に変化する。Ei→j(t) は reversal potential と呼ばれ、対象としているシナプス結合が興奮性なら正（70 mV 程度）、抑制性なら負（-16 mV 程度）の値を取るモデル定数である。シナプス結合間の電気伝導度 gi→j(t) は以下のようにモデル化する。神経伝達物
質とシナプス後ニューロン受容体との結合を化学反応として記述することを考えたとき、最も簡略
すれば、


RT
RT *




(2.8)

のように表すことができる。ここで、R は未結合の受容体、T は神経伝達物質、RT* は既結合の受容体、a および b はそれぞれの反応速度定数である。この関係を微分方程式として記述すれば、dr  T 1  r    r
dt

(2.9)

となる。[ ] は濃度であることを表し、r は神経伝達物質受容体のうち既結合の割合である（したがって、r = [RT*]/([R]+[RT*])）。ここで、神経伝達物質の放出開始・終了が瞬時に行われると仮定する。つまり、シナプス前ニューロンにおける活動電位生成時刻を t0、神経伝達物質の放出終了時刻
を t1 としたとき、

t0  t  t1 : [T]  Tmax
else
: [T ]  0(2.10)とすれば、2 状態について以下のように解析的に導くことができる。
During a pulse (t0 < t < t1)、[T] = Tmax  t  t 0  r t  t 0   r  r t 0   r exp

r


r 

Tmax

Tmax  ,r 

(2.11)

1

Tmax  

After a pulse (t1 < t)、[T] = 0

r t  t1   r t1 exp  t  t0 (2.12)

この r を用いて、

g i j t   g max  r t 

(2.13)

としてシナプス電流を求める。gmax は対象とするシナプス結合間の電気伝導度の最大値を表すモデ
ル定数である。
４）解析結果
脊髄神経系のニューロンモデル・シナプス結合モデル・ネットワークモデルと筋肉内の固有受容器に関するモデルを結合することで、脳・脊髄・骨格筋を統合した解析を行うことが可能となり、脳で生成された運動指令がどのように伝達し、あるいは伝達過程で変化し、最終的な関節運動に至っているかを検討することが可能となる。一対の筋肉から成る仮想的な関節を表すモデルに対し、脳からの運動指令に対応する興奮性刺激を各ニューロンプールに入力することで、関節運動やその時の各ニューロンの活動を解析した。一例として、このモデルに主働筋・拮抗筋に逆位相で興奮性刺激を入力した場合に発揮される筋力を Fig. 2.2 示す。各種介在ニューロンが存在することでα運動ニューロンのみを考えた場合に比べてα運動ニューロンの発火頻度が低下するため（介在ニューロンはα運動ニューロンに対して抑制性のシナプス結合を形成している）、それに応じて発揮され
15

る筋力も下がることがわかる。ただし、発揮される筋力の経時変化については、入力信号の経時変化に対応するものとなっており、振動周期が脊髄の神経ネットワーク内で変化するようなことは起こらなかった。今後定量性も含めてさらに慎重な検討を必要とするものの、関節運動の周期などは
α運動ニューロンへの入力信号の周期に対応している可能性が高いと思われ、パーキンソン病の代表的症例の一つである安静時振戦が見られる状況では、（脳あるいは脊髄内に存在する介在ニューロンによって）α運動ニューロンへ到達する前に振動を生みだすような運動指令が作られている可能性が高いと思われる。

Fig. 2.2: 構築した脊髄神経モデルに周期的な興奮性を入力を与えた場合（左）に発揮される筋力（右上：α運動ニューロンのみの場合，右下：介在ニューロン等まで含めた神経回路ネットワー
クの場合）
．
５）まとめ
これまでに詳細なニューロン発火モデル，ニューロンのサイズ効果などを考慮したニューロン
プールモデル，神経伝達物質の放出・回収過程を考慮したシナプス結合モデルを用い，解剖学に基づき脊髄ネットワークを構成し，さらには筋肉からのフィードバックを含めた解析を可能とする数
値シミュレーションモデルを開発した．これにより，脳神経と骨格筋の収縮モデルを繋いだ統合解析が可能となる．今後は，脳モデル・骨格筋モデルとの統合作業を進め，パーキンソン病などで現れる具体的な症状を再現可能とするよう開発を進めて行く予定である．
参考文献(1) Cisi, R. R. L., Kohn, A. F., Simulation system of spinal cord motor nuclei and associated nerves andmuscles, in a Web-based architecture, Journal of Computational Neuroscience, 25(2008), 520-542.
(2) 2012 年度成果報告書．
(3) Boyd, I. A., The isolated mammalian muscle spindle, Trends in Neurosciences, 3(1980), 258-265.
16

(4) Jami, L., Golgi tendon organs in mammalian skeletal muscle: functional properties and central actions,Physiological Reviews, 72(1992), 623-666.
(5) Mileusnic, M. P., Brown, I. E., Lan, N., and Loeb, G. E., Mathematical Models of Proprioceptors. I.
Control and Transduction in the Muscle Spindle, Journal of Neurophysiology, 96(2006), 1772-1788.
(6) Mileusnic, M. P., and Loeb, G. E., Mathematical Models of Proprioceptors. II. Structure and Function of
the Golgi Tendon Organ, Journal of Neurophysiology, 96(2006), 1789-1802.
(7) Destexhe, A., Mainen, Z. F., Sejnowski, T. J., An efficient method for computing synaptic conductances
based on a kineticmodel of receptor-binding, Neural Computation, 6(1994), 14-18.
(8) Destexhe, A., Mainen, Z. F., Sejnowski, T. J., Synthesis of models for excitable membranes, synaptictransmission and neuromodulation using a common kinetic formalism, Journal of Computational
Neuroscience, 1(1994), 195-230.

17

Ⅲ－１－２－４

マルチスケール骨格筋シミュレータの開発（骨格筋 FEM とサルコメア力学モデル）本プロジェクトではコンピュータシミュレーションによる脳神経疾患の運動機能障害の病態予
測と治療支援を目指して、ヒト全身の神経-筋-骨格系の統合シミュレーションシステムの構築を行なっている。本統合シミュレータは、Fig. 3.1 に示すように、脳神経系シミュレータ、脊髄神経系
シミュレータおよび筋骨格系有限要素法（FEM）シミュレータからなるシステムである。

Fig. 3.1: Integrated simulation system of neuro -musculo-skeletal system for human whole body.
本研究ではこのうち、CT や MRI などの医用画像から再構成した筋骨格系の 3 次元形状モデル
を用いて、生体軟組織の変形、骨格筋の筋活動挙動および骨格の関節運動を再現するための筋骨格
系 FEM シミュレータ(1) (2)の開発を行っている。 Fig.3.2 は本シミュレータによるパーキンソン病振戦（ふるえ）のシミュレーション例である (3)。本解析では、下腿三頭筋、前頸骨筋および下腿の骨格からなる足関節 3 次元モデルを用いて、拮抗する筋への交互の筋活動の入力により、足関節まわりの繰り返し運動、すなわち振戦が再現されている。本解析では筋収縮モデルに現象論的モデル(4)（Hill model）を採用し、筋全体が一様に収縮することを仮定した。

Fig. 3.2: Simulation results of the ankle joint movement generated by the activation of the triceps suraemuscle and tibias anterior from time 0 to 0.2 s.
18

筋収縮のための脳からの司令は、脊髄内で感覚受容器からの情報などと統合され、最終的に脊髄前角にあるα運動ニューロンを通して筋線維に伝達される。各α運動ニューロンはそれが支配する
筋線維群とともに運動単位を形成しており、筋活動はこの運動単位の活動の和として現れる。骨格
筋線維は直径 1 μm ほどの筋原線維の並列集合であり、筋原線維はサルコメアという繰り返し構造から成る。サルコメア内ではアクチンフィラメントとミオシンフィラメントが規則的に配列し、収縮力の担い手と成っている。そこで平成 25 年度は、骨格筋の収縮力発揮の最小単位であるサルコ
メア内部のタンパク質分子の確率的な振る舞いを考慮した、新たなミクロスコピック筋収縮モデルの導入を行った。本収縮モデルはサルコメアを構成するアクチンフィラメントとミオシンフィラメ
ントに対して状態遷移モデルを定義し、それぞれの状態遷移をモンテカルロ法でシミュレートすることでタンパク質フィラメントの挙動を再現し、筋収縮力を算出する手法(5) (6)である。現在、モン
テカルロ法による筋収縮力の算出部分の開発は終了し、今後、筋収縮力とマクロな骨格筋有限要素モデルとのカップリングを行い、タンパク質分子の確率的な挙動から 3 次元骨格筋の収縮や骨格の関節運動を再現する予定である。
１）骨格筋 FEM シミュレータ
骨格筋や腱などの軟組織は超弾性、非圧縮性、大変形、非線形性などの特性を有する。そこで本研究では、非圧縮性超弾性体の 3 次元混合型 FEM プログラムをベースにシミュレータの開発を行っている(1) (2)。本手法の基礎となるひずみエネルギー関数は次式で与えられる(7)。

k1W  W (C )  ( J  1) 2  ( p  p ) 222k
ここで、右辺第 1 項は等容変形（ C  J  C 、 C は右 Cauchy-Green 変形テンソル）によるひずみエ
2
3

ネルギーの寄与分で、超弾性材料構成式により記述される項である。第 2 項は体積変化（ J  det F 、Fは変形こう配テンソル）によるひずみエネルギーの寄与分である。体積弾性率 k は非圧縮性

を微圧縮性として導入したことによる体積変化を拘束するためのパラメータである。第 3 項は混合法の導入によって付加される項で、要素レベルで導入される圧力変数 p と変位場から計算される
静水圧 p の不つり合いを補正する項である。
超弾性体ではひずみエネルギー関数 W の Green-Lagrange ひずみテンソルの成分 Eij による微分、

Sijpas WEijは第 2 Piola-Kirchhoff 応力テンソルの成分 Sijpas に相当する。右肩の pas は筋収縮によらない受動的な変形の影響を表している。
２）タンパク質分子の確率的な振る舞いを考慮した筋収縮モデル(5) (6)骨格筋の収縮力の発揮メカニズムは次のようである。運動ニューロンからの活動電位が筋線維に伝達すると、筋細胞膜に活動電位が発生し、横行小管を通して筋内部に伝達される。その後、活動電位は横行小管に接する筋小胞体に伝達され、筋小胞体からカルシウムイオンが放出される。カル
シウムイオンはアクチンフィラメントに結合し、これによりアクチンフィラメントが活性化すると、ミオシンフィラメントが活性化され、ミオシンヘッドでは ATP が加水分解される。このエネルギーによりミオシンヘッドが首振り運動を行う。この結果、ミオシンアームが引き伸ばされて収縮力が発生し、アクチンフィラメントは引き込まれ、フィラメント間が収縮する。
このようなタンパク質分子の運動は熱ゆらぎのもと、確率的に進行する。そこで本研究では、アク19チンフィラメントおよびミオシンフィラメントに対してそれぞれ状態遷移モデルを定義し、それぞれの状態遷移を直接的にモンテカルロ法でシミュレートすることでタンパク質フィラメントの挙動を再現し、収縮力を算出する手法(5) (6) を導入した。Fig. 3.3 はアクチンフィラメントの Ca2+ 結
合部位およびミオシンフィラメントのミオシンヘッドの状態遷移モデルである。アクチンフィラメ
ントの Ca2+ 結合部位の状態遷移モデルはカルシウムイオンの結合状態と離脱状態を表すモデルとし、結合部位が 1 つであることを仮定した 2 状態モデルと 2 つであることを仮定した 3 状態モデルを採用した（Fig. 3.3(a)）。また、ミオシンヘッドの状態遷移モデルは、ミオシンヘッドの状態としてアクチンフィラメントと離れていている状態、さらに、結合している状態として、弱結合、強結合および首振り状態に分かれる 4 状態モデルを採用した（Fig. 3.3(b)）。各矢印の方向にはそれぞれの状態への遷移率が設定される。

(a)

(b)

Fig. 3.3: Transition state models of (a) Ca2+ binding site of actin filament and (b) myosin head.

３）アクチンフィラメントの Ca2+結合部位の遷移率 (5) (6)Ca2+ 結合部位の結合状態への遷移率 Kon は、まわりの環境における Ca2+濃度[Ca]に依存する。

Kon  K on  [Ca ]ここで、 K on は遷移状態モデルのパラメータである。また、離脱状態への遷移率 Koff は、その結合部位直下に結合しているミオシンヘッドの存在に依存する。

 K 'off
Koff   K off

if there is a binding MH below
otherwise

ここで、 K off , K 'off は遷移状態モデルのパラメータである。ここでは、 K off  K 'off として設定しているため、結合部位直下でミオシンヘッドとの結合が形成されている場合、カルシウムイオンは解離しにくいモデルとなっている。
４）ミオシンフィラメントのミオシンヘッドの遷移率(5) (6)ミオシンヘッドのアクチンとの結合の遷移率 Knp および弱結合および強結合状態から離脱への遷移率 Kpn は次式で与えられる。

Knp  K np  RA  LA n
Kpn  K pn  n
ここで、 K np , K pn は遷移状態モデルのパラメータ、  RA ,  LA はアクチンフィラメントとミオシン
フィラメントのオーバーラップ状態を表すサルコメア長依存性パラメータ、  n および   n はミ
オシンヘッドの結合状態に起因する協調性を表すパラメータである。n は着目するミオシンヘッドの両隣が非結合状態ならば 0 を、1 つが結合状態なら 1 を、両方結合状態の場合には 2 を代入する。
20

ここでは、  80 と設定するため、両隣のうち 1 つが結合状態にあれば、結合の線維率は 80 倍、離脱への線維率は 1/80 倍となる。また、首振り状態から離脱への遷移率 Kpnr は次式で与えられる。

Kpnr   K pn  n  g xb
ここで、  は 1 より小さい定数で、結合状態のなかで首振り状態の結合が最も強く、結合が外れにくいことを表すパラメータである。g xb はミオシンアームのひずみに応じて変化するパラメータ
で、ミオシンアームが適度に張力を発生している状態では結合を維持する効果を表している。
結合状態間（弱結合－強結合－首振り）の遷移率は ATP の加水分解エネルギーが首振り運動によりミオシンアームの弾性エネルギーに変換される形で消費されるとし、平衡状態での状態分布が内部エネルギーに基づく統計力学的な分布則に沿うように決定する。着目するミオシン 1 つの全内部エネルギーはミオシンヘッド内のエネルギーとミオシンアームの伸びによる弾性ポテンシャル
エネルギー Warm の和で与えられると仮定する。
弱結合、強結合および首振り状態におけるミオシンヘッド内に保存されている内部エネルギーをそれぞれ、Up、Upre および Upos とするとき、弱結合から強結合への遷移率 fapp および強結合から弱結合への遷移率 gapp をそれぞれ次式で定義する。

 U U p fapp  capp exp   prek BT 
gapp  capp
ここで、capp は遷移状態モデルのパラメータ、k B はボルツマン定数、T は温度である。同様に、強結合から首振りへの遷移率 hf および首振りから強結合への遷移率 hb をそれぞれ次式で定義する。

 W ( L  x0 )  U pos 
hf  h0C0 exp   arm

k BT


 W ( L  x0 )  U pre 
hb  h0C0 exp   arm

k BT


ここで、h0 , C0 は遷移状態モデルのパラメータ、L はミオシンアームの伸び、 x0 はミオシンヘッ
ドの首振りによるミオシンアームの伸びである。
本研究では各状態の内部エネルギーを ATP 加水分解エネルギー FATP を基準に、次式のように定義する。

Up  FATP
Upre  FATP  U prepost
Upos   pre FATP
ここでは、弱結合から強結合への遷移においては、 U prepost のみ、わずかに内部エネルギーは増加するとし、また、首振り状態ではほぼエネルギーを使い切り、  pre  100 % のエネルギーのみが残っているとする。
弾性ポテンシャルエネルギー Warm のモデルには、次式の線形モデルと非線形モデル 2 種類を用
いた。
21

1
k xb L22
1
,L  0k xb L22knon  linear
Warm( L)   xb2 (exp( a  L)  aL  1)
,  L1  L  0a k xb
 k xb2
 a 2 (exp( a  L1 )  aL1  1)  2 ( L  L1 )  F1 ( L  L1 ) , L   L1
linear( L) 
Warm

ここで、 k xb ,  , L1 は材料定数、 a, F1 はこれらから求まる定数である。
５）モンテカルロ・シミュレーション
以上の状態遷移モデルを用いて、モンテカルロ・シミュレーションを行う。Fig. 3.4 に示すような、複数の Ca2+ 結合部位とミオシンヘッドを持つ半サルコメアモデルを複数設定し、各 Ca2+ 結合部位およびミオシンヘッドの状態をモンテカルロ・シミュレーションにより、確率的に決定する。

Fig. 3.4: Actin binding sites and Myosin heads of a half sarcomere model used in Monte Carlo simulation.
Fig. 3.5 に 1 つの半サルコメアモデルのモンテカルロ・シミュレーションにおける、Ca2+ が結合状態にある Ca2+ 結合部位の割合の時間変化を示す。Fig. 3.5(a) は 2 状態モデルで [Ca] = 1 M 一定の結果である。定常状態での Ca2+ 結合の割合は 0.407 になるように遷移率は設定されている。カ
ルシウムイオンの結合は確率的に決定されるため、時間変化に伴いばらつきは発生するが、ほぼ定常状態の割合が再現されている。Fig3.5(b) は [Ca] = [0:1.0] とカルシウムイオン濃度を線形に増加させたときの Ca2+ が結合状態にある割合の 2 状態モデルと 3 状態モデルの比較である。3 状態モ

1.0

1.0

0.8

0.8

0.6

0.6

[Ca-on][Ca-on]デルでは [Ca] = 0.1 - 0.2 μM あたりで Ca2+ 結合の割合に急激な上昇が確認できる。

0.4
0.2
0.0
0.0

0.4

2 state model
3 state model

0.2

0.2

0.4

0.6

0.8

0.0
0.0

1.0

0.2

0.4

0.6

Time [s]Time [s]

(a)

(b)

0.8

1.0

Fig. 3.5: Concentration of Ca2+-on during Monte Carlo simulation of Ca2+ binding sites. (a) 2 state model.
[Ca]=1.0 μM is constant during the simulation. (b) Concentration of Ca2+-on for 2 state and 3 state model.
[Ca] is linearly increased during the simulation.

22

Fig. 3.6(a) に 1 つの半サルコメアモデルのモンテカルロ・シミュレーションにおける、ミオシンヘッ
ドの結合状態（弱結合－強結合－首振り）の各状態の割合を示す。本解析では、弱結合：強結合：
首振りの割合は、定常状態において、2:1:1 になるように各遷移率を設定した。モンテカルロ・シ
ミュレーションでは、時間変化に伴いばらつきは発生するが、ほぼ定常状態の割合が再現された。
Fig. 3.6(b) は半サルコメアモデルのサンプル数を 1、8 および 32 と変化したときの首振り状態の割合を示した結果である。サンプル数の増加に伴い、時間変化に伴うばらつきは低減し、理論的な割合に近づくことが確認できる。骨格筋シミュレーションにおける適切なサンプル数の決定は今後の検討課題である。
0.5

1.0
Weak binding
Strong binding
MH rotation

0.4

Concentration

Concentration

0.8
0.6
0.4

0.3
0.2
0.1

0.2
0.0

1
8
32

0

1

2

3

0.0

0

1

2

Time [s]Time [s]

(a)

(b)

3

Fig. 3.6: (a) Concentration of weak binding、 strong binding and myosin head rotation during Monte Carlosimulation of myosin heads. (b) Concentration of myosin head rotation for 1, 8 and 32 half sarcomere models
during Monte Carlo simulation of myosin heads.

６）収縮力の算出（有限要素法とのカップリング）(6)
有限要素シミュレーションでは計算時間の問題から時間増分の幅をモテカルロ・シミュレーショ
ンより大きくとる必要がある。そのため、有限要素シミュレーションの 1 増分ステップ内で、複数のモンテカルロ・シミュレーションを行うことになり、そこでクロスブリッジの生成と消滅が繰り返される可能性がある。そこで、有限要素シミュレーションとモンテカルロ・シミュレーションのカップリングでは、有限要素モデルにおける筋原線維要素のアクティブ応力テンソルによる仕事量がその要素に埋め込まれたミオシンアームによる仕事量の総和に等しくなるようにアクティブ応
力テンソルを定義する。
有限要素シミュレーションにおける時間区分 [ T : T  T ] におけるミオシンアームの伸びの平均
値 L は次式で与えられる

L

t 

T
it
it
 XBD 0 Ls   ( LR  XBD  u  t ) 
T 
it


ここで、 t , T はそれぞれ、モンテカルロ・シミュレーションおよび有限要素シミュレーションの時間増分である。 XBD 0T Ls は前時間区分終了時における短縮速度 u による寄与分、 it LR はミオ
シンヘッドの首振り運動による寄与分、it XBD  u  t は現時間区分における短縮速度 u による寄与分である。it はモンテカルロ・シミュレーションステップを表す。
ミオシンアームによって生成された力の平均 FMH ( L ) は次式で記述される。

t  FMH (it L)  T  FMH ( L )it23ここで、 it L は it モンテカルロ・シミュレーションステップにおけるミオシンアームの伸びである。このとき、時間区分 [ T : T  T ] において、ミオシンヘッドが行う仕事量の増分は次式で与えられる。

WMH  FMH 
ここで、

SL0

2

SL0
は半サルコメアの長さ、  は筋線維方向ストレッチの増分である。さらに、半サ
2

ルコメアモデルを構成する単位体積あたりの仕事量の増分は次式で与えられる。
Wact 

SL ns nMH1
 0   FMH (im, is )  
Vmyo  ns 2 is im
SL ns nMH11
 0   FMH (im, is )  f  f : E  S act : E
Vmyo  ns 2 is im

ここで、 Vmyo は半サルコメアモデルの体積、 ns は半サルコメアのサンプル数、 FMH (im, is) は半サ
ルコメア is のミオシンヘッド im が生成する力である。f は自然状態における筋線維方向を表す
ベクトル、 S act , E はそれぞれ、第 2Piola-Kirchhoff アクティブ応力テンソル、Green-Lagrange ひずみ増分テンソルである。以上より、筋活動による収縮力 S act は次式で与えられる。
S act ここで、 C f 

Cf



ff

ns nMH1
  FMH (im, is) である。
Vmyo  ns is im

７）モンテカルロ・シミュレーションによる収縮力の計算例
短縮速度 u およびサルコメア長が与えられている場合には、単位断面積あたりの張力 C f を、有限要素解析を用いることなく計算することができる。Fig. 3.7 に短縮速度 u と単位断面積あたりの張力の関係を示す。ここでは、サルコメア長を 1.9 μm に固定し、一定の短縮速度 u を与え、張力のアンサンブル平均を求めた。ミオシンアームの弾性モデルには線形モデルと非線形モデルを用いて計算を行った。両モデルともに、短縮速度の増加により、筋収縮による張力は減少しており、これは筋線維を用いた実験による、筋線維短縮速度－力関係にも見られる関係と同様である。したがって、タンパク質の確率的な挙動からマクロスコピックな筋線維の短縮速度－力関係が再現できたといえる。
300
Linear model
Non-linear model

Active force [kPa]

250
200
150
100
50
0

0

1

2

3

4

5

6

7

Shortening velocity [m/s]

Fig. 3.7: Relation between shortening velocity and active force.

24８）まとめ筋骨格系のシミュレーションに関して、筋線維から筋肉全体までマルチスケール筋骨格系シミュ
レータの開発を引き続き行なった。骨格筋の収縮力発揮の最小単位であるアクチンフィラメントおよびミオシンフィラメントに対して状態遷移モデルを導入し、それぞれの状態遷移をモンテカルロ
法でシミュレートすることでタンパク質フィラメントの挙動を再現し、収縮力の算出する手法の導入を行った。今後、収縮力をマクロな骨格筋有限要素モデルに組み込み、タンパク質分子の確率的な挙動から骨格筋の収縮や関節運動を再現する予定である。
参考文献(1) Alves. J.L.. Yamamura. N.. Oda. T. and Teodosiu. C.. “Numerical simulation of musculo-skeletal
systems by V-Biomech”. In Proceedings of the 9th International Symposium on Computer Methods in
Biomechanics and Biomedical Engineering. (2010). Valencia. Spain.
(2) Yamamura. N.. Alves. J.L.. Oda. T.. Kinugasa. R. and Takagi. S.. “Significance of fascia for mechanicalbehavior of skeletal muscles -Numerical study using 3D models-”. In Proceedings of the 10th
International Symposium on Computer Methods in Biomechanics and Biomedical Engineering. (2012).
Berlin. Germany.
(3) 山村 直人.清水和弥.高木 周.“パーキンソン病振戦の 3 次元 FEM シミュレーション”.日本機
械学会 2013 年度年次大会講演論文集.（2013.9）.
(4) Johansson. T.. Meier. P. and Blickhan. R.. “A Finite-Element Model for the Mechanical Analysis ofSkeletal Muscles”. Journal of Theoretical Biology. 206 (2000). pp. 131-149.
(5) WASHIO, T., OKADA, J., TAKAHASHI, A., YONEDA, K., KADOOKA, Y., SUGIURA, S. ANDHISADA, T., "MULTISCALE HEART SIMULATION WITH COOPERATIVE STOCHASTICCROSS-BRIDGE DYNAMICS AND CELLULAR STRUCTURES", Multiscale Modeling and
Simulation, 11 (2013), pp. 965–999.
(6) 鷲尾 巧, "マルチスケール•マルチフィジックス心臓シミュレーションに関する研究−サルコメ
ア力学から心筋細胞構造を経て心拍動に至る解析手法の開発と応用−", Ph.D. Thesis,The
University of Tokyo,2013.
(7) Sussman, T. and Bathe, K.-J., “A FE Formulation for Nonlinear Incompressible Elastic and Inelastic
Analysis”, Computer & Structures, 26 (1987), pp. 357-409.

25

Ⅲ－１－２－５心疾患の合理的治療のためのマルチスケール・マルチフィジックス心臓シミュ

レーション（久田俊明・東京大学）
１）研究の背景と目的
生体は，根本的には分子間の相互作用を基本として成り立っているミクロのシステムである．例えば心臓を構成する筋肉は，次のような仕組みで収縮と弛緩を繰り返している．ペースメーカ細胞からの興奮伝播によって心筋細胞が電気的に興奮すると，細胞内の Ca2+濃度が高まり，筋線維内の 2 種のフィラ
メントが間に架橋構造（クロスブリッジ）が形成・解離を繰り返しながら互いにすべることで心筋が収縮する．逆に，心筋細胞の電気的な興奮が収まってくると，細胞内 Ca2+濃度は下がり，フィラメント間の架橋構造が外れ，心筋は弛緩し元の長さに戻る．この一連の筋収縮過程は一般に興奮収縮連関[1]と呼ばれ，モデル化のために生理学の分野では，クロスブリッジの平均濃度や収縮関連分子の平均伸び率を未知数とした時間発展型微分方程式が提案されている．ところが，そうした平均値を変数とした方程式
に対して，クロスブリッジ形成・解離に関わる確率的挙動，近傍の分子との協調性や個々の分子の伸びに依存した遷移率などの非線形性を矛盾なく導入することは非常に困難[2]である．筋収縮を分子のレベルで精緻にモデル化するためには，関連分子の動きを１つ１つモンテカルロ法で計算し，なおかつ，その結果を心筋の有限要素解析と力学的に矛盾なく安定的に結びつける計算手法が必要となる．

２）計算モデル
心筋細胞内では，”サルコメア”と呼ばれる微小な収縮機構が Z 線という構造を挟みながら長軸方向に連なった構造をしている．サルコメア内では，2 重らせん状の細いアクチン・フィラメントペアを 3
本の太いミオシン・フィラメントが取り囲むような周期的な配列が見られる（図 1）．ミオシン・フィラ
メント上にはアクチン・フィラメントと結合するための鍵縄状の構造（先端をミオシン・ヘッド，ひもをミオシン・アームと呼ぶ）が多数存在する．一方，アクチン・フィラメントペアには，ミオシン・ヘッ
ド結合サイトを覆い隠すトロポミオシンと呼ばれるさらに細いひもが巻きつき，周囲には Ca2+結合時にトロポミオシンを動かして結合サイトを顕わにするトロポニン複合体（T,I,C）が点在する．このトロポニン複合体の働きにより，結合サイトが顕わになっていると，ミオシン・ヘッドとのクロスブリッジが形成されやすくなる．また，ある結合サイトにミオシン・ヘッドが結合すると，トロポミオシンがさらに押し退けられて，近傍の結合サイトも露になる．この仕組みがクロスブリッジ形成・解離の協調性を生み出していると考えられている．
（図 2）
モデル化の際には，この構造に加え，結合したミオシン・ヘッドが首振り運動を行うことによってミ
オシン・アームが伸び，収縮力が発生すると仮定する首振り仮説を採用している．
我々のサルコメア・モデルは，ミオシン・フィラメント，アクチン・フィラメントペアで構成された
半サルコメアを対象としている．これは，サルコメアが Z 線を挟んで対称な構造をしているためである．
この半サルコメアにおいて首振り仮説に基づき，ミオシン・フィラメント上のミオシン・ヘッドには 6状態遷移モデルを，アクチン・フィラメント上のトロポニン C／トロポミオシンには 3 状態遷移モデルを採用した（図 1 右）．

26

図1

拍動解析のためのマルチスケールモデルの全体像

これらの状態遷移モデルには，次の 2 つの特徴がある．
１つ目は，ミオシン・ヘッドのクロスブリッジ形成・解離の遷移率が，近傍の分子の状態の影響を受けて変化する協調性（図 2）を導入した点である．具体的には，パラメータ
ヘッドの両隣で のクロスブリッジ数の遷移率はと着目するミオシン・を用いて，クロスブリッジ形成方向（図１赤矢印）に比例した値，クロスブリッジ解離方向（図１青矢印）の遷移率はに比例した値に設定している．これにより， が大きいほどクロスブリッジを形成しやすく，かつ，解離しにくくなる．
２つ目は，マクロ・スケールからの影響によってミオシン・ヘッドの首振りの遷移率が変化する点である．今，エネルギの観点からクロスブリッジ・サイクルを考えてみると， ミオシン・ヘッドは解離時に ATP を供給され，首を振る際には ATP の加水分解エネルギ をアームの歪みエネルギに変換して収縮力を発生させている．つまり，首振りはエネルギ変換過程であり，そこに統計力学的矛盾があってはならない．そのため，首振り確率をミオシン・ヘッドに蓄えられている内部エネルギ とアームの歪みエネルギ

から定められる Boltzmann 因子を基に決定している (図 3)．ここで，アーム 長 は，

ミオシン・ヘッド結合時の初期長，首振りによる伸び，結合後のフィラメント間すべり速度 を用いた次の式で決まる．
(1)
すべり速度 は，マクロ・スケールの連続体としての線維方向 に沿った短縮速度から決まるので，結果的にマクロ・スケールの動きが間接的にアーム長を変化させ，首振り確率をも変化させる．
これらの近傍分子との協調性およびすべり速度依存性を導入した上で個々の状態遷移をモンテカルロ法で計算すると，従来のモデルでは再現が困難であった現象を直接的にシミュレートできる.

27３）並列計算の方法と効果（性能）［ミクロ・メゾ間のマルチスケール解析］
ミクロ・メゾ間のマルチスケール解析においては，時間刻み幅（有限要素解析ののとモンテカルロ法）の大きなギャップと,サルコメア内のクロスブリッジ形成状態と数値細胞モデルの連続体としての剛性との適切な関連づけが問題となる．すなわち，モンテカルロ法においては，各時間ステップにおい図 2 近傍分子間の協調性て「遷移率×図 3 首振りの遷移率とエネルギの関係」が 1 を超えてはならないので， を小さくしなければならない．一方，有限要素解析においては陰解法を適用しているため，各時間ステッ
プにおいて線形方程式を解く必要があり，安定な解析が可能な範囲内で
つまり，ををできるだけ大きくしたい．

のように小さくすることは難しい．
，

例えば,実際の計算ではを採用しており, 500 倍もの開きがある．したがって，モン
テカルロ法の最終ステップにおける状態を有限要素解析で参照するだけでは，の間に起きたクロスブ

リッジ形成・解離の影響を連続体の剛性に反映できず,正確な解析にならない．そこで，時間刻み幅が大きく異なるモンテカルロ法と有限要素解析を適切に

図 4 有限要素法とモンテカルロ法の連成計算

カップリングする連成計算手法を考案した．（図 4）
はじめに，時間刻み幅

で状態遷移モデルに従って心筋要素内の各ミオシン・ヘッドの動きをモンテ

カルロ法で計算し，そのクロスブリッジ持続時間を記録しておく．途中，各時点でのアーム長 が必要になった場合は，を式(1)に代入して計算する．モンテカルロ計算が終了した後，連続体の時刻における変形を決める際には,クロスブリッジ持続時間とにおける変形速度を用いてアーム

長 を再計算し，そこから生じる収縮力を連続体に課す．すなわち， の増分からくる収縮力の変化とにおける変形の増分を結びつけて，剛性行列（連続体の変形増分と応力増分を結びつける行列）を計算している．

28

以上の方法で，連続体の微小変形とミオシン・アーム長の微小変化を結びつけ，任意の微小変形による仕事が連続体側とサルコメア側で一致するように連続体に作用させる応力テンソルいる．これにより，比較的大きな時間刻み幅を定めてでも安定的に計算を進めることができる．

［メゾ・マクロ間のマルチスケール解析］
メゾ・マクロ間のマルチスケール解析においては，基本的なアイデアとして均質化法[3]を採用し，メ
ゾ・スケールには局所的に周期的構造があると仮定した．ここで，メゾ・スケールの変位をロ・スケールの変位を， と をそれぞれのスケールにおける座標とすると，

ルの変形から来る均一な変形とメゾ・スケールでの周期的な変位成分

マクはマクロ・スケーの重ね合わせで定義される．
(2)
そして,仮想仕事の原理に基づく以下の力学的平衡方程式を解く.なお， はマクロ・スケールの参照配置
での構造領域を， はメゾ・スケールの周期的な参照領域をそれぞれ表すとする．
(3)
式(3)の左辺第 1 項はマクロの慣性力，第 2 項はマクロの仮想仕事（メゾの体積平均により定義）， はマ
クロ・スケールの境界に作用する力(心臓の場合は,内壁にかかる血圧など)である.

［2 段階マルチスケール解析］
すでに述べた通り，有限要素解析では陰解法を適用しており，全体の方程式は非線形であるためNewton-Raphson 法を用いて力学的平衡状態が全てのスケールに亘って成立するまで反復する必要がある.図 5 中央の巨大な行列は，その反復ステップで出現する線形方程式の係数行列（マクロの全体行列）の構成と非ゼロ構造を表している．このマクロの全体行列の(2,2)-部分行列（茶）はマクロの全体剛性行列であり，非対角部分行列（赤枠）はマクロとメゾの相互作用を表している．(1,1)-部分行列（青枠）はブロック対角行列であり，各ブロックがメゾの全体行列になっている．また，それらのメゾの全体行列も同様のブロック構造を持っている．
今，メゾ・スケールとマクロ・スケールの間で，解くべき線形方程式が下式のような形をしているとする．（図 5 も参照のこと）(4)
Schur 補行列を用いて，この係数行列を次のようにブロック LU 分解する．
(5)

この時，部分行列 はブロック対角行列になっているための計算はメゾ・ブロック毎に並列に実行可能である．2 段階マルチスケール解析では，各メゾ・ブロックとそこに含まれるミクロ・モデルを 1コアに割り当てて独立に求解させることで，大規模並列計算機においても高効率の計算が可能となっている．残りのマクロ解は，領域分割法で分割・並列化し求解する．この部分の並列度はあまり高くないので，マクロ求解の間，それ以外のメゾ・ミクロ担当プロセスが待ち状態になってしまうが，ミクロ求解，メゾ求解やモンテカルロ部の計算負荷が大きい場合は特に問題とならない．

29

図 5 2 段階マルチスケール解析手法の概要
４）研究成果我々は，サルコメア分子の動きから心臓の拍動までを結びつけたシミュレーションを活用することで，肥大型心筋症という心疾患のメカニズムを解明することを目指した．肥大型心筋症とは，心筋壁が異常に厚くなった結果，左心室の容積減少や拡張障害など，心機能の低下が生じる心疾患である．罹患すると心臓に過負荷が掛かるため，健全に見えたスポーツ選手や若年者が突然死してしまうことが問題となっている．しかし，すでに遺伝子解析によって原因遺伝子が同定されているにも関わらず，未だ病態
形成のメカニズム解明には至っていない．
その肥大型心筋症患者の心臓では，心筋細胞や筋原線維の中で線維方向が乱れる特徴的な構造がよく見られる．そこで，「京」でのシミュレーションでは，このミクロ・レベルの錯綜配列が，45,256 個のメゾ要素を持つマクロ・モデルにおいて拍動性能にどのような影響を与えるのかをシミュレーションで観察してみた．
そのために，メゾ・モデル内（図 5）の線維方向を，数値細胞モデルの長軸方向 に対して角度 の範囲内で適当に乱した方向とした簡易的な錯綜配列モデル（図 6）を考え，拍動シミュレーションを行った．

図 6 簡易的な錯綜配列モデル30

左心室の圧力と容積の履歴を図 7 に示す．この履歴には，弛緩速度と拍出量の減少という肥大型心筋症に似た心機能低下が見られ，乱れ（角度 ）が大きいほど顕著である．前者は，線維方向の乱れから収縮方向が揃わないことですべり速度 が低下し，収縮時のアーム長が伸びて首振り遷移率が小さくなった結果（図 3 の原理）として，クロスブリッジ・サイクルが遅くなったものと考えられる．後者は，乱れによって収縮力の生じる方向が揃わなくなったために，心筋線維方向の収縮が弱まり，心臓全体の収縮能が落ちたと考えられる．
収縮末期の収縮力分布（長軸断面）を図 8 に示す．この図においても，流体部の径から収縮能の低下が確認でき，θ <20°のときは速やかに力が抜けているのに対し，θ <60°のときは外壁から中層にかけて力が残っていることも確認できる．

図 7 左室圧と左室容積の履歴

図 8 収縮力分布（長軸断面）
５）まとめと今後の課題
医学の発展に伴い，関連遺伝子の網羅的スクリーニングによって，原因となる分子レベルの異常が判明したにも関わらず，病態のメカニズムが解明されていない心疾患が残されている．我々はその解明に取り組むため，臓器・細胞・サルコメアまでをモデル化したマルチスケール解析手法を考案し，「京」でその効果を確認した．
今後はそのミクロ・スケールからマクロ・スケールまでを矛盾なく結びつける特徴を活かし，サルコ
メア構成タンパク質の異常から出発して，心筋リモデリングなどの過程を経て，肥大型心筋症などのマクロな病態が形成されていく様子を再現することを目指したい．そのためには，中～大規模の長時間計算を繰り返さねばならず，「京」のような計算環境が必要になると考えている．
31

参考文献
[1] S. Sugiura, T. Washio, A. Hatano, J. Okada, H. Watanabe, and T. Hisada, Multi-scale simulations of cardiac
electrophysiology and mechanics using the University of Tokyo heart simulator, Prog Biophys. Mol. Biol., 110(2012), pp. 380-389.
[2] T. Washio, J. Okada, S. Sugiura, and T. Hisada, Approximation for cooperative interactions of a
spatially-detailed cardiac sarcomere model, Cell. Mol. Bioeng., 5 (2012), pp. 113–125.
[3] T. Washio, J. Okada, A.

Takahashi, K. Yoneda, Y. Kadooka, S. Sugiura, and T. Hisada, Multiscale Heart

Simulation with Cooperative Stochastic Cross-Bridge Dynamics and Cellular Structures, T.Multiscale Modeling &Simulation 2013 11:4, pp. 965-999.

32

Ⅲ－２

後藤

信哉（東海大学）

抗血小板薬薬効予測を可能とする階層統合シミュレータの開発に向けたモデリングⅢ－２－１ 実施計画
本研究では、「予測医療に向けた階層統合シミュレーション」研究の一環として、「予測医療の実現に向けた抗血小板薬薬効予測を可能とする階層統合シミュレータの開発に向けたモデリング」の研究開発を実施する。
また、「予測医療の実現に向けた抗血小板薬薬効予測を可能とする階層統合シミュレーターの開発と大規模計算」の研究を行う上で、関連する研究者と必要な協議等を行うとともに、本格実施に必要な研究体制の整備を行う。
平成２５年度は以下の研究開発を実施する、（１）東京大学・高木グループの「心筋梗塞・脳梗塞のマルチスケールシミュレーション」モ
デルに対して、実証実験に裏打ちされた血小板細胞の活性化と血液凝固系の連携モデルを作成する。
（２）実証実験に裏打ちされた血小板細胞代謝、活性化モデルの精緻化を継続する。
（３）血小板接着における von Willebrand 因子と血小板膜糖蛋白 GPIbαの分子スケールにおける接着モデルと血小板細胞接着モデルの連成を継続する。
（４）活性化血小板表面における凝固系の活性化を取り込んだ血液凝固、フィブリン産生モデ
ルを作成する。抗血小板薬による血小板活性化阻害とフィブリン血栓形成阻害を連成させるモデルを作成する。

Ⅲ－２－２

実施内容（成果）（１） 東京大学・高木グループの「心筋梗塞・脳梗塞のマルチスケールシミュレーション」モデ
ルに組み込むための血小板細胞の活性化と血液凝固系の連携モデルを作成した。血小板細胞の活性化と、活性化血小板細胞周囲における血液凝固系を連携させるための概念的モデルを英文総説論文として発表した（図1:Goto S, et al. Clin Cardiol.
2014）。従来、血小板細胞の活性化と凝固系の活性化は別個の事象として研究されて来た。
本モデルでは活性化に伴う血小板細胞膜の脂質構成の変化に注目した。陰性荷電したリン脂質の細胞表面への発現に応じて、凝固因子が血小板細胞膜上に集積して効率的に凝固系を活性化させるモデルとした。血管壁におけるトロン
ボモデュリンを介するトロンビンの不活化モデルも入れ込
33

んだ。溶解した分子同士の化学反応である血液凝固カスケードと、細胞反応である血小板細胞活性化を連成した概念的モデルを作成した。
本モデルは従来の臨床的経験、心筋梗塞発症予防における各種薬剤による予防効果、細胞スケー
ルの実証的実験成果などを取り込んだ定性的モデルである。今後、高木グループの「心筋梗塞・脳梗塞のマルチスケールシミュレーション」への取込みを目指して、本モデルに取り込んだ各ス
テップの数理モデル化を目指す。
（２）

実証実験の結果を用いて、血小板細胞代謝、活性化の基盤モデルを作成した。汎用性に乏しかった従来のシミュレーターを Systems Biology Markup Language(SBML)に準拠すべく修正した。本モデルでは空間の拡張よりもモデル内にて扱える物質の拡張を目指した。代謝、活性化を担う細胞内の化学反応では、反応定数が未知の場合が多い。血小板細胞では、細胞内の代謝、活性化反応と同時に、血小板細胞表面における凝固反応を組み込むことが可能である。
図 2 に血小板細胞内、膜上、細胞外の化学反応を一括して記載したモデルを示す。
本モデルは将来の論文発表を目指しているため、
本報告書内では反応に寄与する
物質名を読み取れないようにしてある。細胞内／膜上／細胞外を統合したモデルとした。細胞外の反応は血液凝固にかかわる化学反応を取り入れた。この細胞外の化学反応については反応係数が既知である。細胞膜、細胞内の反応については反応係数が未知である。証実験とシミュレーション計算の連携により血小板細胞反応の統合的理解を目指す。
（３）

平成 25 年度には、血小板接着における von Willebrand 因子と血小板膜糖蛋白 GPIbαの分子スケールにおける接着モデルと血小板細胞接着モデルの連成に力を注いだ。血小板細胞膜上の GPIbα、血管壁の von Willebrand 因子ともに蛋白質としての構造は決定されている。
両者の接着には、GPIbαの N 末端部分、von Willebrand 因子の A1 ドメインが寄与することは理解されている。両者が接着したときの接着部位の構造も結晶構造としては確認されている（図 3 右）。本プロジェクトの予算を主体に東海大学にて雇用されている塩崎は、vonWillebrand 因子と GPIbαの接着時の構造から両分子間の重心の距離に対する potential ofmean force を molecular dynamics（MD）の手法を用いて計算した。計算結果を図 3 左の赤点にて示す。計算結果から Morse potential にフィッティングした結果を赤線に示す。この34

potential of mean force の微分値から、両分子間に作用する平均的な力（pN）を計算した（青線）。本計算結果から、von Willebrand 因子と GPIbαの間に平均的に働く力の最大値は 60 pN程度と予測した（図 3 青線の最大値）。
本計算では両分子の重心間距離とpotential of mean
force の関係を計算の対象とした。実際、血小板細胞が GPIbαとvon Willebrand 因子
を介して血管壁に接着するときの両分子の接着により惹起される接着力と分子動力学計算により求めた両分子間の potential of mean force の関係は未知である。両者の関係を示す論理的モデルが必須となる。
塩崎が計算により求めた von Willebrand 因子と GPIbαの接着力は、Kim らが原子間力顕微鏡を用いて求めた値に近い（Kim J, et al. Nature. 2010）。また、本研究協力者である須藤らの原子間力顕微鏡を用いて計測した結果も数十 pN というオーダーとしては近い。しかし、これら原子間力顕微鏡により計測した von Willebrand 因子と GPIbαの接着力は、GPIbαを塗布したカンチレバーが von Willebrand 因子面から剥がれた時に両者に作用していた力であり
potential of mean force から計算した力とは質が異なる。分子動力学計算により算出した
potential of mean force から算出される力と、両分子の接着が剥がれる時の力の関係の連携に確率を媒体とする試みを行い、今後も継続する。
分子動力学計算から算出可能な potential mean
force を、細胞スケールの血小板接着力と連携させるための定量的実証実
験を多数施行した。シミュレーション計算の擦り合わせのためには定量性が必須である。血液をvon Willebrand 因子上に灌流し、血小板細胞を接着させたのち、血小板細胞、von Willebrand 因子を含まない pH 7.4 の等張のバッファーを灌流した。
時間経過とともに接着し35

た血小板が剥離する。単位時間当たりの剥離数はバッファーを灌流する血流条件に依存した（図 4）。von Willebrand 因子上への血小板の接着の初期段階に GPIbαが必須であることは、
如何なる実証実験においても常に示されている。しかし、接着力を規定する因子には時間、vonWillebrand 因子など複数の因子が寄与する。全ての寄与因子を明確にするには至っていない。
今後も条件を変えた実証実験を継続する予定である。
（４）

血小板細胞表面におけるトロンビン受容体 PAR-1 は、凝固系活性化の結果として産生されたトロンビンにより活性化される。PAR-1 受容体を阻害することにより、血小板活性化と凝固系活性化の positive feedback を阻害することにより効率的に心筋梗塞の発症を予防できると想定された。実証実験にて PAR-1 阻害により血小板の接着が部分的に阻害されることを示した。PAR-1 阻害による血小板接着阻害効果は、ADP 受容体 P2Y12 受容体阻害時ほど顕著ではなかった。実証実験結果の論文化を進めている。
東京大学・高木グループの「心筋梗塞・脳梗塞のマルチスケールシミュレーション」では、血小板細胞の活性化が寄与しない接着早期を対象としている。血小板細胞の活性化が寄与する
実証実験の成果を現時点のシミュレーションモデルに組み込むことは困難である。抗血小板薬
による血小板活性化阻害とフィブリン血栓形成阻害を連成させるモデルの妥当性を示唆する臨床試験の結果を論文として発表した（Morrow DA, et al. Stroke 2013）
。

36

Ⅲ－３

中村 仁彦（東京大学）

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

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

Ⅲ－３－２

実施内容（成果）（１）モデル統合インターフェイス「脳神経系の多階層モデル」、「脊髄と四肢筋骨格モデル」、「筋繊維レベルからの骨格筋モデル」
を統合するインターフェイスの開発を行ってきた。東京大学の高木チーム、沖縄科学技術研究基
盤整備機構の銅谷グループと議論を行いながら、我々が開発してきた全身筋骨格シミュレータ、
高木チームの詳細筋有限要素シミュレータ、銅谷グループの脳中枢神経系シミュレータの統合イ
ンターフェイスを設計してきた。統合インターフェイスは、ニューロシミュレータ NEST[R1]、骨格筋有限要素解析シミュレータ V-Biomech[R2]、全身筋骨格シミュレータ sDIMS [R3]、大規模なプロセス間の通信を可能にする MUSIC[R4]から構成される。
前年度から開発してきたインターフェイスを基盤として、神経系・筋骨格系の統合シミュレー
タを開発した(図１)。本シミュレータでは、(a)神経の信号伝達の考慮、(b) 任意形状の筋への対応、(c) 変形や応力分布の詳細解析、(d) 全身統合計算を可能とする。項目(a)の実現のため
膜電位変化を考慮したスパイク発火型の神経モデルを利用し、項目(b)および(c) を実現するた
めに筋の解析には有限要素法を用い、有限要素メッシュとして筋をモデル化している。骨格は剛体を仮定し、MPI による並列化により項目(d) を実現する。統合シミュレータは神経系シミュレー
タと筋骨格系の力学シミュレータの二つに大きく分けられる。後者は前者から受け取った信号をもとに筋活動を算出し、これに基づいて物体間の接触・接続を考慮した筋変形と骨格運動を求める。求めた運動に応じた感覚信号が算出され，神経系シミュレータへ送られる。
統合シミュレータは前述の通り複数のライブラリ・シミュレータが連動して起動する。これら
全てのソフトウェア環境が「京」で利用可能となり、また連動して動作することも確認できた。

37

図１．神経筋骨格シミュレータの概要（２）全身筋骨格モデルの動力学計算アルゴリズムの並列計算
ａ）．複数物体間の接触力の並列計算アルゴリズムの開発
筋肉の体積効果を考慮した全身筋骨格の運動計算において、各筋・骨同士の接触を考慮するため、複数物体間の接触力の並列計算アルゴリズムを開発した。接触力計算自体は有限要素メッ
シュに付与された深さ情報に基づくペナルティ法[R5]に基づく。また接触力導出の際に必要となる干渉検出には BVH(Bounding Volume Hierarchy)[R6] を利用する。図 2 に示すように, 各物体の BVH の上位層を統合した物体全体の BVH を持つプロセスを用意することで、MPI 並列実行を行う。全体の BVH を利用し干渉候補を絞り込むことで、データ通信量を抑えて干渉検出を並列実行できる（図２）。
マルチプロセス干渉計算アルゴリズムを評価するため、複数の有限要素モデルを等間隔で直列に配置して干渉計算に要した時間を測定した(京コンピュータ)。有限要素モデルのサイズは節点数 63000、要素数 59126 である。一つのプロセスにつき一つの有限要素モデルの形状データを担当するものとし、並列実行するプロセス数と同数の有限要素モデルの接触力計算を 1000 回繰り返し実行時間を計測した。プロセス数は 8, 16, 32, 64, 128, 256, 512, 1024, 2048 の 9 種類とした。計測結果を図 3 に示す。1024 プロセスでの実行時間は 8 プロセスでの実行時間の 1.08
倍であり、1024 プロセスまでのスケーラビリティが確認できた。2048 プロセスでの実行時間は
8 プロセスでの実行時間の 1.25 倍となったが、有限要素メッシュ数が少ないため、マスタープ
ロセスの負荷が相対的に増加したためである。ただし最終的な全身シミュレーションにおける筋数は 1000 程度を想定しており、全身統合シミュレーションにおける接触力計算への配分に対して、現状のスケーラビリティで対応可能である。

図２．分散メモリ型接触力並列計算の概要
38

図３．接触力計算の並列性の評価ｂ）．弾性体・剛体の連成運動の並列解析アルゴリズムの開発
筋肉の体積効果を考慮した全身筋骨格の大規模動力学シミュレーションに向けて、弾性体・剛体の連成運動の並列解析アルゴリズムを開発した。骨は剛体、関節は機械ジョイント，筋・腱は超弾性と仮定することで、計算量およびメモリ量を大幅に減少させる。各運動方程式は非線形となり Newton-Raphson 法により求めるが、その過程で必要となる大規模疎行列を係数に持つ線型方程式を並列計算する。筋有限要素ごとに並列計算可能であり、またデータ転送量はその筋に接続されるリンクの個数に依存するため、転送コストを大きく抑えることができる。
弾性体・剛体の連成運動の並列解析アルゴリズムを評価するため、剛体・弾性体を交互に接続
した振り子モデルを利用して計算時間を測定した(Intel Xeon CPU E5-2687W(8 コア) x 2)。一つのプロセスにつき一つの有限要素計算を割り当てた。プロセス数は 4,8,12 として計算時間は表１のようになった。有限要素メッシュ数が少ないため、マスタープロセスの負荷が増加してしまう結果となったが、プロセス数の増加とともに計算時間を削減することが確認できた。
表１．弾性体・剛体の連成計算の並列性の評価c）．弾性体・剛体の連成運動の並列解析アルゴリズムの開発
計算論的神経モデルから人間の身体の物理モデルへの結び付けるための人間の解剖学的幾何
モデルを用いて、詳細な筋骨格モデルの作成を前年度から継続して行ってきた。一方で、解剖学
39幾何データの入手困難な部位を考慮して、ワイヤ型の全身筋骨格モデル[R3]から，有限要素筋モデルを自動生成するアプローチも並行して行った。自動生成の流れは以下の通りである。まず補間によりワイヤから曲線を作成し、立体形状を埋め尽くすよう球を配置する。次に球の位置、間隔などを維持しながら、異なるワイヤに対応する球同士と骨を表すポリゴンとの干渉を避けるよ
う最適化計算を行う。最後に修正された球の位置に基づいて、有限要素メッシュを作成する。ワ
イヤ型全身筋骨格モデルから筋有限要素全身モデルへの変換の様子を図４に示す。自動生成された全身筋有限要素モデルを基盤として、解剖学に基づいた正確な筋形状や筋繊維方向が分かる部
分については、それを反映したモデルへと順次更新していく予定である。

図４．全身ワイヤ筋骨格モデルから全身筋有限要素モデルへの自動変換（３）脊髄の筋支配モデルの開発
前年度では、脊髄神経モデルとワイヤ型筋モデルを接続するインターフェイスを開発した。本
年度ではそれを拡張し、筋有限要素メッシュとの接続インターフェイスを開発した(図 5)。感覚神経は筋紡錘に相当する部位にのみ接続すると仮定し、運動神経は筋収縮力を発揮する領域全てと接続すると仮定する。
文献[R7]、[R8]を基盤として、運動神経信号から筋活動への変換、筋の状態から感覚神経信号への変換を行う。運動神経の信号は、筋線維方向の応力成分に変換する。求めた応力成分は、筋有限要素の応力成分として弾性体解析に組み込まれる。また、感覚神経の信号発生時刻は筋線維方向歪速度から導出する。
筋弾性体と脊髄神経モデルのインターフェイスの評価を行った。上腕二頭筋に 120 の運動神経と 51 の感覚神経を接続し、上腕三頭筋に 240 の運動神経と 103 の感覚神経をランダムに接続し、全運動神経にポアソン分布（平均発火頻度 100Hz)に従って入力し、等尺性収縮のシミュレーショ
ンを行った。シミュレーション中の筋収縮応力、筋歪速度、信号入出力の一例として上腕二頭筋
短頭部のものを図 5 に示す。骨や筋同士で大きく干渉することなく筋収縮が起きており、等尺性40

収縮のシミュレーションが実現できている。筋歪速度の時系列データから分かるように、筋は初めの 0.4 秒で大きく変形し、その後は減衰して平衡状態に近づいていく。感覚神経の発火信号を見ると初めの 0.4 秒で発火頻度が大きく変化しており、運動情報から神経信号情報への変換が行われていることが確認できた．

図５. 筋弾性体と脊髄神経ニューロン間の信号例。 (左上) 運動神経に接続する各有限要素領域の応力, (右上)筋紡錘に相当する有限要素領域の歪速度，(左下)運動神経のスパイク信号列，(右下)筋紡錘のスパイク信号列。

参考文献[R1] M-O. Gewaltig, M. Diesmann. NEST (Neural Simulation Tool), Scholarpedia, 2(4): 1430,
2007.
[R2] 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.
[R3] 山村直人, ALVES Luis, 小田俊明, TEODOSIU Cristian, ヒト骨格筋力学解析のための有限要
素シミュレータの開発. JSME annual meeting: 29-30, 2009.
[R4] 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.
[R5] G. Hirota, S. Fisher and C. Lee. An implicit finite element method for elastic contact.
Proceedings of The Fourteenth Conference on Computer Animation, pp.136-254, 2001.
[R6] G.V.D. Bergen. Efficient collision detection of complex deformable models using AABB
trees. Journal of Graphics Tools, Vol. 4, No. 2, pp. 1–13, 1997.
[R7] R.R,L. Cisi, A.F. Kohn. Simulation system of spinal cord motor nuclei and associatednerves and muscles, in a Web-based architecture. Journal of Computational Neuroscience,25(3): 520-542, 2008.
[R8] A.Prochazka, M.Gorassini.Models of ensemble firing of muscle spindle afferents recorded
during normal locomotion in cats. Journal of Physiology, 507(1), 277-291, 1998.
41

Ⅲ－４

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

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

数値モデルを構築し、それらをつなぐ大規模シミュレーションにより、心筋梗塞、パーキンソン病などの発症プロセスを解明し、治療手法の探索評価を可能にすることをめざしている。その一環として、脳神経系の感覚入力から運動出力までをつなぎ、脳全体の回路から細胞、分子レベルのダイ
ナミクスを含む多階層の脳神経系モデルの構築が不可欠である。
平成２５年度は、大脳基底核を中心とした運動制御系の神経回路モデルを脳神経系シミュレータ
NEST により実装し、「京」上で大規模並列計算を実施する。これにより健常時の運動制御機能の実
現と、パーキンソン病における振戦などの病態の再現に向けて、パラメタ探索とモデル改良を進める。
また、「行動制御の脳神経系の多階層モデルの構築」の研究を行う上で、関連する研究者と必要な協議等を行う。
Ⅲ－４－２

実施内容（成果）（１）大脳基底核モデルの構築
昨年度に自然科学研究機構生理学研究所の橘吉寿助教（当時）に提供していただいたパーキンソン病モデルのサルの神経細胞の活動データ（文献 1）を、ドーパミン枯渇前後の神経活動の特徴の変化について解析した。その結果、ドーパミン枯渇によって、視床下核、淡蒼球内節および外節の神経細胞の一部がパーキンソン病に特徴的な 8-15Hz のベータ帯域振動活動を間欠的に示すことを確認した。そして、8-15Hz のパワースペクトル密度が高い時間帯では、単一細胞がリズミックな
バースト活動を示していることを確認した。
さらに近年の生理学および解剖学の知見を踏まえて、この単一細胞レベルでの 8-15Hz のリズ
ミックなバースト活動を再現するための視床下核−淡蒼球外節の局所神経回路網モデルを提案した。
計算機シミュレーションによって、この神経回路網モデルでは、神経細胞間の結合強度、および淡蒼球外節の自発発火頻度を変更することによりドーパミン枯渇前後（正常時 vs パーキンソン病時）に特徴的な神経活動を再現できることを示した（図１）。
（２）視床–大脳皮質モデルの開発
運動皮質と視床は、大脳基底核とともに、運動に関する情報処理を行っていて、パーキンソン病の発生機構にも関与しているとみられる。そのため、パーキンソン病の病態を再現する全身筋骨格―神経統合シミュレーションに運動皮質と皮質のモデルを導入する必要があり、大脳皮質と視床についての情報処理機構、パーキンソン病における神経活動について、文献調査とモデルの開発を行った。
パーキンソン病の代表的な症状の一つに振戦があるが、その物理的な振動は運動皮質で見られる
10Hz 程度のアルファ周波数帯の脳波と 1:2 の引き込みを示す(文献 3)。アルファ振動の起源は、主に視床と考えられていて、その発生機構も電気生理学的実験でよく調べられ、モデルも提案されて42

いる。Vigayan らの視床のモデル(文献 4)を NEST シミュレータで実装し、シータ～アルファ周波数帯の振動を再現した(図２左)。

図１：大脳基底核回路網モデルでの視床下核神経細胞の活動

運動皮質では、異なる体部位の運動が運動皮質の異なる領域、コラム、クラスターごとに表現されていて、機能的に組織化された神経回路によって実行されていることが分かっている。運動機能とパーキンソン病の病態のモデルによる再現には、異なる回路による異なる motor primitive の実行を再現する必要があるため、運動皮質の結合に関するデータを取り入れ、複数の motor primitiveを表現可能なモデルの構築を目指した。特に、レーザー刺激によって各種神経細胞間の結合の空間的な広がりを調べる実験データを取り入れ、運動皮質の興奮性―抑制性の各種神経細胞間の空間的
結合パタンをモデルに導入した（図２右）。モデルの 2/3 層で、異なる運動を表現すると想定される隣接する 2 つの部位に、同時に信号を与えたとき、空間的な結合の広がりによる側抑制がおき、強い信号を受ける領域の活動だけが起き、弱い信号を受ける領域は抑制される、lateralcompetition が見られた。この現象は、運動皮質の motor primitive の選択の情報処理に働く可能性がある。また、拮抗筋の相互抑制回路として、パーキンソン病の振戦の原因となる可能性もある。
この運動皮質のモデルの規模を拡張し、「京」の 500 ノードを用いて実行試験を行った。実行されたモデルは、4.8mmx 4.8mm の皮質表面に相当し、マウスの運動皮質全体の規模に相当するものである。

43

図２：視床細胞モデルのアルファ周波数帯の Ca2+バースト（左）と3 次元空間で描かれる運動皮質モデルの結合と層構造（右）（３）京コンピュータ上での NEST による 17 億個の神経細胞からなる神経回路シミュレーション
京コンピュータ上で NEST シミュレータによって実行可能な神経回路モデルの規模と動作確認を行った。実行試験は、NEST の開発を主に行っている Juelich Research Centre の Markus Diesmann
氏、Abigail Morrison 氏らのグループ、特に同グループの Moritz Helias 氏、Susanne Kunkel 氏らとの協力により行われた。今回用いた NEST シミュレータは、メモリの消費量をできる限り抑えて実行規模の拡大を可能にしながらも、多様な神経細胞やシナプスの振る舞い、シナプス可塑性などが再現可能になっている。
実行したモデルは、Morrison 氏らが発表した大脳皮質局所回路モデル（文献２）で、スーパーコ
ンピュータによる NEST シミュレータの性能測定で用いられているものである。試験の結果、汎用的な神経回路モデルとしては世界最大規模の 17 億 30000 万個の神経細胞と、10 兆 4000 億個のシナプス結合からなる神経回路シミュレーションに成功した(図３,文献３)。実現した神経回路規模は、
小型のサルの全脳の神経細胞数に匹敵する。パーキンソン病の病態再現モデルを実行する段階では、

図３：「京」で実行した大脳皮質局所神経回路モデルの計算規模と計算時間
44

筋骨格モデルを同時実行する計算リソースの分、神経回路モデルの規模は小さくなるとみられるが、計算リソースが半減したとしても、パーキンソン病の病態を再現する最小限の回路規模を実行するには十分可能であるとことがわかった。
（４）NEST と MUSIC によるモデル統合
これまでに開発した大脳基底核モデル、視床–大脳皮質モデルを、東大グループにより開発された脊髄神経回路−筋骨格系モデルと接続するため、MUSIC (MUltiple SImulation Coordinator) を K上に移植し、並列実行実験を行った。大脳皮質 5A 層のニューロンのスパイク出力が大脳基底核モ
デルの線条体ニューロンに接続され、大脳基底核モデルの淡蒼球内節(GPi)ニューロンの出力が視床ニューロンに接続され、大脳皮質-基底核-視床ループを構成する。また大脳皮質 5B 層の錐体路(PT)ニューロンの出力が脊髄神経回路に送られ、関節の拮抗筋の制御を行う（図４）。神経回路モ
デルは NEST の異なるバージョン、筋骨格系モデルは DIMS による実装であるが、正しく相互通信のもと並列実行が行えることを確認した。
今後、モデル間の結合および各モデルのパラメタを調整し、より現実的な運動制御とパーキンソン病態の再現をめざす。

図４：大脳基底核–視床–大脳皮質–脊髄–筋骨格系の統合モデル

45

参考文献
1) 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
2) Morrison A, Aertsen A, Diesmann M. (2007) Spike-timing-dependent plasticity in balanced
random networks.

Neural Comput., 19(6):1437-67.

3) 理研、OIST、Juelich Research Centre のプレスリリース 「京（けい）」を使い 10 兆個の結合の神経回路のシミュレーションに成功－世界最大の脳神経シミュレーション－(http://www.riken.jp/pr/topics/2013/20130802_2/)

46

Ⅲ－５

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

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

Ⅲ－５－１

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

１） 成果①（発表論文１）
ヒト静止立位姿勢制御に関して、関節の受動的粘弾性が低くかつ能動的神経フィードバッ
ク制御のゲインも小さい（関節がコンプライアントである）という生理学的に妥当な条件の下で、静止立位姿勢を柔軟に安定化する間欠制御仮説の検証を行った。具体的には、我々の研究グループで開発したヒューマン・コンピュータインタフェイスを用いて、拮抗筋の収縮力によって制御される計算機内仮想倒立振子を被験者自身の 2 つの足関節筋（前頸骨筋および腓腹筋）の筋活動によってバランスする心理物理課題（運動学習課題）を実施した。その結果、実験に参加した多くの被験者は、2 つの筋（あるいはどちらか一方の筋）の活動を振子の状態に依存した適切なタイミングで間欠的に休止させることによって、倒立振子の立位
姿勢を達成させていることを明らかにした。さらに、被験者に制御される倒立振子の運動特47

性を様々な評価指標に基づいて定量解析した結果、これらの被験者が採用した倒立振子のバランス戦略は、我々が提唱する姿勢の間欠制御モデルと一致する特性を有していることが示された。この結果は、東大高木・中村グループ、OIST 銅谷グループが構築する大規模並列モデルが実現すべき神経系による身体運動制御のマクロレベルの特性として、制御の間欠性が重要な役割を果たしていることを示唆している。
２） 成果③（発表論文２）
本研究課題の研究協力者である国立病院機構刀根山病院・佐古田院長とともに、すくみ足
症状の無いパーキンソン病患者およびすくみ足症状はあるがパーキンソン病ではない患者の歩行運動計測データの解析と計測された運動を再現する単純な数理モデル化とその数値
シミュレーションを実施した。その結果、すくみ足症状を呈する患者の歩行では、平均的な振舞いとしては左右脚運動の相対位相を１ステップで（ある意味強引に）平均値に修正する
タイプの左右脚協調が見られる反面、協調制御の神経指令はステップ毎に大きく揺らいでいる（位相修正量の確率的変動が大きい）ことが明らかになった。この結果は、パーキンソン病における肢間協調の崩壊において、運動指令の確率的なゆらぎの存在が重要な役割を果たしていることを示唆している。
これら 2 つの成果は、マクロレベルにおけるヒト立位姿勢制御の特性と制御様式を理論的・定量的に明らかにするものであり、本プロジェクトで開発中のマルチスケールなヒト身体運動制御系が再現すべき現象とその背後にある制御メカニズムを提供するものと考えられる。
発表論文リスト
1. Y Asai, S Tateyama, T Nomura. Learning an intermittent control strategy for postural balancing
using an EMG-based human-computer interface, PloS ONE 8 (5), e62956
2. T Tanahashi, T Yamamoto, T Endo, H Fujimura, M Yokoe, H Mochizuki, T Nomura, S Sakoda.
Noisy interlimb coordination can be a main cause of freezing of gait in patients with little to no
Parkinsonism, PloS ONE 8 (12), e84423

48