Ⅲ

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

Ⅲ－１

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

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

Ⅲ－１－１

ミュレーションに向けた血栓シミュレータと冠循環モデルの連成及び筋繊維の集合体としての骨格筋のシミュレーション」の研究開発を実施し、研究開発を統括する。
また、プロジェクトを推進する上で、関連する研究者と必要な協議等を行い、予測医療に向けた階層統合シミュレーション」の研究を統括する。
平成２４年度は、前年度に引き続き，心筋梗塞シミュレーションに向けた血栓シミュレータと冠循環モデルの連成として、現在進行中のグランドチャレンジプログラム「次世代生命体統合シミュレーションソフトウェアの開発」で開発中の血栓シミュレータと、心臓を取り巻く血管である冠動脈のシミュレーションを連成させて、フィブリン網の発達により血栓が凝固して血管を閉塞させるプロセスを再現する手法の開発を進める。本年度は特に、動脈硬化巣に吸着して活性化した血小板により引き起こされる生化学反応を再現し、抗血小板薬として知られるクロピドグ
レルの薬効を再現するのを目指す。また、抗血小板薬として成功しなかった事例として GPIIb,
IIIa の阻害薬の薬効評価のシミュレーションを実施する。
筋繊維の集合体としての骨格筋のシミュレーションの研究課題に関しては、前年度に引き続き、
ヒトのふくらはぎ部分にあたる下腿三頭筋を対象とする。本年度は、運動ニューロンからのシグナルにより筋繊維が収縮し、その集合体として骨格筋全体が収縮する計算モデルを構築し、実際に「京」上で並列計算を始める。

Ⅲ－１－２

実施内容（成果）

本研究は、（Ａ）心筋梗塞のシミュレーションに向けた血栓症のモデリングの部分と、（Ｂ）神経疾患による運動機能破綻のシミュレーションの一部として、推進している筋繊維の集合体としての骨格筋収縮のマルチスケールシミュレーションの部分からなる。以下、それぞれに分けて説明を行う。
（Ａ）「心筋梗塞のシミュレーションに向けた血栓症のモデリング」（１）ソフトウェアの開発・高度化の状況
１）進め方
グランドチャレンジプログラムで開発されたオイラー型定式化による流体・弾性体連成解析
手法(ZZ-EFSI)による赤血球と血小板を含む血流の計算を基盤として、血小板の一次凝集から二次凝集、その後の凝固へと発展する血栓が、血管閉塞に至るまでのプロセスについて、様々の分子種の化学反応を考慮して連成させながら計算する手法の構築を目指している。
血栓形成の初期段階で生じる血小板凝集は、血小板が血管壁表面に粘着していく一次凝集と、血小板が活性化して血小板同士の粘着から発展する二次凝集の二つの段階に分けられる。一次凝集では、血小板表面にある粘着系蛋白質と血管壁に吸着している蛋白質との間の結合が重要な役割を果たすが、血小板と血管壁間に働く力の評価には、動的モンテカルロ法から得られた
血小板表面での粘着系蛋白質の局在や拡散の効果をモデル化して取り入れる。蛋白質間の結合には、赤血球や血小板などの力学的な相互作用が血栓の形成に大きく関与しているものと考えられため、血流との連成計算が必要となる。二次凝集においては、血小板の活性化により生じる内部変化や凝集を促進させる物質の放出等による代謝作用が重要になる。個々の血小板内での化学反応によって生じる物質量の変化と血小板表面から放出される物質量については、反応
速度論に基づいた常微分方程式の連立式によって計算する。これと血小板周りの流れ場の計算と連動させることで、血小板が凝集して生じる血小板血栓の形成過程についての計算を行っていく。
２）状況
血小板凝集が進展していく段階では、血小板自体が活性化を促進する物質を放出し、血小板表面にある受容体と結合することによって活性化連鎖の増幅機構が働いて血小板凝集が促進される。このような血小板の活性化を惹起する物質は血流によって運ばれることから、赤血球や血小板が存在することによる流れ場への影響を反映した血管壁近傍での血流計算が必要になる。そのためオイラー型流体構造計算部分のプログラムについて、「京」において 25,600
ノードを使用して赤血球と血小板を含んだ血流計算の動作確認を実施した。
今後は、惹起物質の一つである ADP(アデノシン二リン酸)について、赤血球や血小板が存在する血管内での血流環境に、血管壁近傍での ADP 濃度分布形成がどのように関係するかについて調べるため、血管壁に粘着した血小板周りの流れ場の計算を行っていく。
（２）研究開発の実施状況
１）はじめに血栓形成の初期段階となる血小板の粘着・凝集の過程においては、血小板の粘着が解離可能

