HPCI計算科学フォーラム
2019年10月9日 東京大学

生体力学シミュレーションによる個別化医療支援～生体現象の理解とリアリティーの追求～

大阪大学 大学院基礎工学研究科
機能創成専攻 生体工学領域
大阪大学国際医工情報センター（兼）和田 成生shigeo@me.es.osaka-u.ac.jp

力学理論を通じた生体現象の理解とその応用診断手法の確立
治療法の策定薬効評価
医療機器の開発

細胞

生体分子

組織

微小血管

マルチスケール解析

細胞膜のＭＤ

血栓症
癌転移

ドラッグデリバリー

計算力学解析
（コンピュターシミュレーション） 末梢循環不全

人体

褥瘡
大動脈瘤
心不全

血管病進行予測喘息
肺癌

大規模並列計算によるミクロ－マクロ流動現象の統合

心筋梗塞
脳動脈瘤

動脈硬化

複合臓器
臓器統合モデルによる病態評価と予測

数値臓器－医用画像統合による診断支援

臓器大血管マルチフィジックス現象理解のための連成シミュレーション
微小循環における細胞流動

Flow

血流と血球細胞の変形との連成

Takeishi et al.,Phys. Rev. E., 2015

毛細血管網における酸素輸送血流と細胞変形と物質輸送の連成

Kobayashi et al.,APS-DFD, 20183

マクローミクロをつなぐ
マルチスケールシミュレーション

せん断流れ場の血球流動

血液の見かけ粘度

Takeishi et al., JFM 872, 818-848, 2019.
4

非平衡分子動力学シミュレーション
Strain of the bilayerthickness

ε z = Lth / Lth0 − 1
Lth :Bilayer thickness

Lth 0 :Initial Bilayer thickness(~4 nm)5 nm

Koshiyama K & Wada S, J. Biomech. 44, 2053-2058, 2011

RuptureAreal straine A = A A0 - 1
A :Area of the simulation box
A0 :Initial area of the simulation box(~38.46 nm2)

膜破断の分子メカニズムの観察
コレステロール無 コレステロール有(40 mol%)

コレステロール含有量の影響

