学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

11NA-27

毛細血管流までを再現するマルチスケール冠循環シミュレーション
岩村 尚(東京大学), 鷲尾 巧(東京大学), 岡田 純一(東京大学),
中島 研吾(東京大学), 片桐 孝洋(東京大学), 久田 俊明(東京大学)
概要 当研究チームが開発している心臓シミュレータでは, 実測データに基づいた冠循環数理モデルを組み込み, さらに冠循環血管系の解剖学的な特徴を利用することで, 均質化法によるマルチスケール解析と同様に, 超並列性を活かした毛細血管レベルまでの冠循環シミュレーションが可能である. 本公募型共同研究においては, 冠循環マルチスケールシミュレーションの高度化と超並列計算に対応した幾つかの最適化を重点的に実施する.

1. 研究の目的と意義いて血管の分岐・集合形態を再現した冠循環数理「ヒトゲノムプロジェクト」の完成をうけて,モデルを作成し, この全 23 階層ネットワークを生命科学の焦点は遺伝情報の解明からタンパク,UT-Heart に組み込んだ. 冠循環血管系の解剖学的細胞内小器官, 細胞, 組織, 臓器, そして個体と特徴を利用することで, 均質化法によるマルチスいう多階層に亘る複雑な生命現象の理解へと移行

ケール解析と類似の数理的アプローチが可能となしつつある. 第三の実験(in silico 実験)として期り, 超並列性を最大限に活かしたシミュレーショ

待されている計算機シミュレーションにおいてもンを行うことが出来る. 最終的にはペタスケール

マルチスケール解析の重要性が増している. 当研

コンピュータ上で精密かつ大規模なシミュレーシ

究チームが開発を進める UT-Heart は, 細胞のイオョンを行い, 冠循環と心拍動の関連, 虚血性心疾

ンチャンネルや収縮タンパクの数理モデルから出患の病理解明につなげる.

発し有限要素法でモデル化された心室の収縮, 血本公募型共同研究においては，酸素の輸送と心液の拍出に至る現象を一貫して再現できる国際的筋組織への供給，Cortassa ら 5)の研究を拡張したにも突出した心臓シミュレータである.

代謝数理モデルの組み込み等，冠循環マルチスケ

UT-Heart では冠循環のシミュレーションも可能

ールシミュレーション手法の高度化を図ると共に，である. 冠循環とは心臓自体に血液を供給する役

計算科学分野の研究者との密接な協力によって，割を担う血管系である. 冠循環系における流れは超並列計算に対応した幾つかの最適化を重点的に心拍動に伴うメカニカルストレスの影響を強く受実施する.

け, 他の臓器には無い特有の血流パターンを示す.
例えば, 心室の収縮期において体循環は順流とな

2. 当拠点公募型共同研究として実施した意義るのに対し, 冠動脈の流れは殆どなく, 一部では

(1) 共同研究を実施した大学名と研究体制逆流すると言われている. しかし in vivo 計測の東京大学困難さから, 特に毛細管を含む微小循環系での挙(2) 共同研究分野

動は現在も明らかではない. 冠循環の異常に起因

超大規模数値計算応用分野

する虚血性心疾患は心疾患による死亡原因の約半(3) 当公募型共同研究ならではという事項など数を占めるため, そのメカニズム解明は大きな医

超並列計算機上での実行に向けた準備として,学的意義を持つ. 当研究室では，Kassab ら 1-4)にサイズを縮小した問題や, 粗い精度のシミュレーよる実測データに基づき毛細血管まで含むすべてションを数多く実施することにより,できるだけのスケール（動脈 11 階層, 静脈 12 階層）にお多くの問題を事前に解決しておく. 本公募型共同1

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

Fig. 1 冠循環マルチスケールモデル

研究においては，計算科学分野の研究者との密接以下に相当する. 心表面を走行する最も太い血管な協力によってそのような準備を実施する.

ネ ッ ト ワ ー ク を CT 画 像 か ら 作 成 し , 続 い て

3. 研究成果の詳細と当初計画の達成状況

Kassab らの計測データに基づき, 心表面冠血管(1) 研究成果の詳細についてから派生する血管の径, 長さ, 分岐を決定し, 同心筋を栄養する冠循環の毛細管レベルまで詳細