な不安定な一次凝集の段階から、粘着が解離しない非可逆的な二次凝集に至ることで血栓の形成が進展していく。血小板が血管壁に粘着して刺激を受けることで活性化され、さらに血小板の活性化を惹起する物質を自ら放出する。このような活性化を惹起する物質が、血小板表面に存在する受容体と結合することによって活性化連鎖の増幅機構が働いて血小板凝集が促進される。ADP は惹起物質の一つであり、光散乱法などによる血小板凝集能の測定実験により、ADP
濃度の違いによって一次凝集から二次凝集に至る過程と血小板凝集が解離する過程に分岐するような挙動を示す知見があり、ADP は血小板凝集において重要な役割を果たすものと考えられている。この ADP との結合を阻害して血小板の活性化を抑止することで血小板凝集を進展させない効果を狙った ADP 阻害薬としてクロピドグレルがあり、臨床でその効果が確認されている。より効果的な投薬治療を行うため、ADP 刺激がもたらす血小板凝集への影響と薬の作用機序についての評価が行えるシミュレーションモデルの構築が求められる。その実現のために，
まず ADP によって血小板が活性化される過程の計算モデル作りを行い、ADP による血小板の活性化増幅機構の挙動について確認を行っていく。
２）ADP 刺激による血小板活性化モデルとシミュレーション
血小板の表面に存在する P2Y12 受容体に ADP が結合すると、その結合の刺激を契機に血小板
内では、Gi(抑制性 GTP 蛋白質)により AC(アデニル酸シクラーゼ)の活性化が抑制、cAMP(環状アデノシン 1 リン酸)濃度が抑制されることで Ca2+(カルシウムイオン)濃度が上昇する。このような種々の化学反応の連鎖の結果として GP IIb/IIIa の活性化が誘導され、血小板同士の強固な粘着を促進することになる。さらに ADP などの活性化惹起物質を含んだ濃染顆粒の放出がおこり、血小板の活性化が増幅されていく。P2Y12 受容体と ADP との結合は解離して再び ADPと結合可能となり、ADP 刺激による血小板の活性化は断続的に生じるものと考えられる。このような血小板内で生じる化学反応は活性化シグナリングと呼ばれているが、実際に関与する物
質の特定や、それらが血小板の変化（血小板の変形なども含めて）にもたらす作用も多岐にわたることから、まだ多くの部分が解明されていない状況である。
(i)ADP 血小板活性化モデル
ADP 刺激による血小板活性化では、P2Y12 受容体が ADP と結合して血小板が活性化することで生じる GP IIb/IIIa 活性化と、ADP 放出(濃染顆粒)による活性化の増幅機構（ポジティブ・
フィードバック機構）が重要であり、この点に着目した次のようなモデルを考える。
 血小板表面に P2Y12 受容体と GP IIb/IIIa を配置する。
 P2Y12 受容体と ADP との結合を契機に、血小板内の Ca2+の貯蔵量を一定量増加させる。
 P2Y12 受容体と ADP との結合は、血小板近傍の ADP 濃度と所定の確率密度によって発生させる。（ここでは指数分布に基づいた結合確率を用いる）
 P2Y12 受容体と ADP との結合は一定時間を経過すると解離する。
 血小板内部に貯蔵された Ca2+濃度が一定量を超えると、GP IIb/IIIa を活性化させるとともに、一定量の ADP を放出する。
ADP 環境下に置かれた個々の血小板は、以上の動作を実行する。(図 2.1 参照)

P2Y12 activation：P>Prob(ADP)

ADP
ADP

P2Y12 receptor

Clopidogrel＋⊿ADP
ADP secretion
Receipted

Released

Timer starts

Timer stops

P2Y12-ADP
reception timer

Blocked

P2Y12-ADP
reception timer

P2Y12-ADP
reception timer

Count upOver Threshold＋⊿Ca2+

Time
counts

Dense
GranuleCa2+Store

GP activation：P>Prob(Ca2+)（Enpty）
GP IIb/IIIa
Activation timer

GP IIb/IIIa
Activation timer

Timer starts

Timer stops
Non-activation

activation

GP IIa/IIIa

図 2.1．活性化シグナリングを簡易化した ADP 血小板活性化モデル(ii)ADP 刺激による血小板活性化のシミュレーション
ADP 血小板活性化モデルを用いたシミュレーションにより、供給する ADP の初期濃度、血小板の配置・血小板数の違いによる血小板活性化の振る舞いについて考察する。GPII b/IIIa は血小板同士の強固な接着をもたらすことから、GP IIb/IIIa の活性化数は血小板の二次凝集に直接的に寄与する指標の一つと考えることができる。そのため、GP IIb/IIIa 総数の中で活性化した数の比率を GP IIb/IIIa 活性化率として、その経時変化をみることで ADP による血小板
凝集過程を比較する。
配置した血小板列の中央位置に一定値の ADP 初期濃度を与えて拡散させる条件とし、拡散す
る ADP 濃度は拡散方程式を解いて求める。実際の血小板表面に存在する P2Y12 受容体は 500～
600 個程度、GP IIb/IIIa は 4 万～8 万個程度あり、これらの蛋白質の分布には局在性がある
ことが示唆されている。ここでは各血小板について、500 個の P2Y12 受容体が 5 個単位で 100
か所に局在し、また 6 万個の GP IIb/IIIa が 200 個単位で 300 か所に局在しているものとして考え、一つの局在箇所内の蛋白質(P2Y12 受容体は 5 個、GP IIb/IIIa は 200 個)は同じ動作をするものとして粗子化した。血小板の配置については、PT1(100 個の血小板を整列配置)、
PT2(50 個の血小板をランダム配置)、PT3(一個分の間隔を空けて 25 個の血小板を整列配置)、
PT4(25 個の血小板をランダム配置)、の 4 つのパタンで比較する。
（図 2.2 参照）○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○(a) PT1

○○
○
○
○
○
○
○○
○○
○
○
○
○
○○
○○
○○○
○○○○
○○
○
○○
○○
○
○○
○○
○
○
○○○○○○
○○

○

○

○

○

○

○○

(b) PT2

○
○

○
○

○
○

○
○

○
○

○

○
○

○

○

○

○

○

○

○

○

○

○

○
○

○

○
○○

○
○
○

○○
○○
○
○

(c) PT3

図 2.2．血小板配置のパタン

○
○
○○○

(d) PT4○