実験結果
二重膜構造から指組膜構造(Lβ-I相）への相転移

Needham & Nunn,
1990

Shigematsu T, et al., Chem. Phys. Lipids, 2014 2 nm
6

流体音響シミュレーション
歯茎摩擦音の再現/u/ /s/ /i/

筋収縮による舌の能動変形シミュレーション

Yoshinaga et al., 2017, J. Accoustic. Soc. Am.

医学・生物学・物理を結ぶ計算力学シミュレーション
人体・臓器レベル

組織・細胞レベル

生体分子レベル

イメージベースト
シミュレーション

マルチスケール／
マルチフィジックス
シミュレーション

分子動力学
シミュレーション

生理機能

生物学的応答

生化学反応

臨床医学

細胞生物学

分子生物学

カラーＭモード・ドップラ画像

OSAKA
UNIVERSITY

プローブ

心尖

心尖

正常例

拡張不全例
速度

血流伝播速度

H

血流L
僧帽弁

大動脈弁

0僧帽弁
超音波ビーム

start

時間
拡張期

end時間

OSAKA

左心室拡張早期における心室内流れ

UNIVERSITYt = 0.10 sec

t = 0.12 sec

t = 0.14 sec

t = 0.16 sec

000079 158

95 19095 19092 184(cm/sec)
t = 0.18 sec

t = 0.20 sec

t = 0.22 sec

t = 0.24 sec

000087 174

77 15471 14266 133

左心室拡張早期における心室内流れOSAKAUNIVERSITYt = 0.13 sec

t = 0.14 sec

000MV

AVzyAWt = 0.12 sec

xPW
18 36

29 58

29 58(cm/sec)
t = 0.15 sec

t = 0.16 sec

t = 0.17 sec

t = 0.18 sec

000036 72

42 84

49 9854 108

血流伝播速度の物理的意義
シミュレーション

僧帽弁大動脈弁

超音波ビーム臨床

velocity(cm/sec)
189

10.0
9.0

170

8.0

151.1

7.0
6.0120

5.090

4.0
3.060

2.030

1.00001.2 2.4
0.6 1.8Time, t (×10-1 sec)

Distance from the mitral orifice

心尖Distance from the entrance of the LV, x (cm)

プローブ

Time

Nakamura, et al., Biomech Model in Mechanobiol, 2003.

リアリティーの追求と物理的理解
リアリステイックな現象の再現• 形のリアリティー• 物理的リアリティー

等価なシンプルなモデルでの解釈

理屈の理解脳動脈瘤コイル塞栓術
マイクロカテーテルの挿入

コイルでのフレーミング

完成図

100 µm
石口, INNERVISION,12(11), 47-53,1997

14コイル充填による滞留効果
瘤内部の流線図
q = 0°瘤内血流の滞留領域体積の評価|v|<0.0001 m/s →滞留領域と定義

q = 45°
Stagnant Volume Ratio (%)

1000.000

0.052

Velocity magnitude (m/s)
0.156
0.104

0.208

0.260

80

現行の基準

0°
9°
18°
27°
36°
45°
54°

60
40200010203040
Packing density (%)

5060

臨床において用いられる基準である20%のコイルを留置した場合でも、瘤内へ直接血流の主流が流入する形状の場合、血流の滞留が不十分な可能性を指摘
Otani et al., Medical & Biological Engineering & Computing 51, 901-910 , 2013

バーチャルコイルモデル
• 伸縮
• 曲げ
• ねじり
• 接触qi
Li ri

ポテンシャルエネルギで表現

コイル
動脈瘤

コイル節点の運動方程式を逐次解法で解析

16

コイルの計算力学モデルsEuler-Bernoulli梁としてモデル化
全長L の3次元上で滑らかな曲線
曲線座標 S動標構（moving frame）

弾性エネルギ
𝑈=ℝ2𝑥𝑧

1% 𝑘' 𝜀 ) + 𝑘+ 𝜅- − 0 𝜅- ) + 𝑘1 𝜅) − 0 𝜅) ) + 𝑘1 𝜅2 − 0 𝜅2 ) 𝑑𝑠
2 6

引張圧縮
𝑦

ねじり

曲げ𝜀: 軸⽅向ひずみ，𝜅- : 梁の軸⽅向回転率，𝜅) , 𝜅2 : 梁の曲率𝜕𝑈 𝜕𝑈 𝜕𝑈 𝜕𝑈
梁の内力 𝐟 = n, 𝐦 = 𝜕𝜀 , 𝜕𝜅 , 𝜕𝜅 , 𝜕𝜅)
2Otani et al., Med. Biol. Eng. Compt., 2017

接触の表現例

・コイル – コイル
・コイル – 血管壁

𝒇+𝒇B

法線方向コイル

線形粘弾性モデルを用いたペナルティ法

𝒇B = −𝑓B 𝒏 − 𝜂B 𝑣B 𝒏接線方向Slip-stickを考慮した一般化Coulomb則に従う摩擦力
・Slip friction
𝒇+ = −𝜇'HIJ 𝑓B𝒗+𝒗+ + 𝜀

n : 法線ベクトル
𝜂 : 粘性係数𝑣 : 相対速度𝑓B : 反発⼒𝜇'HIJ : 動摩擦係数
𝜇'+IKL : 静摩擦係数
𝛿+ : 接線⽅向移動量
kt : 弾性係数（Margolis, 2005）
・Stick friction 𝒇+ ≤ 𝜇'+IKL 𝑓B

𝒇+ = −𝑘+ 𝛿+𝒗+− 𝜂 + 𝒗+𝒗+ + 𝜀

塞栓コイル内血流シミュレーション
流体解析
n NS方程式
n 連続の式

Otani et al., Otani et al., J. Biomech. Eng., 201819

多孔質体モデルとの比較
Coil model

速度Velocity
magnitude (m/s)
0.70
0.60
0.50
0.40
0.30
0.20
0.100.0

せん断速度

Shear rate (s-1)
100＝8060

周辺との速度差40→ 血栓生成に関与

200.0

Porous media model(Darcy’s law + Ergaun’s eq.(1952))

Packing density 27%30%

コイル塞栓術治療支援
CT/MRI
仮想コイルモデル

試行実験

形状抽出

• 手技の検討• コイルの選択• 最適なコイルの設計

個別化医療支援

PC-MRI（血流計測）

ＣＦＤ解析

血流速度

• 治療効果の評価• 予後予測21

医用画像に基づく左心房内血流動態の推定
左心系

正常な左心系の運動（CT画像から三次元構築）

肺静脈・心房細動 (Wolf et al., 1991)
・病的リモデリング (Piccini et al., 2011)
・肺切除 (Yamamoto et al, 2016)