時に心臓有限要素両心室モデルにマッピングした.に取り扱ったシミュレーションは, それら全てを以下, これらをマクロモデルと呼ぶ. 一方, 図パイプ要素で表した際の膨大な計算量により, こ左に示す, 毛細血管を含むレベル 3～-3 のネッれまで行われてこなかった. 我々は冠循環血管系トワークを微小循環モデルと呼ぶ. さらに図中の解剖学的特徴を利用することで, 均質化法によ央のレベル 5～3, -5～-3 の血管には, Kassab

るマルチスケール解析と類似の数理的アプローチに倣って対称性を仮定する 6). 親血管から分岐を用いた, 並列性を生かした冠循環マルチスケーする２つの子血管の血管径, 長さ, その後の分ル解析手法を開発した.

岐形態は全く同一のものを与える. このような3.1. 冠循環マルチスケールモデル

対称性の仮定により, 血管ネットワークの本質

Kassab らの計測結果 1-4)に従い, 毛細血管を的な特性を失うことなく, 計算量を大幅に削減含む全スケールの血管にレベルを与え, 大動脈かすることが可能となる. 以下, 微小循環モデルら発生する太い冠動脈をレベル 11 とし, 分岐をと対称性モデルを合わせてミクロモデルと呼繰り返しながらレベル 1 の細血管, さらにレベル

ぶ.

0 の毛細血管に至ると考える. 毛細血管はレベル

マクロモデルとミクロモデルの接続点は仮想的-1 の細静脈に接続し, 集合を繰り返しながら最終に心筋の圧力節点と同じ位置にあると考える.

的にレベル-12 の太い冠状静脈となり右心房へ接これによりミクロモデルが心筋から受ける外圧を続する.

心筋節点圧 で定義でき, またミクロモデルの接我々が提案する冠循環マルチスケールモデルを続数も心筋節点集中体積 より決定できる

Fig. 1 に示す. 図右に示した 2 つのネットワークので扱いが容易になる. さらに心筋容積に対するモデルは動脈側レベル 6 以上, 静脈側レベル－6

細血管の密度の貫壁性による違いを考慮するため,
2

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

心外膜側で 1, 心内膜側で 2.2 となる血管密度比を心筋節点に与える. ミクロモデル血管の総容積を心臓容積の 8%とすると, 心筋節点に単位体積当たりに含まれるミクロモデル血管容積が与えられる. ミクロモデル 1 個の血管総容積を とすると, マクロモデルの末端点に接続されるミクロモデルの数 は次式で決定できる.
   

節点においては次の流量保存則が成立する.
たFig.
4 ポアズイユ流れを仮定した 1 次元有限要素の取り扱い
だし節点に接続している節点の集合を, とする.

マクロモデルとミクロモデルの動脈側, 静脈側接-  !"  

続部において次式の流量保存則が成立する. Q は

./0

流量を, γは動脈または静脈を表す.
ミクロモデルの血管の式

      

ミクロモデルの血管は心筋を直接栄養する役割

以上の説明図を Fig. 2, Fig. 3 に示す.

を担っており, 血液の輸送に加え血管体積変化による血液の貯蓄や放出を考慮する. 一般に血管要
素の 2 端点で圧力値が異なり要素内でも径が異なるため, Fig. 4 の様に血管要素を 2 つの分割要素とで分けて考える. 分割要素の内圧は* , * で一定と考えると, 血管径は内圧（冠血圧）と外圧12（心筋圧：心筋有限要素の圧力節点の値に相当）の差により次のように決定される. 3は血管のFig. 2 マクロモデルとミクロモデルの接続弾性を, *は基準状態を示す.
6 "7$ !"  45 * !" + 12 " + *6 + 12 $6

分割要素のコンダクタンスは便宜的にポアズイユ流れより得られる値で代用する. よって血管要素のコンダクタンス8 !"は分割要素のコンダクタンス8 !", 8 !"の調和平均として次式のように求められる.
8 !"8 !"$ !"% $ !"%#8 !" 8 !"  8 !" 9:*) $ !"%  $ !"%Fig. 3 与えた血管密度比