図 2.3 は、各配置パタンついて、ADP の初期値を 5 つの条件 (0.1,0.3,0.4,0.5,1.0)で与えた場合の GP IIb/IIIa の活性化率の経時変化を比較した結果例を示したものである。
①ADP 濃度による GP IIb/IIIa 活性化への影響
GP IIb/IIIa の活性化率について、増加継続している場合は血小板凝集が進展していく過程
として、また減少減衰していく場合は血小板凝集が解離していく過程として考えることができる。(a)PT1 について見てみると、ADP 初期値が 0.4,0.5,1.0 の場合には、GP IIb/IIIa 活性化
率は時間経過とともに増加傾向を示しているが、一方、ADP 初期値が 0.1,0.3 の場合では、GP
IIb/IIIa 活性化率は、一旦は増加するものの、時間経過とともに減少している。同様に、他の配置パタンにおいても、GP IIb/IIIa 活性化率の値や継続時間に違いはあるものの、ADP 初
期値の違いによって、増加継続するケース(凝集)と減少減衰するケース(解離)に分岐する挙動を示すことが確認できる。
②血小板活配置による活性化への影響
今回の条件では、ADP 初期値が 0.3，0.4，0.5 の場合において GP IIb/IIIa 活性化率の遷移に顕著な違いが見て取れる。(a)PT1 と(c)PT3 は整列配置で血小板同士の間隔が異なる場合であり、PT1 に比べて血小板同士の間隔が広い PT3 の活性化率が下降傾向にある。(c)PT3 と(d)PT4は同じ血小板個数であるが配置が異なる場合であり、PT4 では血小板が近接している部分もあるが、大きく解離している(血小板 1 個分以上離れている)部分がある。PT3 に比べて PT4 の活性化率が下降傾向にある。このような結果から、血小板がより近接している方が血小板の活性化が増幅される傾向になることが示唆される。

ADP初期値=0.1

ADP初期値=0.3

ADP初期値=0.1

ADP初期値=0.3

ADP初期値=0.4

ADP初期値=0.5

ADP初期値=0.4

ADP初期値=0.5

ADP初期値=1.0

0.8

Rate of activated GPIIb/IIIa

Rate of activated GP IIb/IIIa

ADP初期値=1.0

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

50

100

150

200

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

50

100

Time stepTime step(a) PT1(b) PT2

150

200

ADP初期値=0.1

ADP初期値=0.3

ADP初期値=0.1

ADP初期値=0.3

ADP初期値=0.4

ADP初期値=0.5

ADP初期値=0.4

ADP初期値=0.5

ADP初期値=1.0

0.8

Rate of activated GPIIb/IIIa

Rate of activated GPIIb/IIIa

ADP初期値=1.0

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

50

100

150

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

200

0

50

100

Time stepTime step(c) PT3(d) PT4150200