血流の滞留
左心房

左心耳
血栓形成が頻発

左心室血栓形成脳梗塞 (Wolf et al., Stroke, 1991)
認知症 (Thacker et al., Neurology, 2013)22

4D-CT動画に基づく左心房内の血流動態の再現
4D-CT images at 20 phases(phases 0%RR, 5%RR,…, 95%RR)・・・(a) 旋回流の形成(c) 旋回流の再形成(d) 逆流の発⽣(b) 旋回流の消失

1. 基準となる左房形状を作成(b)

2. 表面形状抽出(a)・・・

物質輸送解析
3. 最尤推定法（EMアルゴリズム）によ→ ⾎流の滞留状態を評価(c)(d)る点群マッチング

Myronenko & Song, IEEE Trans. Patt.
Anal. Machine Intell, 2011

計測結果(Fyrenius et al., 2001)と定性的に一致
Otani T, et al, Ann. Biomed. Eng. 44, 3284–3294, 201623

医用画像に基づく心原性脳塞栓予測
数百例の左心房シネ画像
ジョンズホプキンス大 足利先生より提供

シミュレーション

データ科学

モデル駆動
アプローチ

データ駆動
アプローチ

＋ 流体解析 → 統計分析
形態指標に基づく病態評価

正常心房肥大症患者

物理指標に基づく病態評価24

臨床データを理解するアプローチ
理論的アプローチ

現象論的アプローチ

シミュレーション

実現象

理論モデル（先見的）計測・観察計算条件（仮定）
数値解析
物理現象としての再現• シミュレーションの”ちから”• リアリティーの追求

ビッグデータ

因子抽出（発見的）
支配因子の同定
因果関係（しくみ）

相関分析・統計解析 データ科学

計測 VS シミュレーション
PC-MRI

CFD analysis

どちらを信じますか？

データ同化血流シミュレーション
計測データ

シミュレーション

従来の血流解析

全速度データの誤差の総和を最小化

データ同化のための手法
仮想流体力を導入したフィードバック制御

G wall(Funamato & Hayase, 2005-2013)Ñ×v = 0Ginr (¶t v + v × Ñv) = -Ñp + µÑ v + f2

仮想流体力G out

G wallf = - Kv ( v - v m ): measurement velocity vm
P controller (P制御)

Control variable(操作量，入力)

Process variable(制御量，出力)

Setpoint(目標値)

Proportional gain(比例ゲイン)

fv

vm

Kv

非物理的流体現象

圧力場の補正 ( Ii & Wada, 2014): measurement velocity vm

G wall

物体力を含む流れの支配方程式r (¶t v + v × Ñv) = µÑ 2 v - Ñp + f
f = - Kv ( v - v m )Gin

フィードバックさせる流体力は，圧力勾配によってのみ生じると仮定

G out

G wall

1. Body force is represented by using a scalar potential
f = -Ñf + Ñ ´ ψ ¬ Ñ 2f = -Ñ × f where ¶ nf = 0 at G in , G out , G wall

2. New pressure field is determined by using the scalar potential as theboundary condition-Ñp + f = -Ñ( p + f ) = -ÑP where P = f at G in , G out

時系列PC-MRI速度を用いたデータ同化結果
ü PC-MRI解像度: 0.7 mm
ü CFD解像度: 0.2 mm
ü 周期 T: 1.47 sec
ü レイノルズ数: 100~400
ü ウオマスリ数: ~3
平均速度変化の比較

データ同化による推定結果* Ii et al., 2017, Int J Numer Meth Biomed Eng, e2910.

誤差（最大流速時）

動脈瘤部

主血管部

Hisham MA, et al, MBEC 55, 1605–1619, 2017.

ポスト京重点課題２

個別化・予防医療を支援する統合計算生命科学
課題責任者：宮野 悟（東京大学医科学研究所）
サブ課題A
大量シーケンスによるがんの個性と時間的・空間的多様性・起源の解明
サブ課題責任者 宮野 悟（東京大学医科学研究所・教授）
サブ課題実施 小川誠司（京都大学大学院医学研究科・教授）
サブ課題B
データ同化生体シミュレーションによる個別化医療支援
サブ課題責任者 和田成生（大阪大学大学院基礎工学研究科・教授）
サブ課題実施者 高木 周（東京大学大学院工学系研究科・教授）〃
松田秀雄（大阪大学大学院情報科学研究科・教授）〃
野村泰伸（大阪大学大学院基礎工学研究科・教授）〃
玉川裕夫（大阪大学歯学部附属病院・准教授）
サブ課題C
心臓シミュレーションと分子シミュレーションの融合による基礎医学と臨床医学の架橋
サブ課題責任者 久田俊明（ （株）UT-Heart研究所 ・ 東京大学名誉教授）
サブ課題実施者 永井良三（自治医科大学・学長）
32

