予測医療に向けた階層統合シミュレーション

東京大学
機械工学専攻&大学院工学系研究科
バイオエンジニアリング専攻

高木周

2

階層統合生体力学シミュレーションとその医療応用

シミュレーションによる病態予測と治療支援
方針： ・ミクロ（分子生物学，細胞動態）とマクロ（病態予測，治療法検討）を結びつける．
・マクロ(病態）からミクロ（分子）へトップダウンのアプローチ．
（病態をシミュレートしながら，ミクロ因子の影響を調べる）

階層統合シミュレーションによる新しい予測医療の構築
分子細胞

組織

血小板

血小板

血小板vWFvWFvWF

GP1ｂ

コラーゲンの露出した血管壁

臓器
心筋梗塞

For the Design and Developmentof High Intensity Focused UltrasoundTherapy (HIFU) Device

High Intensity Focused Ultrasound therapy•Prostate cancer•Breast cancer

http://www.prostatecancercentre.co.uk/treatments/hifu.html

HIFU therapy has been developed for the treatment of deeply-placed cancer.

Brain cancer

Liver cancer

http://www.imasonic.com/

Displacement of focal point due to thereflection and refraction of ultrasound atthe interfaces of bones

Control of the focal pointby an array transducer

Recent Development of HIFU Device
Products(1) InSightec(GE) : Exablate 4000
Brain Treatment, Phase Control by Array Transducer(1024 elements),
MR-Guided, (MR-ARFI under development)(1) InSightec: Exablate 4000

(2) Phillips: Treatment of uterus myoma,
MR-Guided(3) Siemens: collaborating withChongqing(重慶) HIFU, MR-Guided(2)Phillips;TemperatureMonitoring by MRI(4) Supersonic Imaging: collaboratingwith Prof. Fink (Inserm), MR-ARFI(4) 256 Channels--------------------------------------------------------------(5) University of Michigan, Profs, C. Cain, B. Fowlkes, Z. Xu
Histotripsy ( Cavitation induced Treatment)(5) Plaque:Utilization of Cavitation, to Control the Focused Area
blown away

HIFU simulator

Comparison of the focal point
Without phase delay

Pmax=1.7MPa
02
Pressure [MPa]With phase delay

Pmax=2.2MPa
02
Pressure [MPa]Experimental Apparatus
StageStage controller
40 mm

Trigger

56ch PZTTransducer

zyxoLabVIEW™
Trigger

61ch Amplifier

Needle hydrophone

Water
Oscilloscope

Acrylic plate (3 mm)
c=2660m/s

Ultrasound through acrylic plate (20deg)
Acrylic plate is arranged with the angle of 20 degrees.
60.3

50.25

4w/o phase shift

y [mm]

30.2

2
10.15

00.1

-1-2-3-15

0.05-12.5-10-7.5-5x [mm]

-2.5

01

60.3

50.25

4w/ phase shift

y [mm]

30.2

2
10.15

00.1-1-2-3-150.05-12.5-10-7.5-5x [mm]-2.501HIFU Ｓｉｍｕｌａｔｉｏｎ for Focus Control for Breast Cancer using MRI data
MRI image

Numerical model

FatWater•••

56ch transducer
Focal distance 100 mm
Frequency 2 MH••••

3D orthogonal mesh
Size 120x110x110 mm
1200x1100x1100 grid
8x4x4 sub domains

Parenchymal
Speed of[m/s]Density[kg/m3]

Fat

1465[1]985[1]

lacteal gland

1547[2]1032[2]connective
tissue1615[1]1090[1]Sound[1] ICU REPORT 61, [2] T. D. Mast,, 2000

w/o control

w/ control

Parkinson’s Disease Simulation
Close collaboration among teams of
K. Doya(OIST),Y. Nakamura(Univ. Tokyo),
T. Nomura(Osaka Univ.),
and S. Takagi (Univ. Tokyo)
12

研究実施体制(筋骨格‐脳神経系 階層統合シミュレーション)
パーキンソン病の階層統合シミュレーションによる病態理解と治療応答予測および運動と脳の関係のモデル化
大脳基底核-皮質-小脳神経回路大規模シミュレーション（OIST－銅谷）
脊髄反射の神経回路モデル（東大-中村，高木）
筋繊維レベルからの筋収縮の有限要素解析（東大―高木）

サルの実験データに基づくパーキン病脳モデル(ミクロモデル）（OIST―銅谷，生理研－橘）患者データに基づくパーキン病モデル（マクロモデル）（阪大―野村刀根山病院ー佐古田）

全身の分布質量筋モデルによる姿勢制御予測（東大―中村）

13

ミクロスケール実験データ：
パーキンソン病モデルサルの脳シグナル (Tachibana et aｌ. (2011)）

マクロスケール実験データ：
パーキンソン病における姿勢反射障害発生機序の仮説モデル (野村（阪大））‐ PLoS One (2009)‐ Hum Mov Sci (2008)‐ J Theor Biol (2012)‐ Math Biosci (2013)

従来仮説
連続的制御モデル

前後(mm)

前後(mm)

左右(mm)

左右(mm)

間欠制御の崩壊？

PSD (rad2/Hz)（シミュレーション）
適切なタイミングで制御OFF
柔かな運動

UPDRS=53(OFF)PSD (rad2/Hz)

新しい仮説
間欠制御モデルパーキンソン病患者

若年健常者

Freq (Hz)Freq (Hz)（シミュレーション）
常に制御ON
高ゲインで安定化
硬い運動今までの実施状況1 (脳神経系シミュレーション) (2011-2014年度)

Diesmann (ユーリッヒ研究所）, 五十嵐（OIST銅谷チーム）ら
2013年7月，ISLiMで開発された脳神経系 シミュレータ NESTを用いて， 「京」上で， 10兆個の結合を持つ大脳皮質局所神経回路のシミュレーションに成功－世界最大の脳神経シミュレーション－
・ 17億3,000万個の神経細胞・ 10兆4,000億個のシナプス結合

マーモセットなどの小型霊長類の全脳規模に匹敵ヒトの脳機能解明に向けた第一歩
プレスリリース http://www.riken.jp/pr/topics/2013/20130802_2/
51616

今までの実施状況2 (2011-2014年度)
筋骨格-神経系階層統合シミュレーション （中村，高木，銅谷）
脊髄神経モデルと筋骨格ワイヤモデルの接続(2) K-Body+HI-Muscle

STIMULATION(1) PyNEST + MUSIC + K-Body

全身骨格モデルと筋有限要素モデルの接続

外力
(刺激)

伸張反射
シミュレーション

実験データを用いた筋・神経のモデルパラメータの同定

Excitatory Synapses
Inhibitory Synapses

下肢の筋活動度推定計算

二頭筋SN発火→IN発火→三頭筋の筋活動抑制

ニューロン番号(サイズ順)

IN
三頭筋SN
二頭筋SN
三頭筋MN

二頭筋MN

運動(MN)・感覚(SN)・介在ニューロン(IN)のスパイク列(サイズの小さいものから順に発火する)

Numerical Method
Suitable for Medical Image Data
(Full Eulerian FSI simulations)Gaehtgens et al.(1980) Blood Cells., 6, 799.

Background
Fluid-Structure coupling analysis of living bodyDiagnostic image (CT, MRI)

Voxel data (Volume Fraction of Constituents)
representing multi-component geometry

without Mesh Generation(Finite-Difference or Finite-Element) Simulationon

Eulerian frame

http://www.rehab.research.va.gov/
jour/02/39/3/ledouxf06.jpg

Full Eulerian approach for Fluid-Structure Interactions
vs.
LagrangianEulerian
How is the two-phase distinguished ?

s1
SolidSolid0.5
FluidFluid0by the boundary of meshby solid volume fraction
How is the solid deformation described?

by the displacement ofmaterial points themselvesby left Cauchy-Green
deformation tensor

・ Sugiyama, Ii et al. (2011) J. Comput . Phys., 230, 596.
・ Ii, Sugiyama et al. (2011) Int. J. Numer. Meth. Fluids, 65, 150.

Basic equations for conservation and constitutive law
  v  0,  m   t v  ( v   ) v    p    σ m ,Cauchy 's stress tensor in mixture formσ m  (1   s )σ f   s σ s ,solid volume fraction
fluid (Newtonian):solid (visco-hyperelasto):σ f  2  f D ,

Mooney-Rivlin model:

1D  v  vT  ,2strain rate

σ s  2 c1B   2 c2 (tr(Β )Β  B  B )  2  s D .

left Cauchy-Green deformation
 xi
B  F  F , Fij ,X j
T・ Sugiyama, Ii et al. (2011) J. Comput . Phys., 230, 596.

Basic equations for kinematics of solid motion
solid volume fraction

 t s  ( v   ) s  0,
left Cauchy-Green deformation tensor

 t B  (v   ) B  L  B  B  L ,TwhereL  v .
T・ Sugiyama, Ii et al. (2011) J. Comput . Phys., 230, 596.

Validation I

Comparison with available numerical data

Gao & Hu (2009)
J. Comput. Phys. 228, 2132.
particle-particle interactions
in a shear flow
full LagrangianParticle-particle interactions in a shear flow
● material points (tracers)(roll-over and bounce-back modes)Movies=f=1, Lx=8, Ly=4, =20, C1=0, C0=80

contour: vorticity (spin)

Vtop=1-3

Fluid
SolidSolid

Vbottom=11

Comparison with Gao & Hu（JCP2009)'s results

Gao & Hu (full Lagrangian) Present (full Eulerian)
NxxNy=1024x512

t=0.2

t=0.8

t=4.0

t=5.6

t=30.0t=33.2Extension to 10-particle system● material points (tracers)Movies=f=1, Lx=8, Ly=4, =20, C1=0, C0=80

contour: vorticity (spin)

Vtop=1-3

Vbottom=11

Feature: full Eulerian approach・makes it easily possible to perform FSI simulations
with complex geometry
Initial voxel data

FSI simulation

without mesh generation/reconstruction

・ Sugiyama, Ii et al. (2010) Comput . Mech., 46, 147.
・ Nagano, Sugiyama et al. (2010) J. Fluid Sci. Tech., 5, 475.
・ Takagi, Sugiyama et al. (2012) J. Appl. Mech., 79, 010911.

The Super-Fast FSI Solver: ZZ-EFSI

ZZ-EFSI achieved the actual
speed of 4.5 PETA FLOPS!!!
Software availablle at
http://www.islim.org/islim-dl_e.html

Development of Multiscale
Thrombosis Simulator

Multiscale modeling of initial stage of thromobosis

Blood flow(continuum mechanics)

Protein(wWF)-Protein(GP1b) binding(stochastic process)•Substance diffusion•Metabolic reaction and activation•Morphology change

Molecular interaction(molecular dynamics)

Blood flow simulation• [0,44]×[0,22]×[0,22] [m], D=20 [m],• 30 RBCs, 10 platelets (Ht ≈ 20 [%])
• <u> ≈ 2 [mm/s]

320x160x160 mesh, MPC-RICC at RIKEN (256 cores)

・ Ii, Sugiyama et al. (2012) J. Biomech. Sci. Eng.,7, 72.

Stochastic Monte Carlo Simulation of
Cell adhesion on vessel wall
R+LkfkrC

R： Receptor
L： Ligand
C： R-L complex
kf： Binding React. Rate Cnst.
kr： Dissosc. React. Rate Cnst.

Receptor－Ligand Bind. Forcef b  Function of  l  l0 

Jianrong Li, et. al(2009)

PlateletVessel Wallll0

Ligand-Receptor Binding Model
Stochastic model with energetic elasticity
Eyring (1935) J. Chem. Phys., 3, 107.
Bell (1978) Science, 200, 618.
Dembo (1988) Proc. R. Soc. Lond. B, 234, 55.
Hammer & Apte (1992) Biophys. J., 63, 35.

P f  1  exp(  k f  t )  R f  formation
Pr  1  exp(  k r  t )  Rr  breakage
Luo et al. (2007) Blood, 109, 603.

forward reaction ratef   p ( l  l0 )Formation

( l  l0 ) 2 k f ( l )  k f 0 exp   ts

2 k bT 

reverse reaction rate

( R f , Rr  [0,1] : Random numbers )Breakage( l  l0 ) 2 k r ( l )  k r 0 ex p  ( p   ts )
2 k bT 

Model Parameter p  10  4 [N/m],  ts  0.9 p [N/m]l0  60 [nm], k r 0  3 [s  1 ]Fox et al. (1988) J. Biol. Chem., 263, 4882.
Arya et al. (2005) Biophys. J., 88, 4391.
Kim et al. (2010) Nature, 466, 992.

Motion of an adhering platelet
dimensionless time t* =  t
 = 1000 (1/s)

 = 3000 (1/s)

 = 6000 (1/s)

Simulated Example : Wall-Adhering Platelet moving in Shear FlowExperiment
Reininger et al., blood, 107(2006) 3537

Simulation(= 3000 s-1)

シミュレーションモデル構築のための実験(後藤信哉(東海大）より提供）
Ht value© Goto(Tokai Univ.)

Number of adhering platelets

シミュレーション用データ取得のためのフローチャンバーを用いた実験

血小板粘着に対する赤血球の影響に関する実験データ（後藤（東海大）ら 2012)

Experimental data by Goto & Tamura.

Adhesion of Platelets in a Shear Flow

Ht = 5.5 % = 800 (s-1) t  200 (ms)

Ht = 15.3 %Ht = 21.9 %Numerical Simulations for the Experiment by Goto et al.
Ht = 0%Ht = 20%

Domain Size: 144mx72mx72m# of grid points: 480x240x24

本年度：実験系と同じ流路サイズのシミュレーション
Ht = 20%

Comp. extent:
400mx100mx100m
Num. grid points: 2,048x512x512
2,048 nodes (16,384 cores) on the K computer

injuredwall
platelet adhesion41Effect of arteriosclerosis size ( Hematocrit Value: 20% )Effect of arteriosclerosis size ( Hematocrit Value: 0% )

抗血小板薬のモデリング
循環器系
・ 抗血小板薬の薬効再現のためのモデリングとシミュレーション

45

実測データを用いた循環器系統合シミュレーション
開発済み

開発中
局所３次元計算

術後予側シミュレーション
Left external illiac artery (Preoperation)
Left posterior tibial artery (Preoperation)
Left external illiac artery (Postoperation)
Left posterior tibial artery (Postoperation)

Blood pressure(mmHg)

140

0D-1D-3D連成モデル

120

神経系モデルの導入

100

8060
45Time(s)Femoral artery (preoperation)
Femoral artery (Postoperation)
Graft25Blood flow(ml/s)

201510504-55Time(s)患者データの活用
およびデータ同化

0D-1D全身血管シミュレーション

血圧計開発

マルチスケール血栓症シミュレーション

課題３．予測医療に向けた階層統合シミュレーション
まとめ
現状・医療応用に向けてソフトウェアの準備ができ，本格計算が始まったところ・ 真の意味での医療応用はこれからが勝負．
医学・医療さらには看護に対して何を創出できるか．

今後の予定，
・超音波診断・治療器設計に向けた計算・血栓シミュレータによる血小板粘着過程の再現と抗血小板薬の薬効評価（P2Y12阻害薬他）
・脳神経系と筋骨格系の統合シミュレーション（パーキン
ソン病の再現）
シミュレーションが支援する新しい医学・医療の創出ヘ…..

Research Achievement (2015)

Won a Distinguished Simulation Award!!!

UT-Heart （Hisada et al.)

謝辞
Principal Investigators
後藤信哉（東海大・医），姫野龍太郎（理研），横田秀夫（理研），松本洋一郎（東大・工），久田俊明（東大・新領域），大島まり（東大・情報学環），劉浩（千葉大・工），天野晃（京大・情），和田成生（阪大・基礎工），小野謙二（東大・生研），松澤照男（北陸先端大・情報セ），山口隆美（東北大・工），野村泰伸（阪大・基礎工），
岡澤重信（広島大・工），中村仁彦（東大・情報理工）
銅谷賢治（沖縄科技大）, Huaxiong Huang (York Univ.)，
長谷部光泉（東海大・医）Main contributors for software development):
杉山和靖（東大），伊井仁志（東大），沖田浩平（日大），塩崎聖治（理研），野田茂穂（理研），石川顕一（東大），梁夫友（理研），山村直人（理研），片岡博之（理研），石峯康浩(理研），K.Lee(理研），X. Gong（上海交通大）, 清水和弥（東大），島本憲夫（東大），世良俊博（阪大），古田琢哉(理研）