図 2.3．血小板配置パタンの違いによる比較
３）まとめ血小板を活性化させる惹起物質の一つである ADP が、血小板表面に存在する P2Y12 受容体と結合すると血小板の活性化反応により血小板同士の強固な接着をもたらす GP IIb/IIIa が活性化する。その際、血小板から ADP が放出され、活性化の連鎖が増幅されて ADP ポジィティブ・
フィードバック機構が働く。この ADP 刺激による血小板活性化モデルを用いたシミュレーショ
ンにより、印加する ADP 濃度の値や血小板の配置・血小板数の違いによって、GP IIb/IIIa の活性化率が増加継続(凝集)と減少減衰(解離)のように分岐挙動となることを確認し、ADP 刺激による血小板活性化の増幅機構の特性を確認した。
ADP ポジティブ・フィードバック機構では、血管壁に粘着した血小板近傍での ADP 濃度分布が、血小板の活性化連鎖の増幅に重要な役割を果たすことが予想される。今後は、赤血球や血小板が存在する環境下での血流計算を行うことにより、血小板から放出される ADP と合わせて、血管壁近傍での ADP 濃度分布形成と血小板活性化との関係性の評価を行うとともに、有限差分法による流体構造連成計算手法を用いた血流計算と血小板内部の活性化の化学反応計算部分とを連動したシミュレーションモデルの構築を行っていく。
参考文献[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.

（Ｂ）「筋繊維の集合体としての骨格筋収縮のマルチスケールシミュレーション」（１）ソフトウェアの開発・高度化の状況
HI-Muscle は脊髄神経系モデルと筋骨格モデルを連成させることで実現する．脊髄神経系モデル
は各ニューロンの膜電位に関する連立常微分方程式とシナプス結合に関する数理モデルから成り，解析手法についてはほぼ確立されている．筋骨格モデルは非線形弾性体に対する 3 次元 FEM
モデルを筋肉に対して応用することで開発を進めているが，計算規模を大きくした際に期待通りの性能を得られていない面があり，改良を進めている途中である．Fig. 1 に 3 次元 FEM 解析の性能に関するテスト計算の結果を示す．並列数が大きくなったときに性能が頭打ちになっており，改良を必要とする．ただし，問題点は把握出来ており，早期の解決に向けて取り組みを進めている途上である．

Fig.1: 筋骨格モデル（3 次元 FEM モデル）のテスト計算結果．

（２）研究開発の実施状況
１）はじめにパーキンソン病のような神経変性疾患は日本ではその多くが特定疾患に指定されており，根本的な治療法も確立されていない．本研究では，スーパーコンピュータ「京」上に疾患状態の人体を再現し，健常者モデルとの比較等を通じて病因究明や新たな治療法の開発を目指す．
人体における運動は，大脳基底核から発せられた運動指令（電気的なスパイク信号列から成る）
が，上位神経系（主に脳）および下位中枢神経系（主に脊髄）における適切な修飾を経て伝達され，筋肉の収縮を引き起こすことで実現されている．そこで，神経変性疾患に伴う運動機能障害を再現するモデルを構築するにあたり，人体モデルを，脳モデル・脊髄神経系モデル・筋骨格系モデルの3 つに大別し，相互接続による統合を図ることを前提とした上でそれぞれのモデルの開発を進めている．本報告では，特に脊髄神経系モデルの開発について述べる．
脊髄は，脳から発せられた運動指令信号を筋肉に伝える役割を担うが，信号を脳から筋肉へと伝達するだけでなく，筋肉の状態（長さ変化や発揮している力）をモニターし，適切な運動を起こす，
あるいは，適切な力発揮を持続するための信号修飾の役割も持つ．そこには数多くの神経細胞が信号を相互に伝達し合う複雑なネットワーク回路が存在する．脳から発せられた信号が筋肉の収縮へと至る過程を正確に再現する上では，このネットワーク回路モデルの構築も重要な要素の一つであり，脳モデルから得られる信号を受ける，また，筋骨格系モデルへと信号を送出するなどの接続性を考慮することが必要である．

Fig. 2: 脊髄神経系ネットワークの構成要素

本研究課題に取り組むにあたり，前年度既報のとおり，まず初めに巨視的モデルによる脳・脊髄・関節運動から成るネットワークモデルを構築し，固縮時の運動ニューロンの活動度などの定性的傾向を把握することが可能となった．簡易的なモデルでありながら運動時の神経細胞の活動状態やそのときの相互作用を検討することができ，第一ステップとしては優れたモデルであった．しかしながら，モデル定数を多く含むため将来的な拡張性に乏しく，また，各神経細胞の活動度が厳密な意
味での物理量と必ずしも対応しない等，脳モデル・筋骨格モデルとの統合という観点からは必ずしも優れているとは言えない面も持っていた．そこで本年度は，コンダクタンスベースのニューロン
モデルに基づいた脊髄神経系モデルの開発に現在取り組んだ．このモデルでは，脊髄神経ネット
ワーク内で伝達される信号は，物理的に意味のある膜電位であるために，脳モデルから得られる運動指令との接続に関しても都合がよい．また，巨視的モデルとは異なり個々の神経細胞の状態（膜電位）を直接解析するため，巨視的モデルに比べてより詳細な筋線維の分布等と接続することができ，今後の拡張性も高いことが期待される．筋骨格モデルは，前年度までにふくらはぎに相当する部分のモデル化に取り組み，ヒトの足首関節の運動を再現するに至った．引き続きモデルの高精度化に取り組むとともに，脊髄神経系モデルとの統合にも取り組み始めた．本報告書ではそれらの成果の概要について述べる．
２）脊髄神経系モデルの構築
現在構築している脊髄神経系モデルの概要を Fig. 2 に示す．脊髄神経は筋線維の収縮を直接支配する  運動ニューロン（図中 M），フィードバック機構の 1 つである筋紡錘を支配する静的  運
動ニューロン（同じく S）および動的  運動ニューロン（D）と，3 種の介在ニューロン（Ia 介在ニューロン：I，Ib 介在ニューロン：X，レンショウ細胞：R）から成る．Ia 介在ニューロンは筋紡錘からのフィードバック信号を伝達し，Ib 介在ニューロンはゴルジ腱器官からのフィードバック信号を伝達する．また，レンショウ細胞は，多くの運動ニューロンに投射することで反回抑制や主働筋・拮抗筋の協調動作などの際に重要な役割を担う．このような神経系モデルを構築するにあたり，
解剖学的知見に基づいたコンダクタンスベースの数理モデルを採用する．たとえば，1 つのニューロンに対しては以下のように定式化される(1)．

CdV  I ion  I Leakdt

(1)

I L  g L VS  EL 

(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 

(3)
(4)

Fig. 3:  運動ニューロンの膜電位の経時変化．

このように 1 つのニューロンの膜電位は 4 つの常微分方程式によって記述され，それらを脊髄内に存在するニューロン数に応じて連立させて解く．
はじめに  運動ニューロンプールモデルを構築した． 運動ニューロンは筋線維を直接支配するニューロンであり， サイズの原理に基づき Rate coding と Recruitment と呼ばれる 2 つの機構によって発揮される筋力の調節機能を担う．Fig. 3 に  運動ニューロンの挙動についての結果の一例として直流電流を印加した際の膜電位の経時変化を示す．ここでは，ヒトの Soleus を支配する 運動ニューロンに関する解剖データに基づき各モデル定数を設定し，900 個の  運動ニューロンから成るプールを用いた．最小・中間・最大の各サイズの  運動ニューロンを抽出し，印加する
電流値を変えたときの挙動の変化を示している．図から確認できるように，印加電流が小さいときには，最小の  運動ニューロンしか発火しておらず，大きなサイズの  運動ニューロンは活動していない．印加電流を徐々に大きくして行くと，より大きなサイズの  運動ニューロンが発火するようになると同時に，既に発火していた小さいサイズの  運動ニューロンの発火周波数が高くなっていることがわかり，Rate coding と Recruitment が正しく再現できていることがわかる．このときに発揮される筋力の時間変化を Fig. 4 に示す．印加電流は脳からから発せられる運動指令に対応し，大きな力発揮に対応する運動指令ほど大きな刺激電流となる． Fig. 4 からわかるとおり，

Rate coding と Recruitment が正しく再現できていることで，刺激電流に対応して発揮される筋力が大きくなっていることがわかる．弱い刺激電流では最終的に発揮される力は小さく，また，その値

Fig. 4: 発揮される筋力の経時変化．

に達するまでに要する時間も長い．それに対して強い刺激電流では，大きな力が短時間で発揮されていることがわかる．
同じような数理モデルを脊髄内に存在する各ニューロンに対して適用し，脊髄神経系ネットワー
クを構築する．たとえば， 運動ニューロン‐レンショウ細胞間の反回抑制効果に関する結果をFig. 5 に示す．レンショウ細胞ならびに関連するシナプス結合の特性についても解剖で得られた結果に基づいたモデル定数を用いている(2,3,4,5)． 運動ニューロンとレンショウ細胞との発火頻度はほぼ線形の関係を持つことが知られているが，本モデルでも同様の結果が得られていることが確認できる．また，レンショウ細胞からの抑制信号によって  運動ニューロンの活動が抑制され，発揮される筋力が最大 25% 程度低下することも確認できた．

Fig. 5: 左：a 運動ニューロンとレンショウ細胞の発火頻度の関係．右：レンショウ細胞による発揮
される筋力の抑制効果（反回抑制）．

３）運動ニューロンの発火と筋骨格系有限要素との連成モデル

運動ニューロンの連続した活動電位は筋線維に伝達され，筋線維は収縮力を発揮する．運動ニューロンの単一の発火に対する筋線維の力学的な応答は単収縮である．このとき，1 つの筋線維が発生する力は次式で与えられる(6)．
f i t  Pi  t 1t Ti 
e
Ti

(5)

ここで，Pi および Ti はそれぞれ，単収縮における収縮力および収縮時間であり，筋線維ごとに異なる値をもつ．本研究ではそれぞれ，次式のモデルを用いた．
Pi  eln 100 
i
n

(6)
1 4.2

1
Ti  90  Pi 

(7)

筋線維の力の応答は運動ニューロンの発火速度に対して非線形の関係にあり，その力の挙動は活動電位の間隔（ISI）の関数として記述され，上式に導入される．活動電位 j に対するモーターユニッ
ト i に所属する筋線維が発揮する力 fi, j は次式で与えられる．


2 T ISI

1 e i j
gi, j Ti ISI j3(8)

最終的に，モーターユニット i の連続的な発火に対する筋線維が発揮する力 Fimotor は次式で与えられる．
Fi motor t    f i , j t  ti , j kj 1

0  t  ti , j(9)Fig. 6: 脊髄神経系モデル・筋骨格有限要素モデル連成解析のテストモデル．

ここでのテスト解析に使用したモデルを Fig. 6 に示す．各短辺（x, z 方向）を 100 等分し，各要
素が 1 本の筋線維に相当すると仮定した．筋線維方向は y 方向に一様とし，1 本の筋線維は長手方向に連続的につながっているものとする．モーターユニット数は 900 とし，各運動ニューロンが支配する筋線維を 3 次元モデル内にランダムに定義することで，3 次元空間内の筋線維と運動ニューロンを連携させた．運動ニューロンの発火は本節前半の脊髄神経系モデル内で述べられている手法
によって解析している．

脊髄神経系との連成モデルによる等尺性収縮シミュレーションの結果の一例を Fig. 7 に示す．
本モデルでは，刺激入力の時間的増加に従い，低閾値の運動ニューロンから順に動員され，それに支配されている筋線維も順に筋活動を始める様子が確認できた．各運動ニューロンが支配する筋線
維は 3 次元空間内でランダムに設定されているため，筋内の筋活動は一様ではなく，ある時刻におFig. 7: 脊髄神経系・筋骨格有限要素モデル連成解
析の結果例．図は端面に掛かる力の方向ベクトルとその大きさ（色）．

Fig. 8: ヒトの足首関節に関する骨格筋モデル．

いては筋活動が生じている筋線維と生じていない筋線維が混在する．筋線維が発揮する力はそれを支配する運動単位によって異なるため，空間的に不均一な分布を持つことが確認できる．現在このモデルを Fig. 8 に示すようなヒトの足首関節に適用することに取り組んでいる．
４）まとめと今後の計画
本稿では，HI-Muscle の現状について報告した．ソフトウェアの開発状況という観点からは，現在十分な性能が得られていない面があるものの，問題点は把握できており早期に解決出来るものと考えられる．研究開発状況という観点からは，フィードバックシステムがそれほど重要ではない状
況（パーキンソン病の安静時振戦を対象とした解析など）については脊髄神経系・筋骨格系の連成に目途がついており，脳から来た運動指令に応じて脊髄におけるニューロンの活動と筋肉の収縮挙動を把握可能である．今後は，筋紡錘・ゴルジ腱器官とそれらが結合する介在ニューロンの実装に取り組むとともに，医療画像データなどに基づいた詳細なヒトの骨格筋モデル（筋線維分布など）の構築に取り組む予定である．
参考文献
(1)

Cisi, R. R. L., Kohn, A. F., Simulation system of spinal cord motor nuclei and associated nerves and muscles, in a Web-basedarchitecture, Journal of Computational Neuroscience, 25(2008), 520-542.

(2)

McCurdy, M. L., Hamm, T. M., Spatial and temporal features of recurrent facilitation among motoneurons innervatingsynergistic muscle of the cat, Journal of Neurophysiology, 72(1994), 227–234.

(3)

Windhorst, U., On the role of recurrent inhibitory feedback in motor control, Progress in Neurobiology, 49(1996), 517–587.

(4)

Cullheim, S., Kellerth, J. O., Morphological-study of axons and recurrent axon collaterals of cat alpha-motoneurones supplying
different functional types of muscle unit. Journal of Physiology, 281(1978), 301–313.

(5)

Sussman, T., Bathe, K.-J., A FE Formulation for Nonlinear Incompressible Elastic and Inelastic Analysis, Computer &Structures, 26(1987), 357-409.

(6)

Fuglevand, A. J., Winter, D., Patla, A. E., Models of recruitment and rate coding organization in motor-unit pools, Journal of
Neurophysiology, 70(1993), 2470-2488.

Ⅲ－２

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

Ⅲ－２－１

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

Ⅲ－２－２ 実施内容（成果）（１）ソフトウェアの開発・高度化の状況
１）モデル統合インターフェイス「脳神経系の多階層モデル」、「脊髄と四肢筋骨格モデル」、「筋繊維レベルからの骨格筋モデ
ル」を統合するインターフェイスの開発を行ってきた。東京大学の高木チーム、沖縄科学技術
研究基盤整備機構の銅谷グループと議論を行いながら、我々が開発してきた全身筋骨格シミュ
レータ、高木チームの詳細筋有限要素シミュレータ、銅谷グループの脳中枢神経系シミュレー
タの統合インターフェイスを設計した。インターフェイスの概要は図１のようになる。
統合インターフェイスは、ニューロシミュレータ NEST[R1]、骨格筋有限要素解析シミュレー
タ V-Biomech[R2]、全身筋骨格シミュレータ sDIMS[R3]、大規模なプロセス間の通信を可能にする MUSIC[R4]から構成される。骨格筋有限要素法計算と全身骨格モデルの剛体力学シミュ
レータの開発環境は、前者が Fortran を後者が C++を基盤とする。これらをつなぐソフトウェアの開発を行った。両環境において相互に参照可能であり、剛体のマクロな力学情報および筋
有限要素全体の境界条件のやり取りを行い、通信情報量の低減による並列計算効率の向上を行っている。一方、筋骨格モデルの入力信号は連続値として表された筋活動度である。これに対して、一般的に神経活動は離散的なスパイク信号として表されることが多い。文献[R5]と[R6]を基盤として、スパイク信号と連続波形との変換を行うソフトウェアを開発した。筋有限要素法が利用する線形方程式のスパースソルバー、および NEST で利用される python インターフェイスを除き、全てのソフトウェア環境が「京」で利用可能な段階となっている。

図１．モデル統合インターフェイスの概念図
２）四肢筋骨格モデルの開発
計算論的神経モデルから人間の身体の物理モデルへ結び付けるための人間の解剖学的幾何モデルとして、ライフサイエンス統合データベースセンターの BodyParts3D [R7]の人間の全身のポリゴンモデルを採用した。このモデルは人間の幾何形状情報として、詳細な骨格、筋、などの運動器の他に、神経系、循環器系、消化器系などを含んでいるものの、機能情報を持たない。このため骨格上の筋腱複合体の端点、求心性・遠心性の運動神経と筋との接合点などの機能情報を解剖学の教科書を見ながら追加する作業を行う必要がある。一方でこのモデルは詳細な筋表面情報を持ち、表面の筋繊維配向については推定可能である。表面上の筋繊維配向から筋繊維配向ベクトル場を設計し、内部の筋繊維配向を補完可能なインターフェイスを作成し、下肢大腿の骨格モデルおよび１７個の筋の有限要素モデルを作成した（図２）。計測データ等による内部筋繊維配向情報も今後は取り込める形で開発を続ける。

図２．下肢大腿部の筋有限要素モデル(左：ポリゴンモデル、中：設計した筋繊維ベクトル場、右：有限要素モデル)（２）研究開発の実施状況
１）体積効果を考慮した四肢筋骨格シミュレーション
前年度から筋肉の体積効果を考慮した全身の筋張力の並列計算手法(図３）の改良を行った。
本手法では、ロボティクスにおける剛体リンク系の力学計算法と、筋形状モデルに対する非線形有限要素解析を同時に解く。筋繊維の三次元的な配向や骨接触部の面圧力などの詳細な力分布を再現できる。筋間の接触力やそれに付随する筋の変形も将来的には統合可能である。逐次二次計画法を応用することで、筋の有限要素計算を剛体リンク系の力学計算と独立に分離することができるため、一つ一つの筋の有限要素計算を並列に計算することができる。また一つの筋の有限要素計算単体において，有限要素解析の並列計算を利用することも可能である。

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

高木チームが開発した骨格筋有限要素シミュレータを改変せずそのまま利用可能なインターフェイスを開発した。微圧縮性を考慮した Mooney-Rivlin 体を利用する計算手法へと変更されるため、弾性体の変位の他に圧力も未知変数となるが、昨年度までに開発した手法と同様の定式化ができることを示し、矛盾のない拡張を行った。
下肢の骨格モデルおよび膝関節回りの 6 つの主要な筋の有限要素モデルを利用して数値実験を行った(図 4)。簡易モデルでは筋数が 6 つのため、6CPU を用いて MPI 並列計算を行った。作
成したモデルを用い、筋活動推定の際の二次計画問題を並列計算と逐次計算の二種類とで実行
して、実行時間を比較した。逐次計算、並列計算での実行時間はそれぞれ 5.03×103 秒，1.067
×104 秒であり，並列計算により逐次計算の約 3.4 倍の実行速度が得られた。また得られた 6
つの筋活動度は、解析対象とする運動中の膝関節の関節トルクを再現できた。これによって並列性が高い計算アルゴリズムであることが確認できた。全身筋骨格モデルにおける高速計算に貢献するものと期待される。筋有限要素法が利用する線形方程式のスパースソルバーを除き、四肢筋骨格シミュレータは「京」で利用可能な段階となっている。

図４．下肢モデルの筋有限要素計算の様子
２）脊髄の筋支配モデルの開発
大規模脊髄神経ネットワークの構築を行い、骨格筋と脊髄神経モデルの接続を行い、神経筋
骨格シミュレータを開発した。脊髄神経ネットワークは NEST を用いて構築し、筋骨格シミュ
レータにはワイヤ駆動型筋骨格モデルを利用した。ワイヤ駆動型筋骨格モデルは後に有限要素
筋骨格モデルに置き換える予定である。文献[R5]と [R6]を基盤として、脊髄神経ネットワー
クと筋骨格モデル間の神経信号のやり取りにおいて、運動ニューロンのスパイク信号から筋活
動度へ、筋長変化から感覚ニューロンのスパイク信号へ変換を行うインターフェイスも同時に開発した。これらは MPI によって実装され並列計算が可能である。
上肢の筋骨格モデルと脊髄神経ネットワークを接続し、伸長反射シミュレーションを行った。
脊髄神経系の設計は文献等を参照しながら作成し、ニューロン数は約 500 程度、ワイヤ数は約300 程度となる。脊髄神経系および筋骨格系で合わせて 4CPU を用いて MPI による並列計算を行った。
並列計算時には、
筋骨格力学計算に 1CPU、NEST での脊髄神経シミュレーションに 2CPU、通信部である MUSIC に 1CPU を割り当てた。逐次計算、並列計算での実行時間はそれぞれ 361
秒、279 秒となった。現状では力学計算と脊髄神経シミュレーションを交互に行って通信を行い、また各計算時間が大きく異なっている。計算効率の向上は課題として残されているが、外乱を加えた際に上腕の筋の伸長反射を再現することができた。NEST で利用される python イン
ターフェイスを除いて、脊髄神経筋骨格シミュレータは利用可能となっている。

図４．上腕の伸長反射シミュレーションの様子

図５．運動(MN)・感覚(SN)・介在(IN)ニューロンのスパイク信号時系列。
赤丸付近で下腕に対して鉛直下方向へ外乱を加えた。
参考文献[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] RRL, Cisi, AF, 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.
[R6] 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.
[R7] BodyParts3D, ライフサイエンス統合データベースセンター,
http://lifesciencedb.jp/bp3d

Ⅲ－３

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

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

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

実施内容（成果）（１）ソフトウェアの開発・高度化の状況
１）このプロジェクトの開発言語である PyNEST の実行に必要な Python が、現時点では「京」でサポートされていないので、「京」に Python の実行環境を構築した。Batch モードではジョブを投げる度に、Python のソースをコピーしてコンパイルする必要があるため、実行までに１時間近くの時間がかかるという問題があったが、共有ディレクトリを使用することによりBatch モードと Interactive モードの両方で Python が使用可能となった。しかしながら、
PyNEST
上で OpenMP が使用できない問題とファイル出力でセグメンテーションエラーが出るという問題が現在でも未解決のまま残っている。
上丘モデル（文献１）と Shouno et al.
モデルを「京」で走らせベンチマークを行った。その結果 1024 コアまでの線形なスケーリングが確認された（図１）。また、どちらのモデルを使った場合にも 1024 コア以上では線形ス
ケーリングが見られないことから、
NEST2.2.1 と「京」の間に何らかの問題があることをつきとめた。
２）実験データと文献による知見を融合できるパラメータ推定法を提案した。この手法ではモデルで再現したい
図１：NEST による上丘モデル（文献１）の「京」実験データのスパイク統計量等の分上でのスケーリング
布を目的分布として作成することにより実験データに含まれる確率的なゆらぎをとらえることができる。そして、NEST 等のシミュレータを用いて、あるパラメータ値
におけるモデルの振る舞いを、注目する統計量の経験分布として記述する。提案手法では、目的分布と経験分布間の KL 偽距離をエネルギーとしたボルツマン分布からメトロポリス法を用いてサンプリングすることにより、目的分布を実現するようなパラメータの分布を求めることができる。この提案手法は、状態遷移カーネルにパラメータに関する事前知識を埋め込むことができるので、文献などでパラメータに関する知見などが得られる場合は、その知見を自然な形で推定に埋め込むことができる。また、この手法は高次元パラメータ空間の探索を効率的に行うために、並列計算との親和性が高いレプリカ交換モンテカルロ法と組み合わせられることも確認した。
（２）研究開発の実施状況
１）大脳基底核モデルの構築
パーキンソン病モデルのサルの神経活動データ（図２）を、自然科学研究機構生理学研究所の橘吉寿助教に提供していただき（文献２）、L-DOPA 投与前後の神経活動の特徴量の変化について解析を行った。その結果、条件によって変化する神経活動データの確率的な振る舞いを分布として捉え直し、その分布を再現するようなモデルの構築やパラメータ推定法が必要である
ことが明らかになった。
投与前

Dyskinesia

on drug

off drug

図２：サルの淡蒼球外節(GPe)のスパイク記録データ (Tachibana et al.)。
L-DOPA の投与により発火パターンが変化する。

その結果を踏まえ、Kumar et al.の視床下核（STN）-淡蒼球外節（GPe）モデル（文献３）とパラメータの分布を自動推定する提案手法を用いて、パーキンソン病に特有な特徴量の再現を試みた（図３）。その結果、ドーパミン枯渇が起こった時に顕著になるベータ波を再現するようなパラメータの分布を自動的に発見することができた。
また、Shouno et al. (2009)モデル（文献４）を拡張し、通常時におけるベータリズム生成の問題を解決した。具体的には、視床下核(STN)から淡蒼球外節(GPe)への興奮性シナプスの短期可塑性のパラメータを変更し、視床下核神経細胞の発火頻度が低い状態における視床下核か
ら淡蒼球外節への興奮性入力を減弱した。その結果、視床下核と淡蒼球外節の間のカップリン
グ強度が低下し、それがベータ派帯域の振動パワーの低下につながった。また、NEST2.0 を使用していた Shouno et al. (2009)モデルを、「京」などの大規模並列計算でのメモリ使用を大幅に改善した NEST2.2 に移植した。

図３：Kumar et al. モデルのシミュレーション結果。視床下核(STN:緑)と淡蒼球外節(GPe:青)の発
火が、ドーパミン枯渇を模した線条体からの抑制入力の変化によりβ帯域での振動が見られる。

参考文献
1) Morén J, Shibata T, Doya K. (2013). The mechanism of saccade motor pattern generation
investigated by a large-scale spiking neuron model of the superior colliculus. PLoS One.
8(2):e57134
2) 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
3) Kumar, A., Cardanobile, S., Rotter, S. and Aertsen, A. (2011). The role of inhibition in
generating and controlling Parkinson’s disease oscillations in the basal ganglia. Front. Syst.
Neurosci. 5:86.
4) 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