全脳循環シミュレーター開発

全脳脳血管モデル全脳血流シミュレーション

脳機能画像

全脳血管モデル脳血流と脳機能との関係解明
n 均質な血流が維持されるメカニズムn 脳血管障害による影響

脳機能障害のメカニズムの解明
n 血流と脳機能との関連性
n 脳機能画像としてどのように見えるのか33

全脳動静脈血管モデルの構築
動脈静脈

１．医用画像を用いた数mm径までの主要血管の構築

２．手順1の末端から脳表面へ向けた数百𝜇𝑚径程度までの脳全体を覆う血管の構築

３．手順2の末端から脳表面に沿って広がる，数十𝜇𝑚径程度の血管の構築

患者個別医用画像データ ＋ 領域制約型マルチレベル血管分岐アルゴリズム
Kitade et al., 2019, CMBE2019 Proc, 1:859-862.

34

主要脳動脈の支配領域
面積比 ACA（青）:MCA（緑）:PCA（赤） = 1.60 : 2.55 : 1

Top view

Bottom view

Front view

Back view
35

脳動静脈血管網モデル

血管密度:8.50 (/𝑚𝑚) )先行研究のデータ
貫通細動脈の直径 (𝜇𝑚)15-240貫通細静脈の直径 (𝜇𝑚)20-125貫通血管の密度 (本/𝑚𝑚) )8.44-15Wahbi K, et al. J. Microcirculation 23, 580-590, 2016

Duvernoy HM, et al, Brain Research
Bulletin 7, 519-579，1981
36

ボクセル型数値流体解析
領域相を表す相関数

中間領域 Gfor x Î W F
ì1ï
f (x) = í(0,1) for x Î G
ï0for x Î W S
î

固体方程式（e.g. 剛体，変形体）
流体方程式（e.g. 非圧縮NS）

相関数fを用いて平均化方程式を導出
f×流体方程式＋(1-f)×固体方程式

e.g. 剛体連成（Weymouth, Yue 2011）
超弾性体連成（Sugiyama et al. 2011）37数値流体解析の超大規模化
n 適合サブドメイン法

形状に沿ったボクセル生成
不要な格子を削減
リソースの有効活用
ノード間負荷バランスの改善
n 格子ボルツマン法

完全陽解法スキーム
局所性の高いアルゴリズム
超大規模計算に適した手法最小血管径：400μm
空間解像度：60 μm
総格子数: 4.5億（内，流体4億)
計算時間：20,000ノード×24時間
実時間：6秒間
3839

個別データに基づく３D脳血流解析
被験者61⼈のMRA⾎管形状データ(Wright SN et al., Neuroimage 82, 170-181, 2013)40

ウィリスの動脈輪付近の速度場
Sample 2

Sample 6

PCA付近の流れの構造に個人差を確認
41

血圧分布
正常脳モデル

脳梗塞モデル（ LMCAの⾎管径を80%減少）

42

まとめ

力学理論に基づく計算シミュレーションは，臨床データから物理的要因に起因する特徴（決定論的要素）を抽出するのに有効な手段である．
計算シミュレーションに立脚した理論的アプローチと，データ科学による現象論的アプローチを融合することで，個別化医療支援に供するより高度な情報を提供できる．

謝 辞
大阪大学

武石 直樹

助教

大阪大学

横山 直人

特任准教授

大阪大学

大谷 智仁助教

大阪大学

野崎 一徳講師

大阪大学

重松 大輝助教

大阪大学 基礎工学研究科
バイオメカニクス研究室

滋賀医科大学渡邊 嘉之 教授
ｼﾞｮﾝｽﾞﾎﾌﾟｷﾝｽ大学 足利 洋志 助教
豊橋科学技術大学 吉永 司

助教

首都大学東京

伊井 仁志

准教授

徳島大学

越山 顕一朗 准教授

神戸大学

今井 陽介

教授

神戸大学

石田 駿介

助教

https://sites.google.com/site/biomechwadalab/