また, 分割要素が溜め込む血液は体積の時間微分

マクロモデルの血管の式で表すことができ（次式第二項）, ミクロモデル

マクロモデルの血管は心臓隅々まで血液を運搬の血管節点においては次式の流量保存則が成立する役割を担っており, Kassab らに倣い血管内流する.

れはポアズイユ流れを仮定する 4). 節点,を持- 8 !" * + * "  -つ要素の節点からへの流量を と表すと次式に./0より求められる. は血管径, は長さ, は粘性,

は節点,の圧力を表す.

./0#) ; $ !"% "(;!

心臓有限要素と血液からなる混合体の基礎式#$% !" * !" + * !""&'( ) 心筋と微小循環系に含まれる血液からなる混合体3

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月の運動方程式は次式で与えられる.
<  = > ? + 1@"

Aは心筋内圧，は混合体密度,<は加速度, ? はCauchy 応力テンソルを示す．体積変化率Bを用いた混合体の体積保存則は以下となる.
B + &" + CDE   
第二項がミクロモデルの体積変化であり, ミクロモデルに蓄えられる血液が増加した分だけ混合体
Fig. 5 冠動脈および冠静脈流量の容積が増加することを表している. DE  はミクロモデル 1 個の初期状態からの容積増加量, nは単位心筋体積当たりに含まれるミクロモデルの数である.
結局, 未知数は心筋の変位 3 成分, 心筋圧力１成分, 血管圧力１成分の計 5 種となり, これらは以上の式を連立して解くことで求められる.
計算手法
マルチスケール解析の考え方を用いて並列計算により解く. 全体のマトリクスは次式で与えられる.
GF HM D*HIK F   K  LM N
J D*G  は血管ミクロモデルより作成される係数マトリクス, J は血管マクロモデル及び心筋モデルより作成される係数マトリクスである. Hはミクロ部とマクロ部を接続する血管の要素マトリクスである.

上式の係数マトリクスは以下の様に分解できる. GF HG HIKFJH

IP GRS
2Q H KKFJO PI
ただし, JO  J + HGRS
2Q H

上式に基づき以下の 3 ステップでミクロ, マクロの増分を計算する.
1. G  D*T2Q  M2Q よりD*T2Q を計算

2. JO D*  M + HD*T2Q よりD*を計算

I
3. D*   D*T2Q + GRS
2Q H D* よりD*  を計

Fig. 6 小血管および毛細血管の流量同士は完全に独立しているため並列計算により効算
以上の様に, マクロのマトリクスに入り込むミク率的な解法が実現できる.

ロの自由度を消去することで, マクロとミクロを3.2. 解析結果分けて計算することが可能となる. さらにミクロ
4

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

Fig.7 マクロ・ミクロプロセスの同時実行

冠動脈への血液流入量および冠静脈からの流出（3）拡張期初期の逆流, が現れている. また貫壁量の経時変化を Fig. 5 に示す. 血液の流出入量は性の影響も再現されており, 収縮力が強い, すな

共に 151.7ml/min となり, 心筋容積に応じた安静

わち心筋内圧の高い心内膜(Endo)側ほど, 動脈の時での血液量の範囲内(147.8～157.4 ml/min)で逆流が強く現れていることがわかる. さらに毛細ある. 動脈では収縮期における 2 相性の血流量の血管内の流れにおいても貫壁性の違いが現れてい減少と拡張期優位な流れが, 静脈では収縮期優位る.

な流れが再現できていることが分かる.

3.3. マクロ・ミクロ計算の一部同時実行

次に, 小動脈(レベル 5), 小静脈(レベル-5)お冠循環のマルチスケール解析においては, マク