Ⅲ－４

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

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

Ⅲ－４－１

から開発が進められている「全身筋骨格モデルにおける筋肉の体積効果の大規模並列化及び人の行動モデルに関する大規模計算」に導入可能な脳神経系による運動制御機構のマクロ機能レベルの数
理モデルを構築する。特に、ヒト静止立位姿勢と二足歩行運動の神経制御メカニズムをモデル化の対象とする。また、運動障害を伴う中枢神経疾患のひとつであるパーキンソン病の病態を、本研究で構築する制御メカニズムの崩壊現象として定量的に説明・再現する。
また、「コンプライアントな生体筋・腱系に駆動されるヒト骨格系の運動制御に対して大脳基底核が果たす役割の数理モデル化」の研究を行う上で、関連する研究者と必要な協議等を行う。
平成２４年度は、（１）東京大学・高木グループの「筋繊維レベルからの骨格筋モデル」、沖縄科
学技術大学院大学学園・銅谷グループの「脳神経系の多階層モデル」、東京大学・中村グループの「脊髄と四肢筋骨格モデル」を統合したインターフェイスの利用を可能にする。
（２）運動に関わる関節の受動的粘弾性が低くかつ能動的神経フィードバック制御のゲインも小さい（関節がコンプライントである）という生理学的に妥当な条件の下で全身機械力学モデルの直立姿勢と歩行を柔軟に安定化する制御機構を数理モデル化する。
Ⅲ－４－２ 実施内容（成果）（１）研究開発の実施状況