よび毛細血管の流量を Fig. 6 に示す. 小動脈ではロ部(心筋, 内腔流体, 冠循環マクロモデル)を担（1）収縮期初期の一瞬の逆流, （2）拡張期初期当するプロセスと冠循環ミクロモデルを担当するの順流のピーク, 及び（3）拡張期優位な流れ, がプロセスを分離することで, 一部の計算を同時実表されている. 一方で小静脈では（1）収縮期初期行することが可能である(以下, それぞれマクロに見られる順流のピーク, （2）収縮期優位な流れ,プロセス, ミクロプロセスと呼ぶ). Fig.7 は一部同時実行時のダイアグラムである. ミクロプロセ
ス側の行列G2Q の計算およびミクロモデルの求解と, マクロプロセス側のU2VQ 計算が同時実行可能である.
T2K(東大)で性能測定した結果を述べる. なお,全て Flat MPI で実行し, マクロ部を担当するプロセス数は 128 で固定した. Fig. 8 は全体の並列性能(ストロングスケーリング)を示している. ミ
クロモデル自由度を 3000, 6000 の 2 種類とし, マクロ部とミクロ部の実行形態を 3 種類とした. 実

Fig. 8 並列性能(ストロングスケーリング)

行形態の default はマクロ部とミクロ部を担当す5

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

るプロセス(128 プロセス)とミクロ部のみを担当するプロセスで構成した場合であり, Fig.7 の同時実行は行えないケースである. split と overlap はマクロ部のみを担当するプロセス(128 プロセス)とミクロ部のみを担当するプロセスで構成したときに, それぞれ同時実行有りと無しに対応したケースを表す. 7152 コアにおける並列化効率は3000 自由度で default:38.7%, split:39.4%,overlap:43.0%, 6000 自由度で overlap:96.0%であった. 6000 自由度においては良好な結果が得られた. また 3000 自由度の default と overlap を比較

Fig. 9 3000 自由度(overlap)での全体とミクロ部計算

すると 11.1%程度の性能向上が見られた. 次に時間

3000 自由度についてであるが, 各コア数での全体梗塞を模擬したシミュレーションを行った. 解析の計算時間とミクロ部の計算時間を Fig. 9 に示す.

結果を Fig. 10 に示す. (1）左心室圧と容積の履歴

これによれば, ミクロ部のスケーラビリティは良では, 狭窄による圧および拍出量の低下が現れ始好だが, プロセス数を固定としたマクロ部の計算

めるまでおおよそ 12x11 拍程度に相当する時間が量が占める割合が大きく, そのため全体のスケー経過していることがわかる. 最終的には拍出量がラビリティが悪化したと考えられる. さらなるス

狭窄発生前の 72ml から 47ml まで落ちている. しケーラビリティ向上には, マクロ部の高速化およかも弛緩末期の容積が正常時の 68%程度まで小さび並列化効率の改善が必要である.

くなっていることがわかる. (2）ATP とリン酸 Pi

3.4. 酸素輸送モデルおよび代謝モデルの導入

濃度の全心室壁における平均値に関しては, まず心臓は休むことなく全身に血液を送り続けていPi の増加が先に顕著になり, ATP の減少がその後る. この仕事のための膨大なエネルギーを調達すに続いている. このケースでは Pi の増加はさほどるため心臓はその重量が体重の 1%以下であるにも拍動に深刻な影響を及ぼさず, ATP の減少がエネ関わらず心拍出量の 4-5%の血液を冠循環に取り入ルギー補給とクロスブリッジの解放を妨げることれ, 酸素や栄養の供給をうけている. この供給のにより心拍動のパフォーマンスを低下させている.

不足は現在の大きな健康問題である狭心症, 心筋実際(3) ATP 消費エネルギーと外部仕事に見られ梗塞などの虚血性心疾患の原因となるが, もう一るように, １拍動トータルでの ATP 消費量はほとつの大きな問題である心不全についてもエネルギ

んど変化しないものの, 収縮期での ATP 消費率がーの枯渇が関わっているとの指摘もあり広く心疾下がっていることがわかる. 弛緩期にはこの消費患の病態の解明にはエネルギー代謝を再現したシ不足を補うようかのように ATP 消費が持続されミュレーションは必須である.ているように見えるが, これは梗塞部においてク今回, 毛細管を含む血流解析が可能となったこ

ロスブリッジの強い結合状態からの解放が弛緩期とにより, 血中の酸素輸送モデル, ATP を産生すにも持続することにより ATP 消費が続くからでる代謝モデル 5), ATP を消費して心筋の収縮力をある. １サイクルトータルでの効率は, ATP のト求める興奮収縮連関モデルを連成させることで,ータルの消費量があまり変化しないにもかかわら冠循環を流れる血流から心筋代謝を経て心臓の血ず外部仕事が減少しているため, 効率が狭窄発生液拍出に至るシミュレーションが可能となった.

前の 22%から 13%まで落ちている. (4) 冠状動脈

そこで左冠動脈の LCX(左回旋枝)と呼ばれる部への血液流入量と冠状静脈からの流出量に関して位の血管径を 1%程度に狭窄させることで, 急性のは, 狭窄発生直後の流量減少に加え, 心拍動のパ
6

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月“Approximation for Cooperative Interactions ofa Spatially-Detailed Cardiac Sarcomere Model,”Cellular and Molecular Bioengineering, vol.5(1),pp.113-126, 2011
(2) 国際会議プロシーディングス(3) 国際会議発表
T. Washio, J. Okada, S. Sugiura and T. Hisada,“Numerical modeling and its solution
techniques for multiphysic and multiscale heart
simulator: UT-heart,” Numerical LinearAlgebra – Algorithms, Applictations, andFig. 10 冠動脈狭窄発生からの諸量の経時変化(横軸のTraining. An NOW-JSPS joint Seminar, Delft

目盛ｘ11 が実際の心臓拍動数を表す)

University of Technology, Delft, The

フォーマンス低下に伴う減少もみられる. これは,Netherlands, Apr.10-13, 2012

収縮時の大動脈圧低下により冠状動脈流入口での(4) 国内会議発表

圧力が低下するからである. この流量減少は, 梗鷲尾巧, 岡田純一, 杉浦清了, 久田俊明, タンパ

塞部位以外での酸素供給量の低下をもたらし, 心ク質レベルの現象と心臓の拍動を結びつけるマル室壁全体に悪影響を及ぼす結果となった.

チスケール解析について, 第 61 回理論応用力学講演会(NCTAM2012), 2012 年 3 月 7 日～9 日(5) その他（特許，プレス発表，著書等）(2) 当初計画の達成状況について

特許 1 件出願済み

マルチスケール冠循環シミュレーションにおけるマクロ部・ミクロ部の同時実行については実

参考文献

装・評価を終えた. 冠循環シミュレーションと酸

1) G. S. Kassab, C. A. Rider, N. J. Tang and Y. B.