平成２４年度は、ヒト静止立位姿勢制御に関して、関節の受動的粘弾性が低くかつ能動的神
経フィードバック制御のゲインも小さい（関節がコンプライアントである）という生理学的に妥当な条件の下で全身機械力学モデルの直立姿勢を柔軟に安定化する制御機構の数理モデル化を実施し、東大高木・中村グループ、OIST 銅谷グループが構築する大規模並列モデルが実現すべき身体運動動態のマクロレベルの特性を明らかにした。以下具体的に成果を述べる。今年度達成した成果は 2 つである。
１）成果１（発表論文１）
この研究では、ヒト静止立位時の身体力学系を下肢と頭部・体幹・上肢に対応する剛体 2 リ
ンクからなる倒立二重振子でモデル化した。この際、2 リンク系の 2 つの関節は足関節および股関節に対応する。立位時の両関節はコンプライアントであるという生理学事実に基づき、このモデルにおける両関節の受動的粘弾性は小さな値を仮定した。このような身体力学系を制御対象とした倒立二重振子の立位姿勢制御の時間遅れを伴う神経フィードバック制御を検討した。特に、フィードバック制御戦略として、従来仮説である持続的比例・微分制御、および我々が提唱している間欠的フィードバック制御による姿勢制御を比較した。その結果、持続的比
例・微分制御器が静止立位を安定化できるゲインパラメータの範囲は極めて限定的であり、従来仮説はヒト脳神経系の姿勢制御戦略として妥当ではないことを示した。一方、間欠的比例・微分制御器は、広い範囲のフィードバックゲインに対して、静止立位を安定化できることを示した。従って、このことは比例ゲインおよび微分ゲインが小さな値であっても、すなわち足および股関節がコンプライントであっても、能動的な神経フィードバック制御を適切なタイミン
グで間欠的に停止・活性化する間欠フィードバック制御を用いれば、ヒト静止立位と同様な柔軟性を実現しつつロバストに直立姿勢を安定化できることを意味する。