素輸送・代謝モデルとの連成については予定通り

Fung, “Morphometry of pig coronary artery終了した. スレッド並列化については, マクロ部

trees,” Am. J. Physiol., vol.265, pp.H350-365,の線形ソルバに現れる Gram-Schmidt 直交化に対1993して施した. 一方, ミクロ部計算に対するハイブ

2) G. S. Kassab and Y. B. Fung, “Topology andリッド並列化は, 対称性モデルの導入の効果もあdimensions of pig coronary capillary network,”

り 200 自由度程度の問題サイズで十分な解析が可Am. J. Physiol., vol.267, pp.H319-325, 1994

能であったため, その必要性はないと判断した.

3) G. S. Kassab, D. H. Lin and Y. B. Fung,“Morphometry of pig coronary venous system,”Am. J. Physiol., vol.267, pp.H2100-2133, 1994

4. 今後の展望
本年度の研究で開発した拍動, 冠循環, 代謝を4) G. S. Kassab, J. Berkley and Y. C. Fung,

連成させたシミュレータをさらに高度化し医学的“Analysis of pig’s coronary arterial blood flow

な知見獲得を目指す. また同時に性能面でもさらwith detailed anatomical data,” Ann. Biomed.に向上させる.

Eng. vol.25, pp.204-217, 1997
5) Cortassa S., Aon M. A., O'Rourke B., Jacques

5. 研究成果リスト

R., Tseng H-J., Marban E. and Winslow R. L., “A(1) 学術論文

computational model integrating

T. Washio, J. Okada, S. Sugiura and T. Hisada,

electrophysiology, contraction and7

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

mitochondrial bioenergetics in the ventricular
myocyte,” Biophys. J., vol.91, pp. 1564–1589,
20066) D. Algranati, G. S. Kassab and Y. Lanir,“Mechanisms of myocardiumcoronary vessel interaction,” Am. J. Physiol.
Heart Circ. Physiol., vol.298, pp.H861-H873,
20108