２）成果２（発表論文２）
この研究では、ヒト静止立位時の重心動揺の生成メカニズムを検討した。イタリアの
Corforto ら(2001)は、ヒト静止立位時の重心動揺と心拍を同時計測することで、重心動揺時に足裏に作用する鉛直方向の床反力および剪断応力を計測し、これらが心臓拍動に起因する循
環動態に同期したゆらぎを示すことを明らかにした。このゆらぎの振幅および力の作用点に関する簡単な力学計算により、心臓拍動が立位姿勢時の足関節に与える力学的摂動の大きさは0.04 Nm 程度であることを明らかにした。我々は、もし立位姿勢が柔軟に安定化させておれば、
この循環動態に起因する極めて微小な摂動であっても、実験的に観測される程度のヒト静止立位時姿勢動揺を再現できるのではないかと考えた。そこで、立位時の身体を単一倒立振子でモデル化し、その神経フィードバック制御（時間遅れを伴うフィードバック制御系）を検討した。
特に、フィードバック制御戦略として、従来仮説である持続的比例・微分制御、および我々が提唱している間欠的フィードバック制御による姿勢制御系を構築し、これらのモデルの足関節に、振幅が 0.04 Nm で周期が 1.2 Hz（心臓の拍動周期）の機械的摂動を加えたときの姿勢動
揺動態を比較した。その結果、持続的制御によって剛直に安定化された立位姿勢制御モデルでは、実験的に計測されている姿勢動揺は再現できなかった。一方、間欠制御によって柔軟性を保ちつつ安定化された立位姿勢制御モデルでは、実験的に計測されている姿勢動揺を再現できた。このことは、ヒト立位姿勢は間欠制御によって達成されおり、立位時の姿勢動揺の起源は、心臓拍動に伴う心循環系の機械力学的作用であることを示唆する。
これら 2 つの成果は、マクロレベルにおけるヒト立位姿勢制御の特性と制御様式を理論的・定量的に明らかにするものであり、本プロジェクトで開発中のマルチスケールなヒト身体運動
制御系が再現すべき現象とその背後にある制御メカニズムを提供するものと考えられる。

発表論文リスト
1. Suzuki, Y., Nomura, T., Casadio, M., Morasso, P. Intermittent control with ankle, hip, and mixed
strategies during quiet standing: A theoretical proposal based on a double inverted pendulum
model (2012) Journal of Theoretical Biology, 310, pp. 55-79.
2. Nomura, T., Oshikawa, S., Suzuki, Y., Kiyono, K., Morasso, P. Modeling human postural sway
using an intermittent control and hemodynamic perturbations (2013) Mathematical Biosciences, .
Article in Press.