平成 29 年度

文部科学省

ポスト「京」で重点的に取り組むべき社会的・科学的課題に関するアプリケーション開発・研究開発

平成 29 年度
「個別化・予防医療を支援する統合計算生命科学」
成果報告書

平成 30 年 5 月 31 日
国立大学法人東京大学
宮野 悟

本報告書は、文部科学省の科学技術試験研究委託事業による委託業務として、国立大学法人東京大学が実施した平成29年度「個別化・予防医療を支援する統合計算生命科学」の成果を取りまとめたものです。

１．委託業務の題目 ...................................................................................................................................1
２．実施機関（代表機関） ........................................................................................................................1
３．委託業務の目的 ...................................................................................................................................1
４．平成 29 年度（報告年度）の実施内容 .................................................................................................2
４－１．実施計画 ...................................................................................................................................2
４－２．実施内容（成果） .....................................................................................................................5
1.

大量シーケンスによるがんの個性と時間的・空間的多様性・起源の解明（サブ課題 A） ..........52.

データ同化生体シミュレーションによる個別化医療支援（サブ課題 B） .................................35

3.

心臓シミュレーションと分子シミュレーションの融合による基礎医学と臨床医学の架橋（サブ課題 C） ............................................................................................................................................83
4.

プロジェクトの総合的推進 ...................................................................................................... 113

４－３．活動（研究会・受賞・書籍等） ...........................................................................................120
４－４．実施体制 ...............................................................................................................................121
別添 1 学会等発表実績
別添 2 実施計画ｉ１．委託業務の題目「個別化・予防医療を支援する統合計算生命科学」

２．実施機関（代表機関）
機関名

国立大学法人東京大学

所在地〒113-8654 東京都文京区本郷七丁目３番１号

ふ り が な
課題
代

責任者

表

氏

みやの

名

生年

西暦 1954 年 12 月 5 日（63 歳）

悟

月日※2018 年 4 月 1 日現在

宮野

所属部署名
連絡先

さとる

医科学研究所
Tel. 03-5449-5615

役職

教授

Fax. 03-5449-5442

E-mail: icls-office@hgc.jp, miyano@ims.u-tokyo.ac.jp

機

さいとう

ふ り が な

関
事務
連絡

氏

名

所属部署名

担当者
連絡先

齋藤
医科学研究所

としや
利弥

研究支援課

外部資金戦略チーム
Tel. 03-5449-5137
E-mail:

役職

主任

Fax. 03-6409-2017

t-gshikin@ims.u-tokyo.ac.jp

３．委託業務の目的
病気は、臓器群の変調という現象として現れるが、その背景には生命の設計図とも呼ばれるゲノムがあり、オミクスと呼ばれるエピゲノム、RNA、タンパク質など多彩な分子が細胞を制御・構成している。
また、細胞には環境や加齢により長い時間をかけて変化していく個々人で異なる細胞コンテクストがあり、そのもとで構成されている臓器の状態も多様である。さらに、その繋がりは人智を超えた複雑さを有している。その理解には、画像や生理データなどを含む高精度臨床データとともに環境・生体・時空間的にゲノムから全身を捉える必要がある。現在、ゲノム解析技術の劇的な革新と高精度計測機器の急速な発展は、健康・医療ビッグデータを生み出そうとしており、この傾向は急激に加速している。こうした背景のもと、本研究の目的は、ポスト「京」によって初めて実現できる「情報の技術」と「物理の原理」の融合により、がんをはじめとして全脳・循環・代謝系など、全身の疾患に対して、ビッグデー
タを活用し、高度の生体階層統合シミュレーションに個体データを同化させる技術、並びに、ライフサ
イエンスにおいてかつてない大規模なデータ解析技術を開発・応用することにより、また、大規模デー
タに基づくアプローチと並行して、分子細胞レベルの研究と臓器個体レベルの研究を融合させ、ミクロとマクロのメカニクスとを関連させて定量的に捉えたシミュレーションモデルを構築することにより、病態の理解と効果的な治療法の探索を行い、その成果を個別化・予防医療へ返す基盤となる統合計算生
命科学を確立することを目的とする。
このため、国立大学法人東京大学を中核機関として、分担機関である国立大学法人京都大学、国立大学法人大阪大学、株式会社 UT-Heart 研究所、学校法人自治医科大学と密接に連携し、研究開発を実施する。
1

４．平成 29 年度（報告年度）の実施内容
４－１．実施計画
平成 29 年度は、本格実施フェーズの 2 年目であり、平成 26～27 年度に実施した調査研究・準備研究
および平成 28 年度の本格実施フェーズ開発結果に基づき、以下に示す重点課題 2 に関するサブ課題（A～C）（1.～3.）について以下の目標を定め、研究開発を実施する。また、プロジェクトの統合的推進 4.
を行う。
1.

大量シーケンスによるがんの個性と時間的・空間的多様性・起源の解明（サブ課題A）
がんはゲノムの変異によって生ずる疾患で、個体別、空間的（腫瘍の場所）、時間的多様性を持っている。がんの起源とその多様性の理解は、がん治療戦略、がんの予防法と超早期発見にイノベーションを起こし、副作用に優しく個人ごとに効き目のよい薬を創出するための戦略上大変重要である。そのため、次の目標を設定し研究を実施する。
(1) 世界の追随を許さない大規模がんオミクスデータ解析システムの開発
平成 29 年度は、大規模がんオミクスデータ解析システムの継続発展（高速化・高感度化）を図り、Genomon に基づくポスト「京」の詳細設計コデザインを実施しつつ、高感度解析において「京」で100 検体/日のパフォーマンスを達成することを目標にして、各機能を確認してシステムを発展させることを実施内容とする。ここで高感度解析とは、通常の SNV 検出や融合遺伝子検出に加え複雑な変異解析(中間的範囲の欠失挿入やウイルス挿入検出などの構造変異検出等)を含む。
(2) がんの起源と集団内多様性を理解するためのシステムの開発
平成 29 年度は、平成 28 年度に獲得した基礎データと構築した基礎となるモデルに基づいたシステ
ムを開発することを目標にし、少なくとも一つのがん種について大量シーケンスデータなどを利用
して、がんに移行しやすい変異パターンの解析を実施する。
(3) 胚細胞バリアントの効率的な解析システムの開発
平成 29 年度は、平成 28 年度に検討したゲノムシーケンスデータに基づいた高感度な胚細胞バリア
ント検出法の検討に基づき、胚細胞バリアントの効率的な解析が行えるシステムを開発することを目標とする。そのために、これまで検討・開発してきた低頻度体細胞変異検出の情報解析技術などをシステム開発に用いる。検出法の開発を継続するとともに、各種がんの大量のエクソームシーケ
ンスデータを用いて正常対照に見られる変異を解析し、がん腫による変異プロファイルの違いを比較する。
2.

データ同化生体シミュレーションによる個別化医療支援（サブ課題B）
ポスト「京」を活用した大規模生体物理シミュレーションと生体計測データとを様々なレベルで同化・融合させることにより、実測データを重視する医療に受け入れられる計算機シミュレータを開発し、シミュレーションで得られる物理情報を活用した個別化医療支援の実
現 を 目 指 す 。 そ の た め 、 平 成 29年 度 は 次 の 目 標 を 設 定 し 研 究 を 実 施 す る 。

2

(1) 個別化医療支援のための生体シミュレータの開発および評価
平成 28 年度に構築した計算モデルに基づいて、個別化医療支援に供する生体シミュレータを開発し、「京」等の大規模計算環境でその性能を評価する。具体的には、全脳循環機能評価のための階層統合シミュレータ、脳血管系における血栓形成シミュレータ、ヒト運動機能解析のための神経筋骨格系統合シミュレータ、構音機能分析のための流体音響連成シミュレータ、微小循環遊走細胞の動態解析のための大規模画像処理ソフトウェアを開発する。
(2) データ同化手法の開発と個別化医療支援への応用
平成 28 年度に開発したデータ同化手法を用いて、脳循環系における血流計測データと血流シミュレーションとを同化させ、脳動脈瘤のコイル塞栓治療の個別化医療支援に応用する。
3.

心臓シミュレーションと分子シミュレーションの融合による基礎医学と臨床医学の架橋（サブ課題C）
タンパク機能（遺伝子）の変化が最終的に臓器レベルにどのような変化を及ぼすかを物理的に解明・予測することができれば、その応用は無限に広がる。ポスト「京」のパワーを活かし、本サブ課題では、世界で初めて心臓シミュレーションと分子シミュレーションを融合させた真のマルチスケールシミュレーションを実現する。これにより計算科学の歴史に新たなマイルストーンを築くと共に、現実の医療・創薬にブレークスルーをもたらす。
(1) 心臓病、特に「心不全」の根本的解明と最適治療を可能とするマルチスケール心臓シミュレータの開発を目的とする。モデル化は心不全研究において得られた遺伝子レベルの変化、発現調節、細胞内情報ネットワークに関する膨大な情報、臨床情報を効果的に取り入れて推進する。
平成 29 年度は、UT-Heart と CafeMol の試験的連成計算を実現する。具体的には「京」で計算可能な程度の小型化モデル、例えばリモデリングを用いて連成計算を実施し、連成手法の検討を行う。
また心不全治療と外科手術の慢性期予測手法の研究に着手する。
(2) 目標 1 は分子レベルの力学と心臓のポンプ機能の関係に焦点を当てたものであるが、同様の技術を
致死性不整脈の病態解明と予防・治療に対して適用する。即ち、重点課題 1 において進められるイ
オンチャンネルの分子シミュレーションと UT-Heart を融合することで、候補物質の分子構造と遺伝子情報のみから不整脈のリスク評価を可能とするシステムの完成を目標とする。
平成 29 年度は、前年度の開発を継続し、最終的にはカリウムイオンチャネルについて、従来の確率モデルを分子シミュレーションから得られた確率モデルに置き換え、UT-Heart による興奮伝播シミュレーションを実現する。
4.

プロジェクトの総合的推進
プロジェクト全体の連携を密としつつ円滑に運営していくため、運営委員会や研究会の開催等、参画各機関の連携・調整にあたる。
特に、プロジェクト全体の進捗状況を確認しつつ計画の合理化を検討し、必要に応じて調査あるいは外部有識者を招聘して意見を聞くなど、プロジェクトの推進に資する。またニュ3

ースレター（電子版）発行、ホームページ運営、シンポジウム開催、人材育成活動などを通じて研究の進捗と連携を推進する。
プロジェクトで得られた成果については、積極的に公表し、今後の展開に資する。

4

４－２．実施内容（成果）
課題全体としての進捗は順調であり、一部、予定よりも進んでいるものがある。ポスト「京」の課題におけるサブ課題 A の超えねばならない目標は大規模シーケンスデータ解析パイプライン Genomon による多数検体のゲノムデータ解析とその高速化に対応することである。そのためには Virtual Grid
Engine が不可欠であることの結論に至った。その開発の中では「京」自身の問題がいくつか判明したがそれを回避する道を見出し、また高速化のための様々な技法を駆使した結果、目標を達成したことを実験結果より検証し、来年度へつなげることができた。また、Genomon に基づくコデザインにおいては、
いくつかの重要な貢献があったと考えている。同時に、開発の途中の段階であるが、がんゲノミクスにおける世界の注目を集める研究成果が得られ始めた。これはサブ課題 A が優れた戦略性を有していることを示している。サブ課題 B は中間評価の結果を受け、構成を全脳循環シミュレータの開発と循環器系
個別化医療支援に向けたデータ同化シミュレータの開発の二つに簡明化することで、データ同化生体シ
ミュレーションによる個別化医療支援のために必要な様々な詳細アプリケーションの位置づけを明確にした。それぞれのアプリケーションの開発の進捗は順調であった。サブ課題内での連携もよい。サブ課
題 C は、心臓シミュレータ UT-Heart と分子シミュレータ CafeMol の融合により、ミクロ・マクロ間の相互作用により病態が進行する心不全の解明と最適治療を可能とするマルチスケール心臓シミュレーシ
ョンを実現することが最大の目標である。分子シミュレーションと心臓シミュレーションの連成法には、計算時間の観点から分子モデルの軽量化が望まれる。そこで首振りメカニズムのエッセンスを損なうことなく自由度が削減できるかどうか検討する為に、ミオシン構造変化の機構分析を行い、キーとなる残基群を調査して成果を得た。そしてサブ課題 A と連携し（心筋の RNA 解析、ゲノム解析等）、また心不
全治療と外科手術の慢性期予測手法の研究に着手した。さらに、もう一つの目標である心毒性スクリー
ニングシステムの開発では、重点課題 1 の寺田 透博士との連携により、心臓シミュレータ UT-Heart と心筋細胞イオンチャネルの分子シミュレーションを融合し、創薬における候補物質の分子構造と遺伝子情報のみから不整脈のリスク評価を可能とするシステムの開発を進めた。新たなチャネル開閉動作の確率モデルを考案し、チャネル分子モデルの特性を端的に反映した連成解析シミュレーションが実現した。
1.

大量シーケンスによるがんの個性と時間的・空間的多様性・起源の解明（サブ課題 A）(1) 世界の追随を許さない大規模がんオミクスデータ解析システムの開発
A) Virtual Grid Engine の開発(a) 概要重点課題 2 のターゲットアプリケーションである Genomon2 の「京」/ポスト「京」上での動作のために、Virtual Grid Engine (VGE)の開発を実施した。VGE は非グリッドエンジン環境において、MPI のプロセス空間を疑似的なグリッドエンジン環境としてユーザアプリケーションに提供するミドルウェアである。
ゲノムシーケンスや RNA シーケンスデータの解析ではパイプラインが一般的であり、本体および内部利用の外部ソフトウェア等の更新が頻繁（~10 回/年）なバイオインフォマティクス分野において、ソフトウェアの個別 MPI 化は発見に至るまでの時間的競争という生命科学的優位性の観点から現実的な対応ではない。グリッドエンジンはポスト「京」へ搭載しないことが決定
しており、「京」を含めた HPCI 計算機群や他のペタフロップス級スーパーコンピュータ上での5柔軟な運用のためには、必要不可欠なソフトウェアである。平成 28 年度に実施した VGE の実装と移植した Genomon に基づき、本年度は実データを用いたパイプライン全体の実行と多数検体解析における高速化を実施した。次節以降で本年度実施内容を順次述べる。
(b) ロバスト性改善「京」割当資源が少なかった頃、東京大学医科学研究所のスーパーコンピュータ「Shirokane」
を様々な面で活用した。その際、正常に計算が実施できないケースを確認した。VGE は、前提とするソフトウェアやライブラリが揃っている以外の環境限定因子は存在せず、Shirokane シ
ステムでの動作不良が設計上は起きないはずである。そのため、原因の究明から本年度の作業を開始した。
VGE は MPI ライブラリを、MPI4PY モジュールを介して利用している。Shirokane には MPI
環境として MPICH、MPICH2、OpenMPI の三つがインストールされている。MPI を利用する場合には、該当プログラム/ソフトウェアをビルドする際に利用した MPI 環境と並列実行の際に使用する MPI 環境をそろえなければならない。さらに、ビルドに際しては選択した MPI ラ
イブラリがバックボーンとして使用する C/Fortran コンパイラも揃っている必要がある。
Shirokane システムにインストール済みのライブラリ/パッケージではこれら前提条件が完全には判別できないものがあるため、今回は MPI ライブラリを含むパッケージ/ライブラリをビルドし直した。用意した環境は OpenMPI（ver.1.10.4）および MPICH（ver.3.2）の二つである。
この二つの MPI 環境を用いてそれぞれ MPI4PY のビルドを実施したところ、問題なくインストールが完了した。
上記環境を用いて、Genomon2 を VGE 環境下で実行スクリプト（図 1.1）を使用したテストを実施した。MPI には OpenMPI を選択し、MPI4PY モジュールも同ライブラリを用いてビルドしたものを指定した。このテストにおいてエラーが確認された。エラーメッセージを図 1.2に示す。次に、MPICH 環境（MPI4PY も対応したものを選択）で同じテストを実行したところ、テストは正常に終了し Genomon は正常に解析結果を出力した。この結果から、VGE 実行
時の問題は VGE そのものではなく、MPI4PY を介して利用する MPI ライブラリに問題があると判断した。
#!/bin/sh –x#$ -S /bin/sh #set shell in UGE#$ -l s_vmem=16G,mem_req=16G#$ -pe mpi 16#$ -cwd

source /etc/bashrcexport HOME=/home/fj_hase#MPIHOME=/home/fj_hase/.local/openmpi_1.10.4
export PATH=$MPIHOME/bin:$PATH
export LD_LIBRARY_PATH=$MPIHOME/lib:${LD_LIBRARY_PATH}
export LD_LIBRARY_PATH=$MPIHOME/lib/openmpi:${LD_LIBRARY_PATH}

6

export LD_RUN_PATH=$MPIHOME/lib:${LD_RUN_PATH}#
export PATH=$HOME/public/python_2.7.11/bin:$PATH
export PYTHONHOME=$HOME/public/python_2.7.11
export LD_LIBRARY_PATH=$PYTHONHOME/lib:${LD_LIBRARY_PATH}
export PATH=.:$HOME/public/bin:$PATH
export LD_LIBRARY_PATH=${HOME}/public/lib:${LD_LIBRARY_PATH}
export LD_LIBRARY_PATH=${HOME}/public/lib64:${LD_LIBRARY_PATH}#
mpirun-n $NSLOTS -machinefile$TMPDIR/machines vge &vge_connect –start
genomon_pipelinedna./sample/sample.csv ./output ./genomon.cfg./dna_task_param.cfg
vge_connect --stop

図 1.1: Shirokane 用 VGE & Genomon 実行用スクリプト
An MPI process has executed an operation involving a call to the "fork()" system call to create a childprocess. Open MPI is currently operating in a condition that could result in memory corruption or othersystem errors; your MPI job may hang, crash, or produce silent data corruption. The use of fork() (orsystem() or other calls that create child processes) is strongly discouraged.

The process that invoked fork was:

Local host: sc336 (PID 6895)
MPI_COMM_WORLD rank: 0

If you are *absolutely sure* that your application will successfully and correctly survive a call to fork(), you
may disable this warning by setting the mpi_warn_on_fork MCA parameter to 0.

図 1.2: テスト実行エラーメッセージ（OpenMPI 環境）
さらなる原因究明のため VGE およびシステムのエラーメッセージを調べたところ、VGE がsegmentation fault で異常終了していることが確認された。図 1.2 のエラーメッセージによると、
OpenMPI ではシステム関数の一つである fork()を用いた子プロセスの生成に問題があるため、その利用を強く抑制していることが分かった。本問題について、OpenMPI 公式サイト FAQ[1]に以下の説明を発見した。

38. Can I use system(), popen(), or fork() in an MPI application that uses the OpenFabrics support?

7

The answer is, unfortunately, complicated. Be sure to also see this FAQ entry as well.

If you have a Linux kernel before version 2.6.16: no. Some distros may provide patches for older
versions (e.g, RHEL4 may someday receive a hotfix).

If you have a version of OFED before v1.2: sort of. Specifically, newer kernels with OFED 1.0 and
OFED 1.1 may generally allow the use of system() and/or the use of fork() as long as the parent does nothing
until the child exits.

If you have a Linux kernel >= v2.6.16 and OFED >= v1.2 and Open MPI >= v1.2.1: yes. Open MPI
v1.2.1 added two MCA values relevant to arbitrary fork() support in Open MPI:
btl_openib_have_fork_support: This is a "read-only" MCA value, meaning that users cannot change
it in the normal ways that MCA parameter values are set. It can be queried via the ompi_info command; itwill have a value of 1 if this installation of Open MPI supports fork(); 0 otherwise.
btl_openib_want_fork_support: This MCA parameter can be used to request conditional, absolute,or no fork() support. The following values are supported:
Negative values: try to enable fork support, but continue even if it is not available.
Zero: Do not try to enable fork support.
Positive values: Try to enable fork support and fail if it is not available.

Hence, you can reliably query OMPI to see if it has support for fork() and force OMPI to abort if you request
fork support and it doesn't have it.
（中略）This is unfortunately no way around this issue; it was intentionally designed into the
OpenFabrics software stack. Please complain to the OpenFabrics Alliance that they should
really fix this problem!
原因の根本は OpenMPI がバックボーンで利用する OpenFabrics ドライバーに存在することが判明した。よって、OpenMPI 環境を利用する際には同通信ドライバーを使用しないようにする
設定（オプション-mca btl ^openib）が必要である。本設定により、OpenMPI 環境においても正常に実行がなされるようになった。
(c) VGE 高速化
H28 年度では、VGE の実行は小規模テストまでであった。本年度は多数検体処理の確認を実施するため、VGE の性能分析を実施した。具体的には、VGE の利用ノード数の増加と VGE へ投入するジョブ数の増加により、VGE のジョブ処理能力がどの程度劣化するかを確認した。
テストコードは 120 秒スリープのみを実行内容に含む genomon_test である。検証条件は四8

つであり、検証 1 ではノード数 200・ジョブ数 2000、検証 2 ではノード数 200・ジョブ数 15000、
検証 3 ではノード数 2000・ジョブ数 2000、検証 4 ではノード数 2000・ジョブ数 15000 である。
「京」で genomon_test を実行すると、必要資源量は 120 秒×ジョブ数である。利用可能計算
資源（ワーカー数）が同時処理可能であるため、上記資源量をワーカ数で除算した値が計算時間の理想値となる。検証 1 の理想実行時間は 1320 秒であり、実際の実行時間は 1358.421 秒であった。検証２の理想実行時間は 9120 秒であり、実際の実行時間は 9246.825 秒であった。検証 3 の理想実行時間は 240 秒であり、実際の実行時間は 251.015 であった。検証 4 では理想の実行時間は 960 秒であり、実際の実行時間は 1117.022 秒であった。
表 1.1: 検証条件

注）ノード数（MPI 並列数）に対し、ワーカはマスターランク一つ分少ない数になる。
図 1.3 は検証 1（199 ワーカ・120 秒スリープ×2000 ジョブ）における各ジョブの開始および終了時刻をプロットしたものである。VGE へ投入されたジョブ数がワーカ数より多いため、各
ワーカは 10 ないし 11 のジョブを処理することになる。
ジョブの第 1 ターン（1 番目～199 番目のジョブ）、第 2 ターン（200 番目～398 番目のジョブ）の結果をそれぞれ図 1.4、図 1.5 に示す。図 1.3 より、2000 ジョブのジョブ開始時間の推移
はほぼ階段状になっていた。図 1.4 より、ジョブ実行順は、立ち上がり部とジョブプール時で安定にジョブが投入され線形に増加する部分（ジョブプール部）の二つに分かれていた。第 1 ターンの立ち上がりは 164 ミリ秒/ジョブで、ジョブプール時は 9.2 ミリ秒/ジョブとなり、1 ジョ
ブを処理する平均トータル時間は、17.8 ミリ秒/ジョブとなった。第 2 ターンは、第 1 ターンのジョブ投入残がなく VGE のジョブ投入処理待ちがなければ、理想的には第 1 ターンと同じようなジョブ投入パターンになるはずである。図 1.5 の第 2 ターンでは、第 1 ターンと同じような

9

傾向になっており、立ち上がりは 164 ミリ秒/ジョブ、ジョブプール時は 11.0 ミリ秒/ジョブで、
1 ジョブを処理する平均トータル時間は 19.4 ミリ秒/ジョブとなった。
図 1.6 は検証 2（199 ワーカ・120 秒スリープｘ15000 ジョブ）のジョブの第 1 ターン（1 番
目～199 番目のジョブ）のジョブ投入時間である。検証 1 と比べてジョブ数を 5 倍にした場合（検証 2）でも、第 1 ターンのジョブ投入傾向は検証 1 と同じようになり、ジョブ投入の立ち上
がりと、ジョブプール部の二つに分かれた。検証 1 と違い、検証 2 での立ち上がりは 842 ミリ
秒/ジョブで、ジョブプール部は 52.8 ミリ秒/ジョブとなり、ジョブ投入性能が低下した。
図 1.7 は検証 3（1999 ワーカ・120 秒スリープ×2000 ジョブ）のジョブの開始時刻と終了時
刻をプロットしたものである。検証 1 との違いは、ワーカ数が約 10 倍になっていることである。
検証 3 のジョブ投入傾向は検証 1 と同じで、ジョブ投入の立ち上がり部と安定したジョブプー
ル部に分かれており、立ち上がり部は 197 ミリ秒/ジョブ、ジョブプール部は 17.6 ミリ秒/ジョ
ブで、1 ジョブを処理する平均トータル時間は 18.5 ミリ秒/ジョブであった。検証 1 と 3 より、
ワーカ数が大きくなってもジョブ投入傾向は、ほぼ同じであり、ジョブの立ち上がり部とジョブプール部に分かれており、1 ジョブの投入時間もほぼ同じであった。
図 1.8 は検証 4（1999 ワーカ・120 秒スリープ×15000 ジョブ）のジョブの開始時刻の第 1
ターンの結果である。ジョブ投入は、ジョブ立ち上がりとジョブプール時の二つに分かれており、それぞれのジョブ投入時間は 930 ミリ秒/ジョブ、63.1 ミリ秒/ジョブであった。検証 2 と同じように、ジョブ数が同じでワーカ数を増やした場合には、ワーカ数増加の影響は少なく、1ジョブ投入時間は同じ傾向であることがわかった。以上の結果を表 1.2 にまとめた。

図 1.3: 検証 1（199 ワーカ・120 秒スリープ×2000 ジョブ）における各ジョブの開始時刻と終了時刻の実行結果

10図 1.4: 検証 1（199 ワーカ・120 秒スリープ×2000 ジョブ）のジョブの開始時刻と終了時刻の実行結果の第 1 ターン（1－199 番目）の結果

図 1.5: 検証 1（199 ワーカ・120 秒スリープ×2000 ジョブ）のジョブの開始時刻と終了時刻の実行結果の第 2 ターン（200－398 番目）の結果

11

図 1.6: 検証 2（199 ワーカ・120 秒スリープ×15000 ジョブ）のジョブの開始時刻と終了時刻の実行結果の第 1 ターン（1－199 番目）の結果

図 1.7: 検証 3（1999 ワーカ・120 秒スリープ×2000 ジョブ）のジョブの開始時刻と終了時刻の実行結果の第 1 ターン（1－1999 番目）の結果

12

図 1.8: 検証 4（1999 ワーカ・120 秒スリープ×15000 ジョブ）のジョブの開始時刻と終了時刻の実行結果の第 1 ターン（1－1999 番目）の結果
表 1.2: 各検証における第 1 ターンのジョブ開始までの VGE オーバーヘッド

13図 1.9: VGE ジョブコントローラの処理フローチャート
VGE のマスタープロセスが実行している、MPI ジョブコントローラの処理フローを図 1.9 に示す。現状では、ワーカの受信前にループ構造があり、ジョブをチェックし（check new job）、
新たなジョブがあれば、ジョブプールへジョブデータを保持する。そして、ジョブを作成し（make new job）、ワーカからの呼び出し信号をチェックし、ワーカからの呼び出しがあればワーカへジョブ送信＆ジョブ結果データの受信を行い、続いてワーカからのジョブ結果データの処理を行う。もし、呼び出しがなければジョブチェックへ戻るループ構造になっている。
このように、現状では新ジョブのチェックを優先する設計となっている。これは、VGE 設計
指針として、VGE と Genomon が同じマスターノードで動く仕様であることから、なるべくマ
スターノードで動作する VGE の負荷を抑えるようにしたためである。
VGE の各処理部にタイマーを入れた場合の性能計測結果を表 1.3 に示す。実行条件は「京」の会話型ジョブで、0 秒スリープ×2000 ジョブである。0 秒スリープジョブのため、計測結果
はほぼ VGE のオーバーヘッドといえる。全体処理時間は 53.44 秒であり、タイマー計測結果より、ワーカ受信前処理では、check new job が 3.88 秒（コスト 7.3％）、make new job が 3.64秒（コスト 6.8%）、check worker が 27.66 秒（コスト 51.76%）であった。ワーカに対するジョブ受送信は 1.18 秒であり、コストは 2.2%であった。ジョブ送信後の処理時間は 8.78 秒でコ

14

ストは 16.4%であった。
このことから、VGE ジョブコントローラにおけるボトルネックは、check
worker 部にあることがわかった。
VGE は Python 言語で記述されている。そのため、「京」と Intel 環境で同じジョブを実行し、「京」と Intel 環境の Python 性能を比較した。図 1.10 は「京」での VGE 性能測定結果であり、計測条件は 40 ワーカ（41 MPI プロセス）で、genomon_test の 120 秒スリープ×2000 ジ
ョブ結果の第 1 ターン（40 ジョブ）のジョブ実行開始時刻のプロットである。ジョブ投入の立
ち上がりは 156.5 ミリ秒/ジョブ、ジョブプール時は 6.8 ミリ秒/ジョブで、1 ジョブを処理する
平均トータル時間は 51.7 ミリ秒/ジョブであった。Intel 環境（東大医科研スパコン Shirokane）
での VGE（バージョン r20170317d、以降 asis 版と表記）の性能結果を図 1.11 に示す。ジョブ投入の立ち上がりは 25.8 ミリ秒/ジョブ、ジョブプール時は 1.1 ミリ秒/ジョブで、1 ジョブを処理する平均トータル時間は 7.9 ミリ秒/ジョブであった。これらの結果をまとめたものを表 1.4に示す。この結果より、「京」と Intel 環境での VGE のジョブ処理の時間は、約 6 倍の差があり、Intel 環境のほうが速いことがわかった。これらの結果は Python と MPI 通信性能によるものであるが、主に Python の処理性能に起因すると想定される。Python の処理能力に「京」とIntel 環境では大きな差があることがわかった。ちなみに Python 版数は、「京」では Pythonver2.7.13 および mpi4py2.0.0（計算ノードインストール済みの環境）であるのに対し、Intel
環境では Python ver2.7.10 および mpi4py2.0.0 であった。
表 1.3: VGE ジョブコントローラ各部の処理時間計測結果15図 1.10: 「京」での VGE オーバーヘッド測定結果

図 1.11: Intel 環境（東大医科研スパコン Shirokane）での VGE オーバーヘッド測定結果16表 1.4: 「京」と Intel 環境（東大医科研スパコン Shirokane）での VGE オーバーヘッド測定結果比較

以上の結果に基づき、VGE の高速化改良を実施した。大規模ジョブ処理の対応で一番大きな問題となっているのは、
ジョブ送信前処理である。改良後の送信前処理フローを図 1.12 に示す。
大規模ジョブ時にジョブの立ち上がり時間が遅い要因はジョブプールサイズが非常に大きくなり、ジョブ情報検索に時間を要することが原因である。VGE のジョブプールの大きさは、高々ワーカ数＋α程度あれば十分である。そこで、ジョブプールサイズを制限する設計が、ジョブ処理高速化に有効である。さらに VGE は Python 共有オブジェクトを用いて、VGE マスターの Genomon ジョブコントローラと MPI ジョブコントローラとの情報共有をしているが、このPython 共有オブジェクトはローカルオブジェクトよりも（共有オブジェクトという制約上）ア
クセス性能が悪い。そのため、MPI ジョブコントローラ側にローカルの配列で 2 次ジョブプー
ルを用意し、MPI ジョブマスターが送信準備・情報探索のためにこの 2 次ジョブプールだけにアクセスするようにすれば、共有オブジェクトへのアクセスが低減され高速化が可能である。
さらに、ジョブ送信後の処理・ワーカ受信前処理で共有オブジェクトのアクセス低減を図る改良を加えた。

17図 1.12: 改良後の VGE ジョブコントローラ処理フロー「京」で改良後の VGE の性能測定を実施した。テストは 47 ワーカ（48 MPI プロセス）、テ
ストプログラムには genomon_test のスリープ時間を 0 秒に変更したものを用いた。実行条件
としてジョブ数を 10,000、20,000、50,000、100,000 とした。計測結果を表 1.5 に示す。10,000
ジョブでの VGE のオーバーヘッド 55.6 秒に対して、20,000 ジョブでは 106.0 秒（10,000 ジョ
ブに対して 1.91 倍）、50,000 ジョブで 267.6 秒(10,000 ジョブに対して 4.81 倍)、100,000 ジョ
ブでは 475.4 秒（10,000 ジョブに対して、8.55 倍）という結果になった。100,000 ジョブの実行に成功し、かつ、オーバーヘッドがジョブ数におよそ比例することが確認できた。
表 1.6 は、「京」での改良後の VGE の実行時間である。計測条件は 1999 ワーカ（2000 MPIプロセス）、genomon_test の 120 秒スリープ、ジョブ数は 100,000 である。この結果、改良後の VGE では 100,000 ジョブの実行が正常に終了し、genomon_test の実行時間は 6752.08 秒でVGE のオーバーヘッドは 632.08 秒となった。この結果から、10 万ジョブでの VGE のジョブ処理性能は平均 6.3 ミリ秒/ジョブとなった。今回、10,000 ジョブ規模でのジョブ処理を確認でき、さらに 100,000 ジョブまで実行できることを確認した。

18

表 1.5: 「京」での改良後 VGE の性能測定結果

5.0

表 1.6: 「京」での VGE（asis 版および改良後）の性能測定結果

19B) コデザイン(a) 概要
本年度も昨年度に引き続き、高速化取り組みと問題点解決のための作業を実施した。具体的な検討項目としては、Genomon パイプライン後半部分のボトルネックである mutation_call の高速化チューニング、Genomon パイプラインの各ステップにおけるメモリ、ファイルサイズ使用量の確認等である。メモリやファイルサイズ等に関しては NDA 関連情報を含むため、本報告では触れることが出来ない。最後に、大規模解析の試行とそれに伴うファイルシステム関連の問題について報告する。
(b) mutation_call の高速化チューニング
mutation_call は Genomon パイプラインの後半部分に相当し、前半のアラインメント処理を経たデータに対して変異を検出する。本部分は計算機資源消費量の観点からは前半部分（特にBWA によるアラインメント部分）と比較して 30 並列程度とごく小さな割合である。しかしながら、解析に要する時間は 10 時間を超え、スループットの観点からは非常に大きな比重を占める部分となっている。
mutation_call 部分の並列数は Genomon 側で指定されており、可変ではない。これは、デー
タの分割区間の前後で検出への影響が入らないよう、生物学的に問題がない区間を選択した結果である。今回、この区間に関して開発者と協議した結果、当初の 30 区間程度から 275 区間への区間数増加が可能となった。並列数がおよそ 9 倍になることから、計算時間の大幅な縮小を期待したが、実際には計算時間はほとんど変化しなかった。そこで、mutation_call 内で最も計算時間の掛かる fisher comparison 部分（表 1.7）について、その内容を詳細に検討した。
表 1.7: mutation_call 部分の計算時間内訳

20図 1.13: fisher comparison のプロファイリング結果
図 1.13 に fisher comparison 部分のプロファイリング結果を示す。メイン関数の処理時間の大部分が Pileup_out 関数部分に集中（91％）していることが確認された。詳細な性能解析を実施するために、本部分を抽出したテストコードを作成し計算時間を測定した。図 1.14 に測定結果を示す。ボトルネック部分の計算について、テストコードではさらに区間分割により詳細な実行時間を測定可能としている。測定では、デフォルトの区間と 10 分の 1 に再分割した区間の2 パターンでの測定を実施した。その結果、大部分の区間が 100 秒前後（デフォルト区間サイズ）に対し、一つの区間だけが特異的に 10,000 秒以上の計算時間を要していることが分かった。
また、この特異点部分は区間サイズを縮小しても計算時間が縮小しない点も判明し、並列化等の高速化を阻害してしまうことも判明した。

図 1.14: テストコードの区間ごとの計算時間（赤丸は特異点）

21

特異点部分は染色体 16 番の塩基番号 46350000～46459999 部分であった。本部分におけるpileup 計算はおよそ 6 分程度で終了しており、計算時間の大部分は pileup 後のループ部分にあることが分かった。本ループはリード長やリード数ではなく、各塩基位置に何本のリードがあるか（depth）に依存した演算となっている。そこで特異点を含む区間の depth を調べた（図1.15）。depth はシーケンサに投入するサンプルの調整具合により多少のばらつきを持つものであるが、平均値から大きく外れることは稀である。今回使用したデータでは、depth の平均値は15～25 程度であり、実際に全データの平均値はその範囲を逸脱していない。しかし、特異点近傍部分には図 1.15 から分かるように最大で 2 万超、1000 以上の塩基位置が区間の 36%程度も含まれていた。depth 値の異常な値に関しては様々な可能性が考えられるが、いずれの場合においてもこのような超高 depth 値を示す部分に関してはデータの信頼性を担保できないため、最終的な変異検出からは除外するべきであるというのが生物学的な判断となる。そこで、パイプラインに修正を加え、depth 値 1000 をスレッショルドとした pileup およびこれ以降の計算打
ち切り処理を追加した。その結果、最大 22,069 秒かかっていた計算を 2,657 秒と 9 倍程度の高速化を達成した。これらを踏まえ、チューニング後の Genomon の解析性能はおよそ 261 検体/日となった。

図 1.15: 特異点近傍における depth 分布(c) 大規模解析の試行と問題点
VGE および Genomon の高速化を踏まえ、本年度後期に大規模解析のテストを試みた。これまで、VGE および Genomon に関しては最少検体数またはテストコードのみの実施であったため、実データを用いて多数検体の解析を行い、基本性能を測定することが目的であった。

22

本テストの試行過程で判明した問題点は二つである。一つはステージングの限界である。当
課題で設定しているデータサイズはポスト「京」実機使用時で 1TB/検体、「京」使用時で 500GB/検体である。今回使用したデータは準備の都合もあり、1 検体のデータサイズはこれらよりも小さな値である。しかし、全データサイズが 4.2TB のジョブがステージングに過大な時間が掛かってしまうために実行できなかった。我々の目標は 1000 検体/日であり、データ総量は単純に考えて 500TB～1PB である。圧縮等を考慮したとしても、僅か数 TB のシステム内データ転送にトラブルを抱えるようでは問題外である。
次に、分散ファイルシステムへの過負荷問題である。高速化チューニングとして、我々は積極的に並列化、分散処理化を施してきた。これにより計算時間が圧縮される一方で、パイプラ
イン各処理でのデータ授受に使用されるファイル数が激増した。オリジナルの Genomon では分
割数が「京」/ポスト「京」より小さく、また、検体ごとにジョブが独立し実行時刻等も分散する傾向にあることから特に問題が生じなかった。しかし、VGE を用いて多数の検体を同時に解析する本課題の運用では、多数のファイルの出現/消去が同時に発生する確率が高く、このような事象が発生した際には分散ファイルシステムのメタデータサーバ（MDS）への負荷が著しく増大し、場合によってはシステムがダウンする。本問題についての抜本的な対応は、MDS を使用しないローカルディスクの使用であるが、「京」には存在しない。
最後に、このような問題の発生に対し、運用側との直接的な連絡手段の開設を希望する。これらの問題の発生は大規模解析時に起こっており、小規模（2000～3000 ノード）では観測されない。大規模ジョブはジョブ実行までの待機時間がとても長く、また、資源消費も大きいため何度もトライすることが出来ない。実行の不具合確認から、ヘルプデスクを介した運用機関への問い合わせを行ったが、メールのみのやり取りの上、間に必ずヘルプデスクを介在させる関係で十分な情報をフィードバックしてもらうことが出来ず、問題の解決に至らなかった。今回の問題は、最終的にはサブワーキンググループにおけるミーティングで議題に取り上げていた
だき、石川 PL を介して直接運用部門の方から詳細な問題個所の情報を提供していただくことで初めて解決した。このような対応は、他の一般課題では不可能であり、今後の対応方法を改めていただくことを希望する。
参考文献[1] https://www.open-mpi.org/faq/?category=openfabrics#ofa-forkC) ネットワーク解析ツール（SiGN）開発の進捗
遺伝子ネットワーク推定ソフトウェア SiGN はノンパラメトリック非線形回帰によるベイジア
ンネットワークを基盤とした SiGN-BN、時系列データを対象として線形状態空間モデルを用いる SiGN-SSM、および L1 正則化法を用いる SiGN-L1 から構成される。そのうち SiGN-BN は、
ベイジアンネットワークの構造推定が非常に困難な問題であることから、適用するデータのサイズ（遺伝子数）に応じてアルゴリズムを使い分ける必要がある。また、その為アルゴリズム
上の研究開発の余地が多く存在する。これまでの研究開発により小規模な遺伝子ネットワーク
推定においては、数学的な意味において厳密な最適解を求めるアルゴリズムの改良により適用

23

可能な遺伝子数が向上しつつある。特に前年度までにおいて開発された TRPOS-2D と呼ばれる
アルゴリズムは、「京」のトーラス・インターコネクト構造を活用した並列最適解構造推定アルゴリズムであり、離散モデルにおいて 36 遺伝子のネットワーク推定を実現している。本年度は前年度から行っている本アルゴリズムの改良が完成し、新たに ParaOS-DC という名称で、著者のウェブサイト ( http://ytlab.jp/paraos/ ) にて公開した。本アルゴリズムは必要な計算ノード間のみの通信に限定することでメモリ消費量を削減したアルゴリズムで、それにより比較的大規模な計算時など、メモリ量が限られた状態において高い並列化効率を実現している。具体的には、離散モデル 33 変数 864 ノード利用時、旧来のアルゴリズム(TRPOS-2D)において 23,919
秒かかっていた計算が、新規アルゴリズムでは 4,429 秒で完了することを確認している。また「京」を用い 20,736 ノード、165,888 コアの使用により離散モデルにおいて 37 変数の最適解の探索をわずか 2 時間 17 分で行えることを実証した。これは制約のない非巡回有向グラフ探索問題で文献上世界最大の問題サイズである。開発したアルゴリズムとこれらの成果を Journal of
Parallel and Distributed Computing において発表した。遺伝子ネットワーク推定に適用可能な非線形回帰モデルを利用した実装もすでに利用可能な状態にあり、「京」互換機およびヒトゲ
ノム解析センター向けの実装が SiGN-BN のウェブページ ( http://sign.hgc.jp/signbn/ ) で公開中である。
(2) がんの起源と集団内多様性を理解するためのシステムの開発
A) Genomon による大量シーケンスデータの解析によるがんの個性と時間的多様性の解明
多様な疾患を含む症候群である骨髄異形成症候群について、この層別化は予後予測因子としても機能するため、個別の患者のリスクに応じた治療の最適化に役立つものである。骨髄異形成症候群は、発症の初期には慢性の経過を呈する前がん病変の性質を持つものの、そのクローン
進化の結果、二次性の急性骨髄性白血病へ進展を起こす。ひとたび急性白血病を発症した場合、治癒は非常に困難であるため、その原因となるクローン進化を明らかとすることは喫緊の課題である。
本研究では、昨年度に大きな成果を上げた、がんの時間的多様性・起源の全般的な理解(Makishima et al. Nat. Genet 2017) を踏まえ、個別の遺伝子の役割に関して詳細な検討を行った。骨髄異形成症候群、そこから派生した二次性急性骨髄性白血病において、大規模オミクス
データの解析を行い、遺伝子変異、遺伝子発現、染色体のコピー数を統合的に解析し、特定の遺伝子が白血病（がん）への進展に関わることを明らかにした。具体的には、まず骨髄幹細胞分画に特異的な遺伝子発現の異常が遺伝子変異と関連して白血病の発症リスクを高めることを明らかにした。本年度はさらに、骨髄異形成症候群について、これを CD34 陽性細胞の発現プ
ロファイルに基づいて、生物学的に性質の異なる二つの群に層別化できることを示した(Shiozawa et al. BLOOD 2017) （図 1.16）。この層別化は予後予測因子としても機能するため、個別の患者のリスクに応じた個別化医療の最適化に役立つと考えている。さらには、骨髄異形成症候群では初めて、SWI/SNF 複合体遺伝子の変異およびコピー数異常を起こし、特定の遺伝
子異常を伴って白血病進展を来すことを発見した (Sakai et al. Leukemia 2018) 。
以上の結果を活用することにより、慢性期の骨髄異形成症候群が将来急性骨髄性白血病へ進24

展するリスクを評価し、白血病の発症を予測することが可能となる。今後は、これまで蓄積したエビデンスの再現性を確認し、臨床的に有用なパネル検査として活用されるよう、前向き試
験などによる検証を行う。このように、骨髄異形成症候群におけるクローン進化の解明ががんの時間的多様性・起源の理解とともに大きく進展した。

図 1.16: 骨髄異形成症候群の発現プロファイルによる層別化
さらに、進化のプロセスを正確にとらえ、個別化・予防医療へと発展させるためには、次項の B)とも関連するが、単細胞レベルで遺伝子の発現と変異を同時に測定する系の開発が不可欠である。この技術を用いて変異が獲得される順序や細胞系列を調べ、クローンの進展の過程を詳細に捉えることを行おうとしている。図 1.17 の例では、TP53 変異をもつ症例において、5qの欠失は巨核芽球と赤芽球のクラスタにみられるのに対して、7q の欠失は赤芽球クラスタのみにみられることを示している。こうした開発により、単細胞を用いた解析を行うことで、高度なクローン構造の解析を目指している。

図 1.17: Single cell 解析によるクローン進化の解明

25B) がんの進化モデルの研究：異なる進化ダイナミクスを統一的に記述する新しいがんの進化シミュレーションモデル(a) はじめにがんは長い間体細胞変異の獲得および変異を獲得したクローンへの自然選択が段階的に起こることによる、単一の正常な細胞から始まる直線的進化により生じる均一な細胞集団であると捉えられてきた。しかしながら、近年、次世代シーケンサを用いた多領域シーケンスにより進化の過程でクローンが分岐することで多数のサブクローンが生み出されていることが明らかになってきた。この腫瘍内不均一性の形成は自然選択による分岐進化（これ以降単に分岐進化とする）により行われる場合と中立進化により行われる場合があることも明らかになっている。さ
らに近年のゲノム解析によって、染色体レベルの大きな変化によって急速な進化が爆発的に起こる断続的進化が起きていることも観察されている。以上、これまでに、直線的進化、自然選
択による分岐進化、中立進化、断続的進化の４種ががんの主な進化ダイナミクスとして提唱されている（図 1.18）。実際はこれらの進化ダイナミクスは、がん種、がんの進展段階によっても異なると考えられるが、どのような条件の違いが進化ダイナミクスの違いを生み出しているかは不明なままである。

図 1.18: ４種のがんの進化ダイナミクス(b) ドライバー変異のみのモデル
そこで我々はがんの進化シミュレーションを用いてこの問題に取り組んだ。はじめにドライバ
ー変異のみを仮定した以下のようシンプルなエージェントベースモデルを構築した。


変異のない一つの正常細胞から単位時間あたり一定の確率 g で細胞分裂を行う（細胞は不死化していると仮定し死亡確率は 0 とする）。



一回の細胞分裂で生まれた二つの娘細胞は n ~Pois(m/2) 個変異を獲得し（すなわち分裂

一回あたり平均 m 個の変異が生まれる）、各細胞は変異を最大 3 個の変異まで蓄積できる。


一変異あたり増殖速度が f 倍になる。すなわち細胞の分裂確率は g ・fN

26(N = Σn: 細胞の持つ変異数)で表せられる。
このモデルを用いてヒトゲノム解析センターのスーパーコンピュータ上で様々な組み合わせパ
ラメータ値で超並列シミュレーションを行い、各シミュレーションから得た一細胞プロファイ
ルを精査することで以下のような結果を得た（図 1.19）。
f が大きい、すなわちドライバー変異が強い場合、最初の変異が入ったクローンが増殖し始め、そのクローンが次の変異を得てさらに増殖速度を増すステップを繰り返すことで、直線的進化が起こる。

f が小さい、すなわちドライバー変異が弱い場合、最初の変異が入ったクローンが増殖し始めるが、そのクローンが次の変異を得る前に、変異を得ていない細胞に次の変異が入ることで、分岐進化が起こる。

図 1.19 :ドライバー変異のみのモデルを用いた解析の結果(c) 中立変異のみのモデル
また次に中立変異のみを仮定した以下のようなエージェントベースモデル構築し、一回の細胞分裂で生まれた二つの娘細胞は n ~Pois(m/2) 個の変異を獲得し、各細胞は最大 1000 個の変異まで蓄積できる。


各変異は増殖確率に影響は与えない。

このモデルを用いた超並列シミュレーションによるパラメータ依存性解析から、以下のような結果を得た（図 1.20）。


変異率 m が大きい場合、の変異が入ると細胞分裂の系譜がトレースされることよってフラクタル様の中立変異のパターンが現れる。

27

図 1.20: 中立変異のみのモデルを用いた解析の結果(d) ドライバー、中立、両変異をもつ複合モデル
さらに、以上の二つのモデルを組み合わせ、ドライバー、中立、両変異をもつ複合モデルを構築し、超並列シミュレーションによるパラメータ依存性解析から、バラメータによって、直線
的進化および分岐進化が起こる、二種のダーウィン進化が起こらない三つの場合および中立進化が起こる場合と起きない場合の組み合わせで大まかに６種類の進化ダイナミクスが観察された（図 1.21）。

28

図 1.21: 複合モデルを用いた解析の結果(e) 断続的変化を再現するモデル
以上により直線的進化、分岐進化および中立進化（の組み合わせ）が再現されたが、最後に、
断続的変化を再現する下記機構を上記のモデルに組み入れた。


増殖のために必要なリソースの限界を環境収容力として設定する。すなわち増殖確率を
g ・fN から g ・fN ・（1 – p/K）に変更する（g: ベースの増殖確率

f: ドライバー遺伝

子の強さ N: ドライバー遺伝子の数: p:集団サイズ K: 環境収容力）。


最大数（3 個）のドライバー遺伝子を獲得した細胞が分裂する際に小さい一定確率 e で環境
収容力を無限にする爆発的ドライバー変異(explosive driver mutation)を導入する。

前年度までの大腸がんにおける多領域シーケンス解析の結果より、早期大腸がんでは分岐進化、進行大腸がんでは中立進化が起きていることが示されているが、以上のモデルを用いて、この大腸がん発生過程における進化ダイナミクスのシフトを再現することを試みた（図 1.22）。このシミュレーション例の前半（水色矢印で示す）ではドライバー変異（変異プロファイル左の青色バンドで示す）がサブクローナルに存在する。次第に環境収容力の効果により増殖が遅くなるが、爆発的ドライバー変異を獲得した細胞（変異プロファイル上の緑色バンドで示す）が出
現し再増殖する。後半（赤色矢印で示す）ではドライバー変異はクローナルにのみ存在し中立
進化で腫瘍内腫瘍内不均一性は形成される。すなわち分岐進化から中立進化のシフトが再現さ
れている。

29

図 1.22: 断続的変化を再現する機構を組み入れたモデルによるシミュレーション(f)まとめ
以上シンプルなエージェントベースモデルを用いたシミュレーションにより、４種ががんの主な進化ダイナミクスを再現することに成功した。またこのモデルを用いたシミュレーションから大腸がん発生過程における分岐進化から中立進化のシフトが断続的変化によって起きていることが示唆された。以上の結果より、本研究により構築されたシミュレーションモデルはがんの起源と集団内多様性を理解するために有用なツールであると考えられる。

(3) 疾患に関連する胚細胞バリアントの効率的な解析システムの開発
平成 29 年度は、これまで検討を行ってきた低アレル頻度変異情報検出手法を発展させ、高感度胚細胞バリアント検出へ向けた、新たなベイズモデル統合法に基づく方法の開発を進めた。
サンプル細胞集団由来の DNA 上に存在する変異の検出は、まずそれらの DNA を短く断片化し、各断片の配列をシークエンサーで読み取り、長さ 100 文字程度の塩基配列情報を含む配列データ（リード）を大量に得た後に、それらのリード群をゲノム標準配列上にアライメントし、各リードがゲ
ノム上のどの位置に対応する DNA 断片から生成されたかを決定することから始まる。多数の細胞由
来の DNA 断片を読み取ることにより、ゲノム上の各位置は複数のリードにより観測される。通常、
胚細胞バリアントにおけるヘテロ接合性変異（ヘテロ変異）の場合、その変異箇所を被覆するリードのうち、変異を含むリードの割合の期待値は 50%となる。しかしながら、発生初期胚細胞での変異に由来するモザイク変異などの場合、変異を含むリードの割合が低くなり、検出が難しくなる。
これはがん組織由来のサンプル中にがん細胞含有率が低い場合や、がん細胞集団中のマイナークローンの変異検出における場合と共通の問題である。
ここで低頻度変異検出の問題を困難にしている大きな要因は、DNA に含まれる変異頻度が塩基の観測エラーの頻度と近づくにつれ、真の変異箇所と観測エラーによる見かけ上の変異箇所を区別することが難しくなることである。
これまで我々は、上記のような低頻度変異を精度よく検出するために観測したリード中に含まれ30

る、変異候補箇所の塩基頻度情報以外の付加的な有用情報を、ベイズモデルの中に取り込むことで、低頻度変異の検出精度を上げる手法の開発を行ってきた。Usuyama et al.,2014 においては、変異候
補箇所の近傍に存在する SNP で規定されるハプロタイプの情報を活用することで真の変異とエラーを区別するために、リードの生成過程を表現するベイズモデルを構成し変異検出を行っている（図1.23 (a)）。また Moriyama et al., 2017 では、ペアエンドリードで同一 DNA セグメントの読み取り
領域にオーバーラップが存在し、その領域内で変異を検出する場合、DNA セグメントを独立に二回観測することができるという事実に着目し、リードの生成過程をベイズモデル化し、それらの観測情報からエラーの場合と真の変異を区別することを目指す手法の提案を行っている（図 1.23 (b)）。
図 1.24 は、実際に観測された、(a)変異候補箇所近傍の SNP、および、(b)重複領域における変異候
補の例である。

図 1.23: リード中に含まれる変異検出のための付加的有用情報源。
（a）変異候補箇所近傍の SNP に規定されるハプロタイプ。（b）ペエンドリードの重複領域における変異シグナル。（c）変異支持リードのストランドバイアス。

31

図 1.24: 実際に観測された例。(a)変異候補近傍のヘテロ SNP。(b)ペアエンドリード重複領域における変異候補。
上述のリード配列に含まれる有用な情報は、変異リードの向きの偏り（図 1.23 (c)）など、他にも考えることが出来る。そこで、それらの複数の情報源を同時にモデルに取り込み組み合わせることができれば単独モデルに比べた性能向上が期待される。しかしながら、これまで変異検出の観点で、複数のベイズモデルを一つのベイズモデルにまとめるモデル化はなされてきていなかった。複数の
ベイズモデルをまとめる代表的な方法としては、ベイズモデル平均化 (Hoeting et al.1999)があげられる。この方法では各モデル内以外にも、複数のモデル間に新たに事前確率を設ける必要があり通常この事前確率の設定は、事前確率自体の解釈が難しく先験的に設定することが容易ではない上、変異検出においては性能に大きく影響を与えてしまうなど、今回の目的においては適用が困難であった。

32そのため、ベイズモデル平均化に基づかないモデル統合のための新規方法論の開発を行った。本手法では利用できる情報源の組み合わせごとにリードデータ集合を分割し(data partitioning)、それぞれの分割データ集合に対して生成モデルを作成する（図 1.25）。各生成モデルの作成には、既存の手法を参照することも出来る。ある変異候補が、変異またはエラーであるかの判断（変異コール）には、二種類の統合モデル（変異モデルおよびエラーモデル）を考え、両統合モデルの周辺尤度の比（Bayes Factor）値を用いる。各周辺尤度の計算には、リードごとに各データ分割空間への所属確率を設定する必要があるが、その先験的な設定には困難を伴う。本手法では変異モデルとエラーモデル間で、これらの所属確率に違いがないという緩やかな仮定を置くことで、当該所属確率を、
変異コールに当たっては陽に設定することを回避でき複数モデルの統合を容易にしている。同手法
は、モザイク変異検出およびがん細胞中の体細胞変異検出、それぞれの問題に対して適用可能である。実際に上記の手法に基づき、実際に低頻度変異の検出において鍵となる情報源を複数（ハプロ
タイプ情報、リード重複、ストランドバイアス）考慮に入れることのできる生成モデルを構築し統合を行った（図 1.26）。結果に関しては、シミュレーションデータと実データを用いた評価を行い論文投稿中である。

図 1.25: データ分割に基づく生成モデル。k : データ分割空間のインデックス。xn : n 番目のペアエ
ンドリード。tn : n 番目のペアエンドリードのデータ分割空間所属指示変数。

33

図 1.26: (a)開発手法に基づき作成したハプロタイプ情報、リード重複、ストランドバイアスの情報を統合した生成モデル。(b)近傍ヘテロ SNP 情報のみが使える場合のリードのパターン。
この疾患に関連する胚細胞バリアントの効率的な解析システムの開発は、ウイルス感染によるがんを始めとする疾患の理解には極めて重要な位置をしめるもので、今後、画期的な成果が得られる
ことが強く期待できる。

34

2.

データ同化生体シミュレーションによる個別化医療支援（サブ課題 B）

本節の構成は以下の通りである。
2.1

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

2.2

循環器系個別化医療支援に向けたデータ同化シミュレータの開発における(5)と(6)2.3

個別生体機能評価のためのデータ同化シミュレータの開発

では、本成果報告書４－１．実施計画の２．
「(1)個別化医療支援のための生体シミュレータの開発および評価」についての成果を述べる。
また2.3

個別生体機能評価のためのデータ同化シミュレータの開発における(1)、(2)、(3)、(4)では、同じく実施計画の２．「(2) データ同化手法の開発と個別化医療支援への応用」についての成果を述べる。
2.1 全脳循環シミュレータの開発(1) 脳大域循環モデルの開発と評価（大阪大学

伊井仁志、和田成生）A) 半仮想血管網における脳血流解析
脳では全身で消費される酸素の約 20%が血液循環により絶えず運ばれ消費されており、この達成には正常範囲内の動脈圧においてその大きさに依らず常に同程度の脳血流量が保たれる自動調節機構が大きく関わっている。本機構において、血管径の能動的制御に起因する血流抵抗の変化が関連するとされているが、複合的なネットワーク構造を有する脳血管網において、脳全域で恒常的に安定な循環場を達成する仕組みは未だ明らかではない。現状、実験観察的に大域・局所脳循環を同時に評価する技術は確立されておらず、また、複雑に連成する生体反応下において、個別に関連要素を評価するのは容易ではない。この問題に取り組むため、前年度に開発したヒト脳血管の医用画像と数理モデルを組み合わせた脳血管網モデルを発展し、半仮想的に構築した脳血管形態において、ボクセル型数値解法に基づく大規模計算アプローチにより脳血流場を数値的な再現を試みたので報告する。
脳血管網モデルの概要は以下の通りである。


ヒト脳の 4D-CTA 画像に対し、画像処理により約 1 mm 径までの血管を 3 次元構築する。
この際、血管の繋がりが途切れている箇所を参考に動脈系と静脈系を分離する。



脳実質に相当すると考えられる脳表面形状を三角形ポリゴンの集合として構築する。



Gödde と Kurz (2001, Development Dynamics, 220, 387)の血管成長数理モデルを参考に、
医用画像より構築された血管末端座標を始点とし、脳表面ポリゴン上において血管セグメ
ントを数理的に成長させる。この際、表面ポリゴンを二層に分離することで、画像構築された血管系および動静脈系が重ならないようにする。



構築された動静脈末端を各層に接合することで閉塞された血管網を構築する。



血管網はセグメントで構成されており、これを設定した半径を持つ円筒とすることで、流体解析に必要となる流体体積率を 3 次元ボクセル格子上で定義する。

本手法により構築された脳表層血管までの脳血管網モデルを図 2.1.1 に示す。赤色で示しているのが動脈、青色が静脈であり、それぞれ医用画像および数理モデルにより構築された血管系を35

示している。医用画像から抽出された動静脈から脳表面に沿う動静脈が構築されている様子が確認できる。

図 2.1.1: 脳血管網モデル
血流解析には前年度より開発を続けているボクセル型流体解析シミュレータ(Voxel based
Brain Blood simulation -Flow: VBB-Flow)を用いた。ここで、計算資源削減および計算効率化のため、解析領域のみをカバーする MPI-box を生成した。図 2.1.2 に主要な脳動脈のみを対象とした MPI-box の様子を示す。これにより、系全体を立方体で覆うことに比べ、約 95%の MPI
ノード数を削減することに成功した。なお一つの MPI-box は 32×32×32 のボクセル格子で構成されている。また、ボトルネックとなるポアソン方程式の求解に幾何的マルチグリッド法の導入することで反復回数の削減を行っており、問題規模に依存するが、実行時間として 10 倍ほどの性能向上を達成している。

図 2.1.2: 脳主要動脈を解析対象とした MPI-box
3 次元流体シミュレータにより、構築された脳血管網における血流場を数値的に再現した。内頚動脈および椎骨動脈のそれぞれ入口に一様流を課し、内頚静脈出口で圧力を固定した。脳大域循環シミュレーション結果に関して、せん断応力分布を図 2.1.3 に示す。ここでは、左側の脳循環の様子を表示している。全脳循環レベルにおいて動脈入口から静脈出口まで流れ場が形成され、壁面せん断応力の値が妥当な範囲で求められた。本解析では、仮想的に脳動脈瘤を設けており、その付近を拡大した結果を同じく示している。大域的な脳循環場において局所的な 3次元血流動態を求めることが出来ている。

36

図 2.1.3: 脳大域循環シミュレーション結果（脳循環全体の壁面せん断応力と脳動脈瘤内の断面速度分布および壁面せん断応力）
図 2.1.4 に脳血管全体での壁面せん断応力の度数分布を示す。分岐血管においては、壁面せん断応力が等しくなるとされる Murray の法則が一般的に成立するとされているが、結果からは
かなりのばらつきを持つことが分かった。また、一部の領域において、壁面せん断応力が非常に大きな値をとる結果となった。これは今後の検討課題であり、血流量の等分配性と併せて検討を行い、生理学的に妥当になるための条件を検討することで脳循環メカニズムの解明につなげていく必要がある。

図 2.1.4: 脳循環における壁面せん断応力の度数分布「京」のシステムにおいて、およそ 16000 ノード(MPI)×8 コア(OpenMP)の並列解析まで達成している。任意形状に対するボクセル型計算手法の特性上、MPI プロセス毎に流体領域となる
ボクセル数が異なり計算コスト負荷にアンバランスが生じるため、並列性能の低下を招く。次年度では、プリ処理およびアルゴリズムを見直すことで、これらの改善を図る予定である。
B) 血管形態モデルの改良
全脳循環を数値的により正しく再現するためには、脳血管網の形態的特徴である、(i)脳表層に沿った経路、(ii)Willis 動脈輪や脳表層での血管吻合など様々なスケールでの複合的な分岐構造、(iii)脳表層での領域ユニット内の動脈・静脈の接合、を考慮した脳血管網モデルが必要となる。

37

これまでに用いてきた脳血管網モデルでは、数理モデルにおける血管径が全て等しくなるため、生理学的に妥当な脳循環の再現には課題があった。そこで、全脳循環場の数値再現に向け、医用画像と血管網数理アルゴリズムを組み合わせ、上記で挙げた形態的特徴(i)-(iii)を満足する新
たな脳血管モデルの構築を目指しアルゴリズムの開発を行ったので報告する。
ヒト脳の医用画像より取得した数 mm 径までの動静脈の末端座標および脳表面形状(Reddan,
2015, https://www.thingiverse.com/thing:801633)を利用し、血管網数理アルゴリズムにより脳表層までの動静脈を構築する。脳表面形状は三角形要素の集合で与えられており、任意に設定した個数分の三角形を無作為に選びだし、これを起点にした領域拡張より脳血流領域ユニットを決定する。この際、各領域ユニットの面積および周長の重み付き平均からなる評価関数を擬
似焼きなまし法により最小化することで、等面積かつ等方的な領域ユニットを構築する。構築された領域ユニットを図 2.1.5 に示す。

図 2.1.5: 脳表面の領域ユニット
これを利用し、Karch ら(1999, Comput Biol Med, 29, 19)の血管生成アルゴリズムを基に、画像処理により得られた数 mm 径血管の末端座標を始点として、脳表層までの血管網を数理的に構築する。Karch の手法では、無作為に加えた端点と既存管路からなる新たな分岐管路が構築されるが、この際、いくつかの仮定のもと経路に依存する管径を評価し、管径路全体の体積が最小となる最適化問題を解くことで分岐管路が構築される。本研究では、最適化問題の簡略化のため、加える端点と最短距離にある管路と分岐管路を形成し、その形成パターンは 4 種類に限定している（図 2.1.6）。加える端点は、無作為に選びだした脳領域ユニット内の脳表面三角形の重心点とする。また、分岐管路の形成には、管路長の範囲などいくつかの制約を設けており、これが満たされない場合は分岐管路の生成を棄却する。分岐経路を作成する際、3 次元空間において脳表面形状の符号付き距離関数を離散的に定義し利用することで、脳表面から一定の範囲にある経路を構築し、さらに脳表面を貫くような分岐経路が提示される場合は線分を細分割することで脳表面に沿うように位置座標を修正し曲線管路を作成する。
静脈系の構築：静脈端点（三角形）は定めた領域ユニットのいずれかに属し、端点間で領域ユニットの重複がないように制約を設ける。全ての静脈端点に経路が作成されるよう、静脈端点数≤領域ユニット数と設定する。
動脈系の構築：脳皮質に穿通する血管網構造として、領域ユニット内では一つの静脈に複数の38

動脈が存在することが知られている。この特徴を考慮し、静脈端点が存在する領域ユニットのみから三角形を無作為に選び動脈端点とし経路を構築する。この際、一つの領域ユニットには設定した個数まで動脈端点を許す。

図 2.1.6: 端点追加時の血管分岐パターン
図 2.1.7 に静脈と動脈それぞれに対し医用画像による血管形状モデルと本アルゴリズムを用いた数理モデルを合わせた結果を示す。全ての血管を可視化すると煩雑になるため、数理モデ
ルは直径 0.4 mm 以上の血管を可視化している。符号付き距離関数の勾配を求め脳形状の内側にある点を脳表面に修正することで、脳表層を覆うように血管網が構築され、これまでに用いた脳血管網モデルと異なり立体的に経路が構築されていることが分かる。また、理想化された流れの支配方程式を適用しているため、経路に応じて血管径が小さくなる様子を再現できている。

図 2.1.7: 新規アルゴリズムにより構築された脳血管網モデル
本アルゴリズムでは、定めた制約により静脈の領域に動脈が一つも合流しない状況が生まれる。その程度を確認するために、静脈末端が存在する各領域ユニット内で動脈/静脈の比に関する度数分布を示す（図 2.1.8）。結果より、動脈が合流しない領域が 1%未満で存在することが分かった。このような領域では、領域内での動静脈の接続ができず全脳循環の解析をすることができないため、血管系の修正が必要となる。また、現在は血管同士で重なりが起きるため、これらの点を解消し、今後は領域ごとの末端からさらに領域内で接合するアルゴリズムを作るこ39

とによって、完全閉塞系の全脳循環モデルの構築を達成する。

図 2.1.8: 静脈末端が存在する各領域ユニット内で動脈/静脈の比に関する度数分布(2) 脳微小循環モデルの研究開発(大阪大学

伊井仁志、武石直樹、和田成生）A) 単純化された毛細血管網における血球分配の検討
全脳循環の機能において、微小循環の役割が注目されている(Blinder et al., 2013, Nat Neurosci,
16, 889)。血液中の固体成分である赤血球は酸素運搬を担うため、脳皮質内の細胞代謝活動・神経系活動との関連より、その分配挙動の理解は重要である。分岐血管における血球分配に関する物理的解釈はこれまでになされているが(Pries et al., 1989, Microvasc Res, 38, 81)、微小循
環スケールでの複雑な毛細血管網における物理的解釈はいまだに統一的な見解がなく、酸素輸送との関連も未解明な点が多い。その理由として、現象における受動的・能動的な寄与が整理されていないことが挙げられる。これに対し、毛細血管網の簡易形態モデルを用いた 3 次元数値シミュレーションを行い、血球分配における力学的な寄与を明らかにし酸素輸送との関連を検討したので報告する。
当研究グループがこれまでに開発してきた直交格子を用いた流体・弾性カプセル力学連成解
析手法(Ii et al., 2012, Commun Comput Phys, 12, 544)を拡張し、血漿と赤血球の相互作用を取り扱う。赤血球膜の力学特性は超弾性膜モデル(Skalak et al., 1973, Biophys J, 13, 245)で表現し、血漿および赤血球内部はニュートン流体に従うとする。この際、赤血球内部には高粘度の液相を仮定する。毛細血管網の形態構造を簡易的に模擬した流路モデルを作成し、入口に速度、出口に圧力を与える（図 2.1.9）。ここで、入口すぐの領域において、赤血球を一定時間毎に与えることで、赤血球が流路全域に行き渡った際の平均ヘマトクリット(HT,whole)を制御する。
流路網の中央部に流路 St、Sm、Sb を定義し、各流路における排出ヘマトクリット(HD=QRBC/Q)を評価する。ここで、Q および QRBC はそれぞれの流路における全流量と血球流量である。本解析では、各流路への血球分配に着目する。上下方向の対称性を考えると、流路 St と Sb への血球分配は時間平均で等価となり、流路 Sm との関係は解析条件に依ると考えられる。そこで血球分配指標(ΔHD*)を定義し評価する。

40

図 2.1.9: 毛細血管網の簡易流路モデル
図 2.1.10 に、異なる HT,whole において、十分に赤血球が流路に行き渡った際の血球分配の様子
およびΔHD*の結果を示す。定性的にはどの条件においても流路全域に血球が分配している様子が分かるが、HT,whole を大きくするにつれΔHD*の中央値およびばらつきが小さくなっていることが確認できる。図 2.1.11 に流路 St、Sm、Sb における流量と圧力の比である血管抵抗(R=dP/Q)と局所的な排出ヘマトクリットの関係を示す。排出ヘマトクリット HD の増加に従い、血管抵抗がばらつきを持ちながら上昇している様子が確認できる。ここで、直円管において十分発達した流れ場に対し実験・解析を行った Pries ら(1989, Microvasc Res, 38, 81)の経験則から外れてきており、ネットワーク構造においては、単純に経験則を適用することが出来ない可能性が示された。結果より、本解析が対象としているストークス流れでは、分岐管における流量分配は下流側の流動抵抗によりほぼ決定され、これを担うのは下流側流路での局所的な赤血球数であると解釈できる。流路内の赤血球が一時的に増大することで流動抵抗が上昇し他の分岐管に血球が分配されやすくなるが、血球が分配されると同様の現象が起こるため、一定値を基準に時間変動をともなう結果となる。よって、流路全体の赤血球数が多くなる(Ht が高くなる)ほど、血球による流動抵抗変化がより働き、流路毎の血球分配量に差がなったと考えられる。これらの結果は、毛細血管毎の血球量を一定に保つような調節機構が受動的に働くことを示しており、組織への恒常的な酸素供給と密接に関連することが示唆された。

図 2.1.10: 異なる全体ヘマトクリット(HT,whole)と血球分配(ΔHD*)の関係

41

図 2.1.11: 排出ヘマトクリット(HD)と流動抵抗(R)の関係
開発中のボクセル型流体・固体解析プログラム(VBB-Flow)を拡張することで派生アプリ
VBB-Capsule-Eulerian を開発した。領域分割による MPI 並列において、各領域に含まれる計
測点を分割コアに割り当て通信を行った。本アルゴリズムは完全オイラー型定式化になっているため、直交格子の並列化の利点を最大限活用することができる。これまでに、「京」のシステムにおいて、200 ノード(MPI)×8 コア(OpenMP)のハイブリッド並列を達成しており、次年度
以降、前年度までに開発した Immersed

Boundary 法を用いた手法と並行して開発を進め演算効率の比較などを行っていく。それと並行して解析系の大規模化を行い、実際の脳微小循環形
態を用いた解析に適用していく。
B) 血球流動と酸素輸送の連成解析手法の開発
生体の恒常性を担う血液中の酸素輸送は、血球流動に関する力学と血球内部タンパク(ヘモグロビン)の(脱)酸素化に関する生理学が連成したマルチフィジックス問題を構成する。これまで、血球流動の詳細な力学解析は行われてきたが、血球流動と組織への酸素輸送を連成した解析は未だ統一的なアプローチがなく、微小血管網における血球流動場と組織内の酸素消費場との物
理的関連および生理学的意義は未だ解明されていない。そこで、血球流動による酸素輸送の詳細を明らかにすることを目的に、赤血球・血漿・組織の各領域における酸素輸送を直交格子上で統一的に取り扱う混合化方程式モデルを提案し、ヘモグロビンの生理学反応を考慮した新たな酸素輸送の数理モデルを構築したので報告する。
赤血球内部・血漿・組織の各領域における酸素濃度に関する反応拡散方程式(Lücker et al.,2015, Am J Phys - Heart Circ Phys, 308, 206)に対して、各領域界面における酸素分圧と流束の連続条件を適用し、相関数を用いて混合化を行うことで、ヘモグロビンの酸素飽和度・脱酸素化率および組織領域での酸素消費率を考慮した混合酸素濃度に関する反応拡散方程式が導出する。ここで、溶解する酸素濃度と分圧の関係には Henry 則を適用している。
提案手法の妥当性検証を行った。球対称系の計算領域を考え、血球・血漿およびそれらを覆う組織層を考える。領域端にはノイマン境界条件を課し、初期条件として、血球内部のみに酸素濃度(C=4 µlO2/cm3)を与える。溶解度係数α、拡散係数 D は先行研究(Lücker et al., 2015, Am
J Phys – Heart Circ Phys, 308, 206)に基づいて決定した。図 2.1.12 に解析系、異なる時間(t=1042

ms、100 ms)での動径方向の酸素濃度分布を示す。ここで、空間解像度Δx を変化させた際の結果と球対称を仮定して求めた参考解と比較している。空間解像度を上げるにつれ結果が参考解に収束していく様子が分かる。また、Δx~0.2 程度あれば十分な精度で求解可能であることが分かった。

図 2.1.12: 酸素輸送解析手法の妥当性検証
次に、前年度に計算された直円管内を流動する赤血球形状を用いて酸素輸送に関する予備的な計算を行った。なお、本解析では力学系とのカップリングは行っていないため、血管外側の組織領域に一様速度を課すことで、血球群が変形せず移流していく様子を再現している。ある
時刻における酸素濃度を図 2.1.13 に示す。組織領域での酸素濃度は、0.1 s 以内にヘモグロビン
から放出された酸素によって急上昇し、1.0 s 後には酸素は組織内で消費され、再び初期濃度に近い値まで戻った。一方、赤血球領域の酸素濃度は初期状態から約 10%にまで減少した。これ
らより、無酸素状態での酸素供給の過渡的状況の検討を行うことができた。今後は、これまでに開発している力学連成手法とカップリングすることで、血球分配が酸素輸送にどのように関与するか検討していく。

図 2.1.13: 直円管内の酸素濃度分布の予備計算結果
2.2 循環器系個別化医療支援に向けたデータ同化シミュレータの開発(1) 脳循環系における血流同化手法の開発（大阪大学

伊井仁志、和田成生）A) フィードバック法による血流データ同化手法の開発と評価
位相コントラスト MRI(Phase-contrast magnetic resonance imaging: PC-MRI)は時空間に分布する血流速度ベクトルを非侵襲的に取得できる計測方法である。近年、データ同化を用いることで計測速度を数値流体力学(CFD)解析に反映させる試みがなされている。しかし、時系列データに対する非定常解析は、定式化に由来する計算負荷の高さから未だに発展段階である。こ43

れに関し、昨年度までに開発を行った、フィードバック制御における入力による影響を圧力境
界条件 に組み込ん だ物理的整 合性を満足 するアプロ ーチ (physically consistent feedback
control based data assimilation: PFC-DA)を適用し、時系列 PC-MRI 計測データを用いた解析を行いアプローチの実現可能性を示したので報告する。
流体の支配方程式には非定常非圧縮性 Navier-Stokes 方程式を採用する。血流は出入口をともなう閉空間の流れであるため、入口・出口に圧力境界条件を与え、この空間分布を推定する。
物理的整合性を満足するフィードバック制御に基づくデータ同化手法(PFC-DA)では、計測とCFD の速度差に由来するフィードバック力をスカラーポテンシャルにより近似的に表現し圧力境界値と関連付ける。ここで、計測速度と CFD 速度との偏差をフィードバック制御の一つである PD 制御の概念に基づき外力項を評価する。これをソース項とするスカラーポテンシャルのポアソン方程式を解くことで、入口・出口の圧力値として採用する。詳細は前年度の研究結果を学術論文(Ii et al., 2017, 34, e2910)として報告している。これに加え、本アプローチでは CFD
解析で用いる時間刻みで計測データが定義される必要があるため、時系列データを参照解として利用するために、血流の周期性を利用しフーリエ級数展開により全計測点における速度の時間変化を予め定めておく。
内頸動脈を対象とし、CT 装置により取得し DSA 画像から画像処理ソフト AMIRA5.5.0(FEI
VSG, ZIB,Germany)を用いて血管形状を抽出し体積率関数(VOF 関数)を構築した。MRI 装置
により取得された PC-MRI 画像から計測速度ベクトルを構築した。速度エンコーディングは 1
m/s と設定した。また、CT 計測と MRI 計測の手技の違いにより位置ずれが生じるため、剛体の並進・回転にともなう 6 自由度成分を離散的に調査し、形状を評価するパラメータを最小にする値を採用することで位置ずれを修正した。 DSA 画像の解像度は 0.1 mm、PC-MRI 画像の解像度はおよそ 0.5 mm であった。PC-MRI 計測では一心拍(周期 T=1.47 s)あたり 15 フレー
ムの速度情報が得られるため、7 次の展開モードでフーリエ級数展開により連続時系列データを作成した。なお医用計測は大阪大学が定める所定の倫理手続きを経て大阪大学医学部により行われたのち、匿名化を行った画像が提供された。CFD 解析のボクセル格子幅は 0.2 mm に設定
し、最大流速において CFL 条件がおよそ 0.1 となるように時間刻みを設定した。
図 2.2.1 に、収縮期ピーク(Peak Systole)および拡張期末期(End Diastole)における計測点での PC-MRI 速度、推定速度およびそれらの偏差を示す。結果より、主血管において比較的大きな速度差があるものの、推定速度場は大域的な計測速度場を再現した。図 2.2.2 に推定速度および計測速度の空間平均値の時間変化を示す。速度の空間平均量に関する推定結果は計測速度の時間変化を良く再現した。さらに、図 2.2.3 に収縮期における全 CFD 格子点での推定速度ベク
トル、これを用いて評価された壁面せん断応力を示す。結果より、少ない計測データを用いた場合でも、CFD で採用した高解像度場において非物理的な振動なく速度場が推定されている様子が確認でき、それらにより評価される壁面せん断応力も生理学的に妥当な範囲で得られることを示した。計測データとの速度差の平均値は一周期を通じ 30%ほどであった。直円管を用いた妥当性検証の結果より、誤差を含む計測速度との差 30%程度でも実際の解との差は 10%以下である状況が存在することを確認している。よって、本解析は妥当な結果である可能性が高い
が、計測との速度差は時相および速度絶対値により異なるため、様々な計測ノイズを模擬する

44

ことで妥当性の検証を引き続き行い推定手法の精度検証を引き続き行う必要がある。

図 2.2.1: 収縮期（上図）および拡張期（下図）における計測点での PC-MRI 速度(Simulation)、
推定速度(PC-MRI)およびそれらの偏差(Difference)。

図 2.2.2: 計測速度(PC-MRI)および推定速度(Assimilated)に関する空間平均値の時間推移

図 2.2.3: 収縮期における全 CFD 格子点での推定速度ベクトル（左図）、これを用いて評価された壁面せん断応力（右図）。

45

本手法により、既存研究で指摘されているような、血流の時間変化に起因するような血行力
学因子を評価することが可能となった。図 2.2.4 に脳動脈の時間平均された壁面せん断応力および Oscillatory Shear Index(OSI)を示す。これにより、複数の個別血流に対する解析より脳動脈瘤のリスク評価と血流動態の関連を検討することが可能となったため、次年度以降取り掛かっていく。

図 2.2.4: 脳動脈瘤の時間平均された壁面せん断応力(Time-average WSS)および Oscillatory
Shear Index(OSI)
B) 変分法に基づく血流データ同化手法に関する研究開発
クモ膜下出血の主要因である脳動脈瘤の破裂およびその発生・成長には、壁面せん断応力(WSS)など血行力学因子の関連が指摘されている。非侵襲的な血流計測法の一つである PC-MRI により得られた計測速度を利用したデータ同化血流解析は有用な手段であるが、いくつか課題もある。PC-MRI 撮影時には、最大血流速をカバーできるよう Velocity encoding(VENC)が設定されるが、この VENC 値と計測ノイズは密接に関連するため、血流速度が遅い動脈瘤領域に撮影範囲を限定し VENC 値を小さく設定することで計測ノイズを下げることができる。データ同化
においては、データノイズが推定精度に影響を与えることから、脳動脈瘤領域のみを解析対象とすることで、瘤内の血流場を高精度に推定することが可能になると考えられる。また瘤内のみを解析対象とすることで、計算コストの削減が可能となり、より多くの個別データを用いた解析に繋げることができる。この解析においては、瘤内への流入・流出部の境界条件を推定する必要が出てくるが、これまでに我々が提案しているフィードバック手法および圧力境界条件のみの推定では、空間連続的に変化している流れ場への境界条件の自由度が低いため推定精度が低下するのは自明である。これを解決するため、変分法によるデータ同化アプローチの新たな定式化を提案したその妥当性を検証することで、脳動脈領域のみを対象としたデータ同化血流解析の実現可能性を示したので報告する。
流体領域において、連続の式と定常非圧縮性 Navier-Stokes 方程式により流れ場を記述する。
境界条件として、出入口境界の流速、壁面で no-slip 条件を与える。モデル速度および計測速度の偏差に関する空間的な L2 ノルムからなる評価関数を定義する。流体の支配方程式を拘束条件とし、評価関数が最小となるように出入口境界の流速ベクトルを制御変数とし推定する。
Lagrange の未定乗数法により汎関数を導入し、その第一変分より停留条件をまとめることで、元の流体方程式に付随する強形式で記述された随伴方程式および制御変数に関する関連式が導出される。制御変数に関して最急降下法などの反復手法を適用することにより、流体方程式および随伴方程式を独立して求めることが可能となるため、通常の変分法で現れる大規模な連立46

方程式の求解を回避することができる。また、流体解析には当研究グループが開発している
VBB-Flow を用いる。
手法の妥当性検証として、脳動脈瘤を含む実形状血管モデルに対し、母血管に流入口でRe~300 となるよう一様流速を与えて順解析を行い、定義した瘤領域において順解析結果を間引いて抽出し瘤内の参照速度を取得し、瘤形状においてデータ同化解析を行った(図 2.2.5)。つまり、速度を推定する出入口境界は瘤形状の底面となる。参照解の解像度は PC-MRI を想定して0.5 mm とし、正規分布に従う速度ノイズを、ノイズレベル(SNR)を変化させ参照解に加えた。
流体計算の解像度は 0.1 mm とした。図 2.2.6 では、瘤断面における順解析およびデータ同化解析による速度絶対値と出入口境界(底面)から 0.2mm 上方で評価した壁面せん断応力(WSS)の箱
ひげ図を示している。瘤内速度は空間分布も含め定性的に良く一致し、また参照速度との正規化誤差ノルムはノイズを含まない疑似計測データを用いた場合、約 6%となり定量的な一致も確認できた。さらに、速度場より評価される WSS の分布もノイズレベルを高くした場合においても良い一致を示した。

図 2.2.5: 母血管を含む順解析系(a)と抽出した動脈のデータ同化解析系(b)および疑似計測速度(c)図 2.2.6: 各ノイズレベルでの疑似計測データを用いたデータ同化結果
次に、実際の PC-MRI の計測データを用いた解析を行った。速度ベクトルや形状との位置合
わせはフィードバック手法による解析と同じである。図 2.2.7 にデータ同化解、PC-MRI および

47

その誤差ベクトルを示す。両者は定性的に良い一致を示した。ここで、誤差ベクトルの正規化された L2 ノルムは 34%程度であり、解析対象は多少異なるもののフィードバック手法による結果と同程度であった。提案した変分法によるアプローチは、速度ベクトルの境界条件を推定しているため、スカラー値である圧力を推定するより自由度が高く、より高精度に速度場の推定が行えると考えられる。それにもかかわらず同程度の推定精度となったことから、計測データの誤差およびノイズそのものに何らかの問題点があることが示唆された。これは今後の検討課
題であるが、変分法に基づく新たな血流データ同化手法の提案により、瘤内部のみを対象とした解析の実現可能性について示すことができた。

図 2.2.7: PC-MRI データを用いたデータ同化結果(2) 血管造影シミュレータの開発（大阪大学

伊井仁志、和田成生）

数値シミュレーション結果に対し擬似的に計測データを再現することで、シミュレータの評価、シミュレーション結果がどのように計測されるかの評価、また計測データを用いたデータ同化シミュ
レーションが可能となる。今年度は、血管造影計測に着目し、血液と造影剤の混相流モデルを用いた造影剤注入の数値シミュレーション手法を開発した。直円管モデルを用いて造影剤動態の基礎的検討を行い、種々の解析条件における数値解析結果と実測データとの関連を循環機能評価の観点から検討したので、これらに関して報告する。
血管造影は、脳血管障害の診断における狭窄部の特定や動脈瘤の評価や冠動脈疾患の診断における狭窄部の特定などに用いられており、循環機能の評価手法とし広く採用されている。血管造影では、注入された造影剤が血液と混合して血管内を流れ X 線により撮影されるが、血管内の造影剤の分布に偏りが生じた場合や造影剤が流れない領域が存在する場合、血管形状を正確に捉えることが困難となり診断精度が低下する恐れがある。そのため、注入した造影剤を素早く血液と混合させ、血管内に均一に分散させることが望まれる。管内の造影剤の注入現象は管内噴流であり、注入された造影剤が周囲の血液を引きずりながら流れていく。単相流れにおいて管内噴流は古くから研究されているが、血液と造影剤の二相成分がある中で流体力学的な無次元パラメータが混合挙動に与える影響については十分に検討されていない。本研究では、血液と造影剤の混相流モデルを用いた造影剤注入の数値シミュレーションを行い、直円管モデルを用いて様々な条件と混合挙動の関連を明らかにする。
血液および造影剤の混相流れとして混相流モデルを用いる。このとき、等エントロピー、低 Mach数を仮定した非定常微圧縮性流れを考える。支配方程式は、Navier-Stokes 方程式と微圧縮性項を加えた連続の式である。離散化には Fractional-Step 法を拡張した手法(Ii, Wada, 2017, J Biomech, 50,

48234)を用いる。血液と造影剤で体積平均した流体の密度および粘度を用いることで、混相流れの定式化を行う。本解析で用いるパラメータについて述べる。血液は密度 1.05 kg/m3、粘性係数 3.5×10-3
Pa･s とする。また造影剤は、ヨード性造影剤と生理食塩水の混合物を使用するとし、
密度 1.2 kg/m3、
粘性係数 2.2×10-3 Pa･s とする。また、本解析における混合速度は O(100) m/s であり、音速 a は液体中では O(103) m/s であるため、Mach 数は十分に小さい。ここで、仮に音速を O(102) m/s としても Mach 数は十分に小さく、微圧縮性の性質は変わらないと考えられる。そこで本研究では、離散化により得られる連立方程式の求解における計算コスト削減のため a = 500 m/s とする。また、本解
析では重力の影響を考慮しない。時間刻み幅を 5.0×10-6 s とし 2.0 s 間の造影剤注入を再現する。
本解析に使用した形状モデルを図 2.2.8 に示す。血管を直円管に単純化し、入口にカテーテルを模した形状を作成する。直円管は内径 5 mm、長さ 112 mm とし、カテーテルの形状を外径 2.5 mm、内径 1.5 mm、長さ 2.0 mm とする。本形状モデルを含む領域を一辺が 0.1 mm の直交格子で分割する。
図 2.2.8 に示す左側断面において、カテーテル外部を血液の入口境界、カテーテル内部を造影剤の入口境界とし、それぞれの入口境界に異なる境界条件を与えることで、血流に対する造影剤の注入を表現する。また右側断面を出口境界とする。

図 2.2.8: 造影剤注入シミュレーションの計算系：直円管と設置されたカテーテル先端の様子
実際の血管造影では造影剤の輸送が 3 次元的に視覚化されるのではなく、X 線透視撮影を行うことで 2 次元の血管造影画像が取得される。そこで、造影剤注入シミュレーションで得られた結果から、図 2.2.9 に示すような X 線による x-y 平面への平行投影を考え、2 次元の血管造影画像を疑似的に取得する。

図 2.2.9: 血管造影における疑似投影データ作成の概念図
血流の Reynolds 数(Reb)および造影剤注入の Reynolds 数(Rec)を変化させ 20 通りの解析を行った。
ここで、血管内への造影剤の注入現象は管内噴流と見なせることから、噴流(造影剤)と周囲流れ(血49

流)の運動量に関する無次元パラメータである Craya-Curtet 数(Ct)を用いて、造影剤の分散状態と逆流の有無を整理した。結果を図 2.2.10 に示す。ここで、図中の印の違いは異なる分散状態を示しており、×は造影剤が分散せず噴流のまま流れ去った状態、△は血流部に逆流が生じず造影剤が管断面全体に分散した状態、○は血流部に逆流が生じ造影剤が全体に分散した状態である。図 2.2.11 に
主流方向の速度のカラーマップを示す。これより、逆流の有無が確認できる。Rec = 818 では、Reb(および Ct)の値に依らず造影剤が分散しなかった。しかし、流れの散逸により噴流のコア部は注入部からの距離に伴い管内に広がっていくことが期待される。ここで減速区間の中間点を造影剤が分散する位置と考え、造影剤の分散が確認された 15 通りの解析についてその座標点を調べたところ、Rec
が大きくなるほど概ね値が小さくなる傾向が見られ、造影剤が分散するまでの距離が短くなることが確認された。造影剤が分散した状態に着目すると、本解析においての Ct 値が 0.39 以下では全ての場合で逆流が生じており、0.6 以上では逆流が生じていないことがわかる。このことから、逆流の有
無は Ct の値によって分類できると考えられ、逆流が生じるための Ct の閾値は 0.4 ~ 0.6 程度であると示唆された。単相流において逆流が生じるための Ct の閾値は 0.7 ~ 0.98 と報告されており(Barchilon, Curtet, 1964, J Basic Eng, 86, 777; Revuelta et al., 2004, Phys Fluids, 16. 208)、本解析で得られた値とは異なっているため、二相流による影響が関連していると考えられる。

図 2.2.10: 造影剤分散に関して Reynolds 数(Rec)と Craya-Curtet 数(Ct)の関係

図 2.2.11: 異なる解析条件において、ある時間で平均をとった主流速度のカラーマップ

50

図 2.2.12: 造影剤注入の様子(左：ボリュームレンダリング、右：疑似投影データ)
主要脳動脈の一部を用いた血流解析において、血管造影シミュレータにより得られた結果を図2.2.12 に示す。内頚動脈から注入された造影剤が混合しながら脳血管に流入していっており、実際の血管造影画像に近い造影剤の注入の様子が再現されている。さらに、大域脳血管モデルを用いた血流解析において、標準モデルおよび狭窄モデルを用いた造影剤の投影図を図 2.2.13 に示す。血管の一部に狭窄を設けることにより変化した流動場が血管造影画像に反映されている様子を再現することができた。

図 2.2.13: 大域脳血管モデルにおける、標準モデルおよび狭窄モデルを用いた解析結果の疑似血管投
影画像
本研究により、流動場が造影剤流動および分散に与える影響を検討することができた。これらは、実際の医用計測において物理現象がどのように取得され評価されるかに密接に関連しており、数値シミュレーション結果の解釈および循環機能評価に向けたポスト処理としての研究開発への第一歩を示すことができた。
(3) 脳動脈瘤コイル塞栓治療の治療支援にむけた計算力学フレームワーク構築（大阪大学

大谷智仁、

和田成生）
コイル塞栓術は脳動脈瘤の破裂を予防する血管内治療法の一つであり、カテーテルを通して金属製のコイルを瘤内へ留置し、瘤への流入血流を抑制することで破裂を未然に防ぐ。本研究では、コイ
ル留置による脳動脈瘤への治療機序を明らかにし、患者個別に適切なコイルや留置手順を臨床へと提供する治療支援を目指している。平成 28 年度には、患者個別の瘤形状するコイル挿入の数値シミ

51ュレーションにより、コイル留置後の瘤形状を構築し、瘤内の流体領域を Volume fraction (VOF)
法により直交格子中に投影して、有限差分法による血流の数値流体計算を行うフレームワークを構築した。壁面境界条件には BDI 法（Weymouth and Yue, J. Compt. Phys., 2011）を用い、Fractional
step 法に基づく分離解法により非圧縮性 Navier-Stokes 方程式と連続の式を連立して解く。本手法
は直交格子法を用いたことで数値流体計算の前処理を大幅に自動化することに成功し、これにより、多数の症例や様々なコイル分布形状を用いての大規模な血流解析の実施が現実的に可能となった。
本年度では、昨年度に構築したコイル塞栓時の瘤内血流を安定かつ高精度に解析するための数値
流体計算フレームワークを用い、コイルの局所的な充填率の違いが瘤内血流の滞留度合に与える影響を検討した（Otani et al., J. Biomech. Eng., 2018）。平成 28 年度に構築したフレームワークを用いて、実際の患者の CT 画像から構築した 2 例の脳動脈瘤形状（図 2.2.14）に対し、コイルの局所的な充填率の違いが瘤内血流の滞留度合に与える影響を調べた。入口境界に生理学的に妥当な流量を設定し、出口境界の圧力を一定とした。ここで計算の簡略化のため、入口境界の速度は一定とした。コイルの挿入条件を変えて挿入シミュレーションを行い、各症例について 50 例の異なるコイル
分布形状を得た。ここで挿入するコイルは全例について同一のものとし、コイルの充填率が臨床において適切な値の範囲 (>25%) であるよう設定した。全例についてそれぞれ血流解析を行い、コイ
ル留置前後における瘤全域の血流の運動エネルギー変化率をそれぞれ取得した。このとき、瘤と母血管との境界面から 0.5 mm ずつ区間を分割し、0-0.5、 0.5-1.0、 1.0-1.5 mm の各領域におけるコイルの充填率を局所的なコイル充填率として定義し、運動エネルギー変化率と比較した。
局所的なコイル充填率に対する瘤内全体における血流の運動エネルギー変化率を図 2.2.15 に示す。
瘤全域のコイル充填率を 27%とした場合でも、局所的なコイル充填率は 5%から 40%程度ばらついた。母血管からの距離が 0-0.5 mm の領域において、局所的なコイル充填率と、瘤全域の血流の運動
エネルギー変化率には有意な負の相関が生じ（図 2.2.15 上段）、この傾向は、瘤のネック部から中心部に近づくにつれて小さくなった（図 2.2.15 中段・下段）。この結果は、コイルが密に充填されている場合でも、ネック領域におけるコイルの留置が不十分であれば、瘤内の運動エネルギーが十分に減少せず、瘤全域の塞栓が達成されない可能性を示唆している。

図 2.2.14: 計算対象とした瘤形状

52

図 2.2.15: 瘤と母血管との境界面から 0-0.5 mm（上段）、0.5-1.0 mm（中段）、1.0-1.5 mm（下段）の区間における局所的コイル充填率と、コイル塞栓前後における瘤内全体の血流の運動
エネルギー変化との関係。
(4) 左心房内血流解析のための流体構造連成計算手法の構築（大阪大学

大谷智仁、和田成生）

本研究では、心原性脳梗塞の主要因である左心房内血栓症の患者別リスク評価にむけ、患者個々の状態を反映した心房内血流の数値流体計算手法の構築を目指している。左心房は自ら能動的に伸縮することで左心室の補助ポンプとして機能しており、血流計算にあたって左心房形状だけでなく動態を考慮した患者個別解析が求められる。壁面の移動境界問題を安定的に扱うための計算手法の構築を目指し、本年度では脳血流解析用に構築してきた直交格子法による数値流体計算ソルバを拡張するとともに、CT 画像データに基づく実際の患者の左心房内血流解析への応用を試みた。
本研究にて提案する移動境界問題に対する直交格子法ソルバの構築にあたり、血流解析用の直交格子法ソルバで壁面境界の扱いとして用いている BDI 法における VOF 関数と固体領域の既知の速53度が時々刻々変化する場合を考える。代表者らの先行研究（Otani et al., Ann. Biomed. Eng., 2016）において、対象とする左心房の壁面は三角形要素の集合として与えられるとともに、壁面の速度場
は複数の正弦波の重ね合わせとして計算され、三角形要素の各節点に定義される。これらの情報を
直交格子内へ投影するアルゴリズムを構築した。まず、MPI 並列化 Fast Marching Method（Yangand Stern, J. Compt. Phys., 2017）を用いて、壁面の位置をゼロとする符号付距離関数（レベルセット関数）を直交格子を構成する全節点に設定した。この関数を変換して直交格子内に VOF 関数を構築した。さらに、三角形要素の各節点に定義された速度場を直交格子中に外挿するにあたり、符号付距離関数の空間勾配（壁の法線方向）を移流速度とし、壁面速度について移流方程式を解くことで固体領域内の速度場を設定した。この処理を CFD 計算の各時間ステップで繰り返し行った。
CFD 計算では非圧縮 N.S.式と連続の式を有限差分法により解いた。
テスト計算のため、実際の患者の CT 画像から 1 心拍中における左房の形状および速度場を取得した。肺静脈断面の境界条件として圧力一定（=0 Pa）を設定するとともに、僧帽弁の閉鎖時（左室収縮期）には僧帽弁領域を壁面境界と設定し、開放時には患者個別に計測した左房から左室への流出流量を与えた。
「京」を用いた 64 ノード×8 コア並列計算により、 3 心周期分の左心房内血流計算を行った。3 心周期目における左心房内血流のスナップショットを図 2.2.16 に示す。結果は生理学的知見と定性的に一致しており、時々刻々変化する左房の複雑形状および速度場を直交格子法により安定的に取り扱えることを示した。今後は構築した手法を応用し、左心房の画像解析から左房内血流解析までを直交格子中で実施するシームレスな解析システムの構築により、左心房内における血栓形成リスクの患者別評価にむけた大規模臨床データ解析を実現する。また、本手法は、サブ課
題 C で開発している心臓シミュレータと連立させることにより、患者個別の医用画像に基づく高度な診断支援を可能にするものである。

図 2.2.16: １心拍中における左心房内の血流動態のスナップショット

54

(5) 血栓形成シミュレータの開発（東京大学

清水和弥、高木周）

本研究は、心疾患・脳血管疾患の一因となっている血栓症を対象としたものである。血栓症を再現可能な数値シミュレーションツールを開発し、原因究明、疾患予測・予防、治療支援、薬効評価などに貢献することを最終的な目的としている。血栓形成は多くの因子が複雑に絡み合い連鎖的に起こる事象であるが、その機序は、血小板の粘着・活性化過程を支配要因とした血小板塊の形成過程と、凝固系の連鎖反応によるフィブリンの生成過程に大別される。血小板塊の形成過程は、血管内皮の損傷部位への血小板粘着を契機とするが、昨年度までの検討で、血小板粘着現象は流体力学的要因のみでおおむね説明可能あることを明らかにした。すなわち、血小板が血流内を受動的に運動（移動）する過程において、たまたま損傷部位に接触すると VWF-GPIb alpha 結合が形成され血小板が粘着する。血小板が損傷部位に到達するかどうかは周囲流体の速度分布に依存する。赤血球の存在は血流を乱す要因となるため、血小板粘着数はヘマトクリットに依存し、高ヘマトクリット
条件ほど血小板粘着は起こりやすい。このような現象はこれまでにも示唆されていたが、実験的な検証では生化学的作用と流体力学的作用とを切り分けることが難しく、数値シミュレーションによる検証を通じて初めて明確となった。
今年度の取り組みとして、血小板の粘着力に着目した検討を行った。血小板粘着は VWF-GPIb
alpha 結合によって実現・維持されている。血小板表面上には血小板 1 個あたり 15,000 個程度の
GPIb alpha が存在することが知られており、実験的な測定および分子動力学計算などから 1 組のVWF-GPIb alpha 結合の力は 50 pN 前後であることが示されている[1][2][3]。一方で、平均的な血流速度から血管内壁近傍の血流速度分布を推定すると、個々の血小板は概算で 100 pN 前後の流体力
学的な力を受けていると見積もられる。これらのことは、 VWF-GPIb alpha 結合が数組程度あれば血小板粘着を維持出来ることを示唆するものであるが、血小板表面上に存在する GPIb alpha 数との間に大きな乖離があることが新たな疑問点として浮かぶ。そこで、まず粘着している血小板が受ける流体力学的な力を数値シミュレーションから検証することとした。解析手法は、昨年度までと同様、流体構造連成解析手法[4][5]に改良を加え Immersed Boundary 法[6]を用いた流体構造連成
解析手法とした。構造物となる赤血球・血小板は超弾性体として扱っている[7][8]。さらに、血小板粘着を取り扱うために必要な VWF-GPIb alpha の結合に関する力学モデルを導入している。この力学モデルでは、遷移状態理論に基づいた反応速度定数[9][10]から確率的に各時間ステップにおける結合・解離を判定する。プログラムは OpenMP - MPI ハイブリッド並列プログラムとなっており、MPI 並列には V-Sphere ライブラリ[11]を利用している。結果の一例として、図 2.2.17 に粘着血小板の状態（瞬時値）とそのときの個々の血小板における VWF-GPIb alpha 結合の総和を示す。また、
図 2.2.18 に粘着血小板 1 個当たりの VWF-GPIb alpha 結合の数とそれらの結合が発揮している結合力の総和を示す。粘着血小板の位置はほぼ動かないので、結合力は血小板が流れから受ける力と等しいとみなすことができる。
血小板 1 個当たりの結合力を見ると、50 ～ 250 pN、平均で 130 pN 程度であり、おおむね事前の見積と合致する結果となっている。このときに血小板 1 個当たりに形成されている結合数は 50～ 350、平均で 200 程度となった。先に述べたように結合 1 個当たりの結合力が 50 pN 程度であることを考慮すると結合数が非常に多いと思われる。流体から受ける力が数百 pN 程度であることを合わせて考えると、結合のうち大部分は実際には力を発揮していないことを示している。この55

ことは、現在用いている数理モデルにおいて解離の速度定数が小さ過ぎる可能性を示唆するものであるが、一方で、実際の血流中においても本計算結果が示すように、VWF-GPIb alpha 結合のうち大部分は結合を形成してはいるものの有効には機能していない可能性も示唆している。VWF-GPIb
alpha 結合は血小板活性化の契機となる重要な現象であるため、血小板活性化機序を明らかにするためにもこの点に関する調査は重要である。実際の血流中における血小板粘着時の VWF-GPIb
alpha 結合の挙動がどのようなものか詳細が明らかでないために現時点で結論付けることはできな
いものの、数理モデルの改良も含めて今後検証しなければならない点であることを示している。

図 2.2.17: 粘着血小板の結合力検証解析の結果（𝛾1500s ）。粘着血小板の様子（瞬時状態、右上）およびその血小板における VWF-GPIb alpha 結合力の総和（左下）。

図 2.2.18: 粘着血小板 1 個当たりの VWF-GPIb alpha 結合数（左）と結合力（右）血小板粘着・活性化に関する解析と合わせて、凝固系連鎖反応に関する研究も開始した。よく知られているように、血栓形成は血小板の粘着・活性化過程と凝固系連鎖反応によるトロンビンおよびフィブリン生成過程に大別され、薬剤による対処療法も抗血小板剤を用いた血小板粘着・活性化を阻害する方法と、抗凝固剤を用いたトロンビン・フィブリン生成を阻害する方法とに大別される。
しかし血栓形成においては、状況によって起こり易い難いはあるものの、両者は同時に起こり進行して行くものである。したがって、本来どちらか一方だけを考慮しただけでは血栓症への対処としては不十分である。凝固系の連鎖反応は 12 種のタンパク・酵素などの因子が連鎖的に反応する（図2.2.19）。

56

図 2.2.19: 凝固系連鎖反応の概要図（文献 [12][13][14] などを基に作成）
経路は複数あるものの最終的な生成物はトロンビンおよびトロンビンの働きによって生成されるフィブリンであり、血栓成長においてはこれらの濃度が重要となる。この連鎖反応は、血管内皮の損
傷部位や粘着血小板の極近傍および表面において爆発的に進行する。フィブリンはフィブリンモノ
マーが自発的に多数重合することで形成される繊維状のタンパク質であり、大きなフィブリンポリ
マーが形成されるとそれが網の役割を果たし粘着血小板塊を覆うことで大きく強固な血栓へと成長する。したがって、凝固系の解析モデルとして、特にフィブリンモノマーが重合しポリマーへと成長する過程に着目し、フィブリン重合過程とそれへの寄与が大きな反応過程（膜表面上でのプロト
ロンビンからトロンビンへの活性化やフィブリノーゲンからフィブリンポリマーの生成過程など）を考慮した数値モデルとした。解析は dCME (discrete Chemical Master Equation) [15]で記述される各酵素・タンパクなどの反応過程を時間発展的に解く方法を採用した。凝固系反応に関する研究
は、化学反応を扱う際に用いられる手法を応用した決定論的な反応モデルに基づいた解析例は見られる[16]ものの、確率的な表現に基づいた解析例は少ない。生体内で起こる現象は本来確率的な要素を含むため、本研究では確率的な扱いが可能な数理モデルを採用した。結果の一例を図 2.2.20 に示す。解析はフィブリンモノマー n 個が存在する状態を初期状態とし、フィブリン n-mer が生成されるまでを dCME モデルに基づいて確率論的に解いている（図 2.2.20 は n = 64 の場合）。フィブ
リン生成は結合過程のみであり、解離過程は存在しない。さらに、フィブリンの重合度が大きくなるほど結合サイトが増えるため、結合速度は速くなる。一方で、フィブリンの分解はプラスミンと呼ばれる酵素が担うがその効果はそれほど大きくない。したがって、フィブリノーゲンからフィブ
リンモノマーが生成され血管内壁の損傷部位周辺に多数存在すると、移流などによってその空間から離脱するか、抗凝固剤などを外部から加えない限りフィブリンポリマーが生成される。図 2.2.20
はそのような結果を表しており、本解析モデルでは、0.7 s 程度経過すると 64-mer が形成される可能性が高いことを示している。今後、これまで開発してきた赤血球・血小板を含む大規模血流解析
モデルとの融合を図ることで、移流や薬剤などによる阻害効果を定量的により評価し、フィブリン
ひいては血栓が成長するかどうかを検討可能な解析モデルへと発展することが期待される。

57

図 2.2.20: dCME モデルを用いたフィブリン重合過程解析結果の例
参考文献[1] Kim, J., et al., Nature 466(2010), 992-997.
[2] Tobimatsu, H., et al., J Atheroscler Thromb, 22(2015), 1091-1099.
[3] Shiozaki, S., et al., J Atheroscler Thromb, 23(2016), 455-464.
[4] Sugiyama, K. et al., J Comput Phys, 230(2011), 596-627.
[5] Ii, S. et al., J Biomech Sci Eng, 7(2012), 72-83.
[6] C. S. Peskin, J Comput Phys, 10(1972), 252-271.
[7] R. Skalak, et al., Biophys J, 13(1973), 245-264.
[8] Modeling and Simulations of Capsules and Biological Cells, Edited by C. Pozrikidis, Chapman &Hall/CRC, 2001.
[9] S. Shiozaki, et al., J Biomech Sci Eng, 7(2012), 275-283.
[10] J. F. Marko et al., Macromolecules, 28(1995), 8579-8770.
[11] Ono, K., et al., Lecture Notes in Comput. Sci. Engrg., 67(2009), 9-18.
[12] Davie, E.W. and Ratnoff O.D., Science, 145(1964), 1310-1312.
[13] Hoffman, M. and Monroe, D.M., Hematol Oncol Clin N Am, 21(2007), 1-11.
[14] Goto, S., et al., Circulation Journal, 79(2015), 1871-1881.
[15] Cao, Y., and Liang, J., BMC Systems Biology, 2(2008).
[16] Guy, R.D., Mathematical Medicine and Biology, 24(2007), 111-130.
(6) 微小循環遊走細胞解析プラットフォーム（大阪大学松田秀雄）

毛細血管等の微小循環系での細胞の動態を非侵襲で長時間観察できる蛍光顕微鏡で得られる膨大な時系列画像データを高精度かつ高速に処理するための、大規模ライブセルイメージング解析技術を「京」を使って開発している。
サブ課題 B で取り組んでいる全脳循環代謝シミュレータでは、脳血管障害や脳機能障害の解析プ
ラットフォームの構築を目指しているが、主要な脳血管障害の一つである脳血栓症の発症メカニズ

58

ムの解明では、慢性的に進行する動脈硬化とそれに伴う血栓の形成の過程のシミュレーションが重要と考えられる。
そこで、本研究では、前年度に引き続き、生体内部を精密かつ経時的に観察できる二光子励起顕微鏡を使って、微小環境でのライブイメージングデータを取得し、動脈硬化の進行に伴う血管内部の変化をイメージング解析により解明することで、脳血栓のシミュレーションのデータ同化に必要な計測データを取得することを目指している。
近年、動脈硬化は炎症性病態であることが明らかとなりつつある（P. Duewell et al., Nature, 2010など）ことから、動脈硬化の進行を炎症を誘発する刺激剤（LPS、リポ多糖）の投与（以下、炎症刺激と略す）による血管内での炎症進行の経過でモデル化している。炎症刺激後の血管内皮の状態
変化を、緑色蛍光で標識した白血球（炎症時に血管内皮の細胞に接着・浸潤することが知られている）の動きの変化として観察している。
これまでに取得したライブイメージングデータは、マウス生体の毛細血管や骨髄腔等の微小環境中を二光子励起顕微鏡で撮影した経時観察動画像であり、毎秒最大 30 フレーム、各フレーム 512×512～2048×2048 ピクセル、24 ビットカラーであり、1 回の観察ごとに 60～2000 フレームの画像を取得した。生体の部位や条件を変えて 200 回以上撮像している。
平成 29 年度には新たに脂肪組織内での免疫細胞の動きと、骨組織内での領域識別について、イメー
ジングデータを取得して解析したので報告する。
脂肪組織のイメージング解析の意義としては、体重の増加により肥満状態が進むと、やがてメタ
ボリックシンドロームに陥るが、これは体内で起こる現象が問題と考えられている。肥満状態が進むと、体を守るはずの免疫細胞が“暴走”し始めると言われている。
「免疫細胞の暴走」とは、メタ
ボの人の脂肪細胞がある種の信号伝達物質を放出し、それを受け取った免疫細胞が活性化して、自らも同じ物質を放出することで免疫細胞の活性化が拡がることだと考えられている。
実際に、マウスに高脂肪食を与えたときに脂肪組織内での免疫細胞（マクロファージ）の活動を
二光子励起顕微鏡で観察した例を図 2.2.21 に示す。図 2.2.21 は観察した画像のフレームの一つで、青い領域が脂肪細胞で、緑色に光るのが免疫細胞、赤い筋が微小な血管を表している。図 2.2.22 (a)および(b)は、図 2.2.21 の黄色い四角で囲まれた領域を、隣接したフレームからそれぞれ抽出した画像である。図 2.2.22 (a)と(b)を見ると、細胞の形状が丸くなく、非常に伸びた形をしているが、これは脂肪細胞と脂肪細胞の間の隙間にはまりこんでいると思われる（特に大きな変形をしている細胞を黄色で囲っている）。高脂肪食ではなく通常の食事を与えたマウスでは、免疫細胞は動かずにじっとしているが、高脂肪食を与えてしばらくすると活発に動き出し複雑に形状を変えていく。
図 2.2.22(a)と(b)に対して深層マッチングを適用した結果を図 2.2.23 に示す。深層マッチングにより、細胞の動きだけでなく、形状の変化も変位ベクトル（白い矢印）の大きさとして検出されていることがわかる。このように脂肪組織中の免疫細胞の変化を二光子励起顕微鏡で観察し、深層マッ
チングを適用することで、脂肪組織での免疫細胞の活性化のレベル、すなわちメタボの進行度合い
が定量的に求められることが期待される。

59

図 2.2.21: 脂肪組織内での免疫細胞の観察画像（画像提供：大阪大学医学系研究科 石井優教授、水野紘樹助教）

図 2.2.22: 図 2.2.21 の一部の区画（黄色い四角）部分での細胞の形状変化（特に変化が激しい部分
を黄色い枠で囲っている）

60

図 2.2.23: 図 2.2.22 の画像に対する深層マッチングの適用結果「京」での深層マッチングの計算では、プログラムは C++言語でコーディングされており、それを「京」に移植している。畳み込み演算の部分はスレッド並列版の BLAS（富士通社の数値計算ラ
イブラリ中の SSL2BLAMP）を使用し、さらに OpenMP によりスレッド並列を実装している。また、東京工業大学の秋山研究室で開発された MPIDP を使用して、経時観察した動画像を 100 フレ
ーム程度のブロックに分割して、各ブロックの深層マッチングを MPI によるノード並列で処理するようにしており、スレッド並列とノード並列を組み合わせたハイブリッド並列処理となっている。
「京」での実行時間は、1 ノード 8 スレッド並列で実行したとき、512×512 画素の画像で隣接 2 フレーム間の 1 回のマッチングを計算するのに約 13.6 秒かかっている。動画像のフレームを「京」の各ノードに分配して MPI によるノード並列処理を実装して、フレーム数を増やしていくことで「京」の全ノードを使った並列処理を達成できた。1 ノード当たりマッチングを 100 回実行するように固定した Weak Scaling で、「京」全ノードを使った並列処理では 約 829 万フレームという膨大な量のマッチングを約 41 分で実行できた。マッチング１回あたりの処理時間は約 24.6 秒となっている。
利用している二光子励起顕微鏡では 30 フレーム/秒で画像を取得できるので、「京」の約 738 ノード（=30×24.6）を使えば、画像の取得速度と同程度の処理速度で深層マッチングの計算ができることなる。実際には、画像ファイルの「京」への転送に時間がかかるため、計算前に「京」に移しておく必要がある。免疫細胞の"暴走"は、病態としては炎症という形で現れるが、炎症の進行は急性で 3～4 時間以上持続するため、フレーム数にすると急性で 32.4 万～42.3 万フレームの計算が必要である。さらに慢性の炎症は 1 週間以上かけて徐々に進行するとされている。常に 30 フレーム/秒で観察する必要はないが、一定の間隔で週単位の時間をかけて画像を取得する必要があり、膨大な量のフレームを解析する必要がある。また、本稿では 2 次元の画像しか示していないが、二光子励起顕微鏡は XY の 2 次元だけでなく深さ方向の Z 軸で焦点を設定できるため、得られる画像は 3 次元

61

成熟
骨組織
血管

骨髄腔

図 2.2.24（上）: 骨組織の観察画像の例（提供 大阪大学医学系研究科 石井優教授）
図 2.2.25（下）: 時刻経過による骨髄腔領域への血液の染み出し（左から右に時刻経過）＋時間軸の 4 次元画像となる。3 次元画像では 2 次元画像と比べて膨大な数のフレームが生成される
ため、より大きな計算資源が必要となる。さらに、１例だけの解析ではなく、将来的に数千から数万例のコホートで処理することも考えられる。以上を合わせて考えると、ポスト「京」規模の計算資源の利用が必要になると考えられる。今後、intravital imaging の画像データのスパコンによる解析が進んでいくのではないかと期待される。
次に骨組織内での領域識別の例を示す。骨粗鬆症や関節リウマチなどの疾病には骨の代謝に関わる破骨細胞や骨芽細胞の働きが関与していることが知られているが、そのメカニズムは完全に解明されていない。硬い石灰質に覆われた骨内部を生きたまま観測することは困難であったが、二光子励起顕微鏡の登場など、顕微鏡技術の進展により生きたままでのイメージングが可能となった。しかしながら、大規模なデータに対して手作業で解析を行うことは困難であるため、計算機による解析手法が求められる。しかし、二光子励起顕微鏡は、蛍光物質を介して対象を間接的に観測する方式であるため、複雑な撮影環境においては複数の組織が同一の色で観測される、境界が不明瞭になる、染色できない組織は観測できないなどの問題があるため、単純な画像処理技術を用いるだけでは解析は困難である。
そこで、骨組織のライブイメージング画像に対し、骨髄腔領域を領域分割する画像解析手法を開発した（本手法は、骨組織以外の同様の性質を持つ二光子励起顕微鏡画像への応用も可能である）。
図 2.2.24 はマウス頭蓋部の骨組織内の血管に緑色の蛍光物質 FITC (Fluorescein Isothiocyanate)を投与して二光子励起顕微鏡で撮影したものであり、血管が緑色で、成熟骨組織が二次高調波により青色で観測されている。骨髄腔の領域は染色できないため、直接観測することは不可能である。
しかしながら、血液が骨髄領域に染み出すため、骨髄領域の色は時間とともに緑色に染まっていく
ことから骨髄腔領域を判断できる（図 2.2.25）。また、骨髄腔が直接確認できないだけでなく、血液中を流れる比較的大きい組織により血管領域においても暗い領域があること、血管領域においても色ムラが存在するなど、非常に複雑な画像となっている。骨髄腔は骨の代謝のメカニズム解明に非62常に重要な役割を果たすことが知られているため、骨髄腔領域を領域分割の対象とした。
図 2.2.24 に示す骨組織の画像を対象として、骨の代謝に重要な役割を果たしていると考えられている骨髄腔領域を領域分割する画像解析手法を開発した。骨の代謝のメカニズムを明らかにすることにより、骨粗鬆症や関節リウマチなどの疾病の解明に繋がることが期待できる。また本手法は、骨組織以外の同様の性質を持つ二光子励起顕微鏡画像への応用も可能である。

データセット

1

2

3

4

用

真値
（手動で作成）

用

真値を いてテクスチャの学習

分類器を いてテクスチャ解析

グラフカット（テクスチャ解析結果と空間連続性を考慮）
領域分割
結果出力

分類器
テクスチャのみによる
ピクセルごとの認識結果

図 2.2.26 開発した領域分割手法の概要
図 2.2.26 に本手法の概要を示す。骨髄腔領域と血管領域が同じ緑色で観測されること、また色ムラなどの問題により、輝度情報を中心とする方法による解析は困難である。そこで本研究では、血管中を流れる際の画像の模様が組織ごとに異なることに着目し、テクスチャ解析をベースとした解析手法を開発した。各観測画像について手作業により作成された骨髄領域の画像（真値）を用意し、骨髄領域とそうでない領域のテクスチャ特徴量を、機械学習手法の一つである SVM (Support Vector
Machine)を用いて学習する。これにより、学習結果として出力される分類器を用いて未知のデータに対してテクスチャ解析を行うことが可能になる。テクスチャ解析結果はあくまで模様の観点から判別を行うもので、輝度値そのものや組織の空間の連続性などは一切考慮しない。

データセット 1

データセット 2

データセット 3

データセット 4

図 2.2.27: 撮影を行った骨髄のデータセット（スケールバー : 50 μm）

63そこで、テクスチャ解析結果を、空間の連続性などを考慮して領域分割を行うことができるグラフ
カット法に統合する。グラフカット法は画像解析分野で広く用いられている手法で、与えられたモデルに対して最適解を出力することを保証するだけでなく、拡張性が高く、テクスチャ解析結果を統合するようなモデルを構築する場合でも応用が容易である。
テクスチャ解析には広く用いられている 2 次元ウェーブレット変換を利用し、ウェーブレット係数を特徴量として解析を行う。ウェーブレット変換を行う際には、マザーウェーブレットやスケー
ルパラメータと呼ばれるパラメータを指定する必要がある。また、SVM を用いる際にも学習が必要なパラメータがある。そこで本研究では、これらのパラメータの探索も同時に行い、得られた中で最も高い精度を達成するパラメータを用いて分類器の出力を行った。
また、グラフカット法では画素ごとに骨髄腔らしさのモデル化を行うデータ項と、隣り合う画素
同士の空間の連続性を考慮する平滑化項の和で与えられるエネルギー関数を定義し、最小化することで領域分割結果を得る。グラフカット法のデータ項 𝑔 𝑋
𝑔 𝑋

0

𝑒

𝜆𝑔 𝑣

𝑔 𝑋

1

𝑒

𝜆 1

ただし、𝑔 𝑣
ここで、𝑔 𝑋1を次に示す。
(2.2.1)𝑔 𝑣（骨髄腔）

0 （それ以外の領域）(2.2.2)(2.2.3)0 は画素 𝑣 を背景とした時に発生するエネルギーであり、𝑔 𝑋領域とした時に発生するエネルギーである。𝐼 は画素 𝑣 の輝度値 (0

𝐼

1 は骨髄

1) であり、𝜌 は画素の輝度値の変化に対してデータ項を調整するパラメータである。𝜆 はテクスチャ解析の結果をどの程度重視するかを決定するパラメータである。一般的なグラフカット法では輝度のみを考慮するのに対し、本手法ではテクスチャ解析結果の項 𝜆𝑔 𝑣 および 𝜆 1

𝑔 𝑣を加える拡張を行った。

なお、平滑化項については従来法と同様のものを用いている。
データは二光子励起顕微鏡を用いて撮影を行ったものを使用する。通常、二光子励起顕微鏡は 800〜1,000𝜇𝑚 程度まで到達が可能であるが、骨組織を撮影する場合には骨基質に含まれるリン酸カル
シウムにより励起光を散乱させるため、200𝜇𝑚 程度の深度しか到達できない。そのため、骨基質が比較的薄く、骨髄腔までの距離が近いマウスの頭頂骨の撮影を行う。データは 2 次元の時系列画像
として取得し、骨髄腔の領域分割行うことが目的であるため緑チャンネルのみを抜き出して解析を行う。機械学習により多様な画像の解析を行うため、画像の解像度、データのフレーム数、骨髄領域の数、明るさ、スケールなどが異なる 4 種類のデータを用意し、それぞれについて学習および評価を行うための真値データを手作業により作成した。データセットの概要を図 2.2.27 と表 2.2.1 に示す。評価の際には 3 種類のデータセットを用いて学習を行い、残る一つのデータセットに対して
領域分割結果を求める。従来法を一般的なグラフカット法として、開発した領域分割法（提案手法）と比較したものを図 2.2.28 に示す。また、定量評価として認識結果の F 値（再現率と適合率の調和平均）で比較した（表 2.2.2）。その結果、本手法が概ね良い結果が得られていることを確認した。
今後の課題として、本解析技術を用いて骨髄領域の生物学的な働きについて検証していく必要がある。

64表 2.2.1: データセットの詳細
データ

解像度

セット

フレー

骨髄腔

ム数

の数

1

752×7527

6

2

1,059×1,000

15

6

3

958×1,000

22

4

4

1,001×1,0001312

表 2.2.2: F 値による定量評価
データ

従来法

提案手法

1

0.227

0.219

2

0.608

0.615

3

0.655

0.745

4

0.660

0.728

セット

データセット

従来法

提案手法

図 2.2.28: 本領域分割技術（提案手法）と従来法の比較

2.3 個別生体機能評価のためのデータ同化シミュレータの開発（大阪大学

野崎和徳、吉永司、玉川裕夫、

和田成生）(1) 発話障害の診断および治療を支援するための大規模生体シミュレータの開発
A) はじめに
本年度は、発話障害の診断の根拠となり得る物理的メカニズムに関して新たな知見を明らかにした。特に摩擦音の発生メカニズムについて、異なる二つの摩擦音/s/や/sh/を空力音響学的に解き明かすことにより、それぞれの音素を構成するために必要な要素の推定に成功した。これらの成果は、音響分野で最もインパクトファクタ－が高い国際誌である Journal of AcousticalSociety of America に掲載され、さらに流体分野で最もインパクトファクタ－の高い国際誌の一つである Physics of Fluids に採録された。以下ではこれらの業績に関して主となる内容に関して述べる。
B) 大規模空力音響シミュレーションを発話障害に適用する意義
言語聴覚士や歯科医師による治療において、口腔形状の特性や先天性疾患等による舌機能の障害のため、摩擦音の構音（註１）に必要な狭窄流路を適切に形成できない[1-2]ことが、課題となっている。そのため、発音時の口腔内でどのように気流から音が発生しているのかをシミュ
レーションを用いて明らかにすることにより、発話障害に関する言語聴覚士が行うリハビリの支援等に活用することが期待できる。
摩擦音は、日本語のサ行を発音する際に用いられる/s/や/sh/等の音である。この摩擦音を発音する際には、口腔の前方部で上顎と舌先端により狭窄流路を形成し、狭窄流路から発生したジェット流から音が発生することが知られている[3]。これまで、摩擦音の発生メカニズムを調べ65

るため、様々な手法が用いられてきた。Shadle[4]は、口腔形状を円管にするという単純化したモデルを提案し、狭窄流路の形状や位置が発生音の特性に与える影響を調べた。Howe ら[5]は、口腔内の特徴的な形状を表す４点の断面積と、前後方向の長さから 1 次元の解析モデルを提案
し、発生音のスペクトルを予測した。しかし、これらのモデルでは、口腔内にどのような流れが発生し、音を発生しているのかについて議論されていない。
発音時の口腔内の気流については、これまで流体力学シミュレーションにより調べられてきた。Nozaki[6]は、CT 画像から取得した/s/発音時の口腔形状に対して、非圧縮性流体のラージ
エディーシミュレーション (Large Eddy Simulation， LES) を行うことにより、ジェット流の乱流から発生する音源を可視化した。しかし、非圧縮性の仮定を入れてシミュレーションを行ったため、音の伝播に関しては考慮されていない。口腔内のジェット流による乱流と音源の発生、またその音源からの音の伝播を全て考慮してシミュレーションを行うには、一般的な流体シミュレーションに比べて非常に計算コストが大きくなる。そのため、大規模シミュレーシ
ョンを行うことは必須である。
本報告ではまず、CT 画像から取得した口腔形状（実形状）に対して圧縮性流体の LES を行うことで、ジェット流の発生と音の伝播を同時に計算し、気流と音の可視化により発音のメカ
ニズムに関して述べ、次に、解剖学的な形状のうち、空力音響学的に摩擦音発生に必要な要素を推定するため、単純化した口腔形状を用いて行った実験とシミュレーションの結果を報告する。
C) 研究手法の詳細
まず、摩擦音/s/及び/sh/発音時の口腔形状を、CT画像より抽出した。CT画像の被験者は、発話障害のない42歳男性（被検者A）の日本語話者である。摩擦音/s/はコーンビームCT (CBMercuRay，日立メディコテクノロジー)、/sh/は320-row Area Detector CT (Aquilion one，東芝メディカル) により撮像し、口腔領域を含む画像を取得した。この時、/s/は9.6秒間持続して発音している間に撮像し、/sh/は「うすいみそしる」と発音する際の1コマ(0.05秒間)として撮像した。そのため、/s/の画像の解像度は0.1 × 0.1 × 0.1 mm ボクセル、/sh/の解像度は0.488 × 0.488
× 0.5 mmボクセルとなっている。この画像よりソフトウェアitk-SNAP[7]を用いて輝度値が空気である領域を抽出し、流路を形成した。また、計算コストの削減のため、咽頭部 (口唇から後方に約6.4 cm) より後方の流路は計算領域から除いた。次に、流路内にて流体解析を行うため、計算格子を構築した。摩擦音/s/及び/sh/発音時の口腔流路壁面の計算格子を図2.3.1に示す。ジェット流による乱流の渦を解像するため、狭窄から口唇にかけて計算格子が細かくなるよう設定した。/sh/の格子表面に凹凸が見られるのは、CT画像の解像度が低いためである。格子数は/s/で約2.3×107、/sh/で約3.5×107となっている。
次に、本研究で作製した口腔の単純形状モデルを図 2.3.2 に示す。口腔形状を矩形管で表現し、寸法は/s/を発音する被験者 A の医用画像[8]から決定した。外管は被験者の口唇の大きさ(8 × 25mm)から決定し、厚さ 1 mm 高さ 4 mm の上下前歯をモデル出口から 10 mm の位置に設置した。上下前歯の隙間は 1 mm となっている。管内には、高さ 1.3 mm 幅 8 mm の狭窄を持つ舌を設置した。舌先端と下前歯の距離を LCA とし、/s/と/sh/の口腔形状の違いを表すため LCA= 066

mm 及び 7 mm となるよう舌を位置させた。モデル出口には顔の面を摸したバッフルを設置した。シミュレーションにおける格子数は約 4.0×107 となっている。

図 2.3.1: Computational grids of realistic vocal tract geometry of /s/ (a) and /sh/ (b).

図 2.3.2: Simplified vocal tract model of sibilant fricatives /s/ and /sh/.
発音時の口腔内に発生する流れ場と音場を計算するため、圧縮性流体の3次元Navier-Stokes方程式と状態方程式
𝜕𝜌̅
𝜕𝑡

𝜕𝜑
，
𝜕𝑥

𝜕𝜑
𝜕𝑡

𝜕𝜑 𝑢
𝜕𝑥

𝜕𝑝̅
𝜕𝑥

𝜕𝜎
，
𝜕𝑥

𝜕𝜌̅ 𝐸
𝜕𝑡

𝜕𝜑 𝐸
𝜕𝑥

𝜕𝑝̅ 𝑢
𝜕𝑥

𝜕𝜎 𝑢
𝜕𝑥

(2.3.1)

𝜕𝑞
，
𝜕𝑥
𝑝̅

𝛾

1 𝜌̅ 𝑒̅，𝑒̅𝐶𝑇(2.3.2)(2.3.3)(2.3.4)を離散的に解いた。ここで、𝑡は時刻、𝜌̅ は密度、𝑢 は𝑥 軸(i = 1， 2， 3)方向の速度、𝜑
流束、𝑝̅ は圧力、𝐸

1⁄2 |𝑢 |𝜌𝑢 は𝑒̅ は全エネルギー、𝑒̅ は内部エネルギー、𝑇は温度である。𝜎 及67

び𝑞 は粘性応力テンソルと熱流束であり、
𝜎

2𝜌̅ 𝜈

𝜈

𝑠
1
𝛿 𝑠
3

𝑞
と計算する。この時、動粘度𝜈

𝜌̅ 𝛾𝐶 𝛼𝜕𝑇𝜕𝑥

𝛼と熱拡散率𝛼(2.3.5)(2.3.6)

はLESのサブグリッドスケール (Sub Grid

Scale， SGS) での渦粘性と熱拡散を表しており、本研究では、Fureby[9]の一方程式モデルを用いて決定した。𝑠 は変形速度テンソルであり、𝑠1⁄2 𝜕𝑢 ⁄𝜕𝑥𝜕𝑢 ⁄𝜕𝑥 と計算する。また、気体に関する定数である動粘度𝜈、比熱𝐶 、比熱比𝛾、熱拡散率𝛼は、摂氏20度の空気の値を用いた。以上の式を有限体積法数値計算ソフトウェアOpenFOAM 2.3.1を用いて口腔流路内の格子上で計算を行った。計算の境界条件としては、入口部に一様流速、出口部に無反射境界、口腔壁面に断熱ノンスリップ条件を課した。入口流量は被験者の呼気流量である400 cm3/sとし、単純形状では300 cm3/sとした。
「京」を用いた大規模解析において、各口腔モデルの格子を使用するノード数の数に分割し、MPIによる並列計算を行った。本研究では、各ケースにおいて385ノード並列で約81時間計算を行った。
単純モデルを用いた物理実験では、モデルをアクリル板により作製し、コンプレッサ(YC-4RS，Yaezaki)から空気を流入することにより音を発生させた。発生した音はマイクロフォン(Type4939， Bruel & Kjaer)により遠方 30 cm（地点 A）においてサンプル周波数 44．1 kHz で計測した。また、被験者 A の発音する/s/及び/sh/を同じ地点 A で計測した。計測した音は離散フーリエ変換により周波数スペクトルに変換した。
D) 結果と考察
実形状を用いたシミュレーションにより得られた結果の精度を検証するため、遠方場での圧力𝑝̅と被験者の発音の計測結果を比較した。図2.3.3に音圧スペクトルのシミュレーション結果と、被験者の摩擦音のスペクトルの実測値を示す。この時、シミュレーションで得られた圧力及び実測値の圧力は、基準値p0 = 20×10-6 Paを用いて音圧レベルとした。また、シミュレーション
では口唇から7 cmでの音圧、被験者の発音はマイクロフォンで計測した口唇から9 cmでの音圧を示しており、被験者の発音は15回の計測の標準偏差をエラーバーで示している。/s/及び/sh/のシミュレーションともに、摩擦音の特徴である高周波数域のブロードバンドノイズを再現した。また、/s/と/sh/の違いである、4 kHzから5 kHzのピーク周波数の違いを再現した。シミュ
レーション結果において、/s/の9 kHz及び/sh/の6.5 kHzに被験者の発音には見られないピークが発生した。このようなピーク形状は、周期的な音として聞こえる可能性がある。この原因としては狭窄等の狭めにおいてCT画像から流路を抽出する際に壁面が実際の形状とずれたため、実際の口腔内には発生しない周期的な渦が発生した影響と考えられる。
シミュレーションにより得られた速度場の瞬時値を図2.3.4に示す。流速|𝑢 |の値は流路内の368次元の計算格子点上に得られるため、流速が1 m/s以上の領域をボリュームレンダリングにより表示した。また、口腔壁面を茶色で表示している。/s/の口腔形状では、狭窄流路から発生したジェット流が上下前歯で大きく乱れる様子がわかる。/sh/の形状では、舌と下前歯の空間でジェット流が乱れ始め、/s/の形状に比べ、横方向に広く流れが進んでいく様子がわかる。図2.3.5に乱流の渦構造を可視化した様子を示す。この時、渦構造を表示するため、変形速度テンソル𝑠 の第２普遍量Q = 109の等値面を表示している。この渦構造によって、/s/では上下前歯で流れが乱れ始めているのに対し、/sh/では舌と下前歯の空間で流れが乱れている様子が確認された。これ
らの結果より、/s/と/sh/で舌の形状により流れの構造が変化した様子を可視化することができた。

図 2.3.3: Sound spectrum of /s/ (a) and /sh/ (b) predicted by the realistic vocal tract geometry.
図2.3.6にボリュームレンダリングした速度場と圧力𝑝̅ = 105 Paの等値面を示す。図内の圧力の等値面は複数存在するが、同一時刻の値であり、主に音波の波長の間隔で存在する。圧力の等値面は時間が進むに連れて遠方場(図右側)へ移動した。この可視化により、上下前歯の流れの乱れから音の粗密波が遠方へ伝播していく様子を可視化できた。

69

図2.3.4: Instantaneous flow velocity in the vocal tract geometry of /s/ (a) and /sh/ (b).

図2.3.5: Iso-surface of second invariant of velocity gradient tensor in the realistic vocal tract
geometry of /s/ (a) and /sh/ (b).

70

図2.3.6: Instantaneous flow velocity and iso-surfaces of pressure value p = 105 (Pa) in the
vocal tract of /s/ (a) and /sh/ (b).
単純形状を用いた実験により得られたスペクトルを図 2.3.7 に示す。図 2.3.7(a)にはモデル
LCA= 0 mm の音と被験者の/s/のスペクトルを比較している。モデルから発生した音は被験者 Aの/s/の 4.7 kHz のピークを再現しており、4.7 kHz 以上の高周波数域の音量の差は 5 dB 程度と非常に似た傾向を示している。図 2.3.7(b)にはモデル LCA= 7 mm の音と被験者 A の/sh/のスペクトルを比較している。モデルから発生した音のピークは 3.5 kHz となり、被験者 A の/sh/と同様に図 2.3.6(a)と比べてピーク周波数は約 1.2 kHz 低下したが、10 kHz 以上の高周波数域で10 dB 以上の差が見られた。モデル内の舌を 7 mm 後ろに移動することにより、被験者 A の/s/と/sh/の周波数特性の違いを表現することができたが、より/sh/のスペクトルを精度良く再現す
るには、舌の位置だけでなく、医用画像の観察により報告されている舌先端の形状の変化等もモデルで表現する必要があることが示唆された。
図 2.3.8 にシミュレーションにより得られた正中矢状面の速度分布を示す。LCA= 0 mm では、
狭窄から発生したジェット流は前歯に衝突することにより大きく乱れ、上下前歯の隙間を通過した。LCA= 7 mm では、狭窄から出た流れは前歯と舌の間の空間で大きく旋回し、舌先端付近
で大きく乱れている。LCA= 0 mm 及び 7 mm の発生音のピーク周波数(4.7 kHz 及び 3.5 kHz)における圧力振幅分布を図 2.3.9 に示す。流れに駆動された圧力振動は LCA= 0 mm で前歯付近に発生しているのに対し、LCA= 7 mm では前歯から舌先端まで前後方向に広く分布した。この
圧力分布の違いがモデル内の定在波の波長を変化させ、特徴となるピークの周波数を変化させたと考えられる。

71

図 2.3.7: Spectrum of sound measured at 30 cm from the model. The spectra were comparedwith subject’s /s/ (a) and /sh/ (b).

図 2.3.8: Velocity magnitudedistribution at mid-sagittal plane of the model. (a) LCA = 0 mm,
and (b) LCA = 7 mm.

図 2.3.9: Pressure amplitude distribution at mid-sagittal plane of the model. (a) LCA = 0 mm
at 4.7 kHz, and (b) LCA = 7 mm at 3.5 kHz.
E) まとめ本研究では、摩擦音/s/及び/sh/の発生メカニズムを調べるため、空力音響シミュレーションを行った。シミュレーションにより流れから音が伝播していく様子を確認でき、/s/と/sh/のジェット流の渦構造や音伝播の様子の違いを観察することができた。また、単純モデルを構築すること
により、/s/と/sh/の発音の区別には舌の前後の移動による流れと音場の共鳴を変化させることが必要であることが明らかとなった。今後はシミュレーションを、発音障害を持つ被験者の口腔形状に適応するだけでなく、結果に関して言語聴覚士と議論を行い、リハビリの支援の応用可能か議論を行う。
参考文献[1] C. Runte, M. Lawerino, D. Dirksen, F. Bollmann, A. Lamprecht-Dinnesen, and E. Seifert“The influence of maxillary central incisor position in complete dentures on /s/ sound production,”
The Journal of Prosthetic Dentistry, 85, 485-495, 2001.
[2] S. Wood, J. Wishart, W. Hardcastle, J. Cleland, and C. Timmins “The use ofelectroparatography (EPG) in the assessment and treatment of motor speech disorders inchildrenwith

Down’ssyndrome:

Evidence

from

two

case

studies,”Developmental

Neurorehabilitation, 12, 66-75, 2009.
[3] G. Fant “Acoustic theory of speech production,” Mouton, The Hague, Paris, 328 pp., 1960.

72

[4] C. H. Shadle, “The acoustics of fricative consonants,” Ph. D. thesis, Massachusetts Instituteof Technology, Cambridge, MA, 1985.
[5] M. Howe and R. McGowan, “Aeroacoustics of [s],” Proceedings of the Royal Society A, vol. 461
1005-1028, 2005.
[6] K. Nozaki “Numerical simulation of sibilant [s] using the real geometry of a human vocal
tract,” High Performance Computing on Vector Systems 2010, Springer, Berlin, Heidelberg,
137-148, 2010.
[7] http://www.itksnap.org[8] Nozaki, K., Yoshinaga, T., Wada, S. “Sibilant /s/ simulator based on computer tomography
images and dental casts,” Journal of Dental Research 93(2), (2014), 207-211.
[9] C. Fureby “On subgrid scale modeling in large eddy simulations of compressible fluid flow,”
Phys. Fluid., 8(5), 1301-1311, 1996.
註１

口から声を出す場合、すなわち口腔内で音の調整がなされることを「構音」、口腔だけでなく喉頭部声門付近でも音の調整がなされることを「発音」、「発音」に関わる筋骨格系を制御し、機能を実現することを「発話」と呼ぶ。
(2) 構音メカニズムの解明に向けた発話時の筋線維収縮力の逆推定手法の研究開発（大阪大学

伊井仁

志、吉永司、和田成生、野崎和徳、玉川裕夫）
舌は調音器官の一つであり、舌内に分布している筋線維を収縮させることにより舌を能動的に変形させ、口腔形状を変化させることで調音を行っている。筋線維は舌内部で細かく配向しており、発話の際はこのような筋線維の収縮を、主に舌下神経を介して大脳皮質中枢から制御することで、舌形状を機敏かつ複雑に変形させている。発話における舌変形メカニズムの調査は古くから行われているが、筋線維の複雑さのため詳細な解明には至っていない。これまで実験の観点から発話時の舌
変形メカニズムの調査が広く行われてきたが、内舌筋に関する記述は少ない。内舌筋は舌の内部に複雑かつ微細に配向しており、発話時はこのような複数の筋線維が収縮することで舌形状を変化させている。そのため、実験から十分なデータを計測することは困難である。この問題に取り組むため、本研究では、舌の 3 次元変形モデルと医用画像より得られる舌形状を一致させることで発話時の舌筋線維収縮力を逆推定する手法に関して研究開発を行ったので報告する。
阪大学歯学部付属病院より提供された、安静位の 1 名の男性日本語話者の 3D-CT 画像から舌形状を抽出することで、3 次元舌モデルを作製した。画像は大阪大学大学院歯学研究科・歯学部付属病院
倫理委員会で承認を受けた。ボリュームデータは、冠状スライス 320 枚、ピクセル数 512×512、ボクセルサイズ 0.488×0.488×0.5 mm である。画像処理ソフト Amira(Visage Imaging Inc.)を用いて形状の抽出を行い、解析ソフト ICEM CFD(ANSYS)により四面体一次要素で構成されるボリュー
ムメッシュを作製した。節点数は 49,043、要素数は 293,841 である。筋線維は、前部オトガイ舌筋(GGa)、中部オトガイ舌筋(GGm)、後部オトガイ舌筋(GGp)、舌骨舌筋(HG)、茎突舌筋(Sty)の五つの外舌筋と、上縦舌筋(IL)、下縦舌筋(SL)、横舌筋(Trans)、垂直舌筋(Vert)の四つの内舌筋の計 9 つに分けて配向させた。各筋線維の配置と配向方向は解剖学データ(Miyawaki, 1974, Ann Bull RILP,73

8, 23; Takemoto, 2001, J Speech Lang Hear Res, 44, 95)に基づいて決定し、各メッシュに対し筋線
維走行を定義した。茎突舌筋に関しては、CT 画像から茎突舌筋の外舌部の形状を抽出することができなかったため、茎突舌筋の内舌部の筋線維のみを配向させた。また、オトガイ舌骨筋や顎舌骨筋
といった口底筋は考慮していない。境界条件として舌モデルの底面に変位ゼロの固定条件を与えた。
これらの様子を図 2.3.10 に示す。

図 2.3.10: 医用画像から舌形状の再構築と舌筋線維配向の様子(GGa: 前部オトガイ舌筋、GGm: 中部オトガイ舌筋、GGp: 後部オトガイ舌筋、HG: 舌骨舌筋、Sty:
茎突舌筋、IL: 上縦舌筋、SL: 下縦舌筋、Trans: 横舌筋、Vert: 垂直舌筋)
初期配置で表記された静的力学平衡式を、一次要素を用いた Galerkin 有限要素法により離散化する。ここで、図 2.3.10 で定義している舌筋線維群にそれぞれ能動的に収縮応力を与えることで変形が駆動される。提案した舌変形モデルにより求められる舌形状と発話時に取得された舌形状に関する最小化問題を設定することで、筋線維応力を逆推定する。本研究で使用する舌形状に関する初期形状と目標形状(発話時の変形形状)は、それぞれ異なる画像データより得られるため、両者の対応関係 を一意 に決 定する こと はでき ない 。そこ で、 目標形 状に 関する 距離 関数 (Signed Distance
Function: SDF)を利用することで、モデルおよび計測より取得した舌表面の距離を最小にする最適化問題を考える。SDF は SDF 生成ライブラリ V-SDFlib (VCAD、RIKEN)より計算した。これにより任意の点において目標とする形状表面との最近傍距離を離散的に評価できる。定式化においては、
SDF の持つレベルセット関数という特徴より、変位に関する移流方程式をラグランジュ的に評価し、反復的に最適化問題を収束させる。
Stavness ら(2012, J Acoust Soc Am, 109, 2974)が推定した舌筋線維活性化を基にした値を収縮応力として与え順解析を行い舌の 3 次元形状を得て、これを逆推定の目標形状として利用し、推定された収縮応力と比較することで提案した逆解析手法の妥当性検証を行った。収縮応力は、舌先端の前方突出(forward protrusion)、上方向の曲げ(upward bending)をそれぞれ表現し与えている。図2.3.11 に、初期形状、逆解析によって得られた舌形状と目標形状および推定された筋線維収縮力を74示す。また、順解析で与えた収縮応力(reference)および推定された収縮応力(estimated)の比較を図2.3.12 で示す。本手法により非常に高い精度で収縮応力が推定されており、異なる複数の筋線維応力の組み合わせにおいても、高い推定精度を保っていることから、本手法のロバスト性を示すことができた。

図 2.3.11: 妥当性検証問題における、舌の初期形状(i)、舌先端の前方突出(ii)および上方向の曲げ(ii)に関する推定結果。カラーマップは筋線維収縮応力を示す。

図 2.3.12: 舌先端の前方突出(a)および上方向の曲げ(b)の解析条件における筋線維収縮応力の推定結果(estimated)および参照解(reference)との比較。
次に、実際の発話時に医用計測された舌形状を用いて解析を行った。二つの母音/u/を発話している舌形状を、大阪大学歯学部付属病院より提供された「うすいみそしる」と発話する被験者の 3D-CT
画像から抽出した。被験者は安静位の際と同じ 1 名の男性日本語話者である。CT 画像は全部で 76
フレームあり、フレーム間の時間ステップは 0.05 s である。各フレームにおけるボリュームデータ
は、冠状スライス 320 枚、ピクセル数 512×512、ボクセルサイズ 0.488×0.488×0.5 mm である。
この際、最も形状が大きく異なっていた、初めの母音/ʉ/(First /ʉ/)と最後の母音/ɯ/(Last /ɯ/)を発話している、5 フレーム目と 70 フレーム目のボリュームデータをそれぞれ使用した。推定結果に関して、正中矢状面、冠状面、横断面における舌の輪郭を図 2.3.13 に示す。初めの母音/ʉ/を目標形状とした場合では、逆解析によって得られた形状は各断面で目標形状と定性的に一致した。一方で、最後の母音/ɯ/を目標形状とした場合では、舌先端部や側面部において両者の形状に若干の差異が生じた。舌背と舌咽頭表面のみに着目すると、両者の形状は各断面で定性的に一致した。推定された筋75

線維応力は、どちらも負の応力が発生した。これは筋線維が伸張する方向に応力が発生したことを示している。本研究では安静位の舌を初期形状として使用したが、その時の舌は応力ゼロ状態ではなく、安静位の形状を保つために収縮応力が作用していたと考えられる。そのため、これらの負の応力は初期状態ですでに収縮していた筋線維が弛緩することで発生したと見なせる。

図 2.3.13: 初めの母音/ʉ/と最後の母音/ɯ/における、異なる断面での舌輪郭の比較。点線（黒）: 初
期形状、破線（赤）: 発話時の目標形状、実線(青): 推定された形状。
本研究では、3 次元的に形状を一致させる逆解析手法の提案および検証を行った。手法の検証の結果、逆解析より得られた形状は目標形状と 3 次元的によく一致した。筋線維応力は高い精度で値を推定することができ、応力の分布も一致していた。また、実際の母音/u/の形状に対して逆解析を行うことで、異なる画像から得られる形状においても逆解析を行うことができることを示した。推定された形状と目標形状にはわずかに誤差が残ったが、筋線維応力は過去の研究で得られた筋活動と同様の傾向を示した。以上のことから、本研究で提案した手法は、3 次元的に形状を一致させると共に、発話時のような複数の舌筋線維応力を推定するのに十分妥当であると考える。これらをまとめた成果は学術雑誌に掲載済みである(Koike et al., 2017, J Biomech, 64, 69)。一方で、本アプローチ
においては、物理モデルが拘束条件となるため、推定精度の改善に向け、より個別の舌変形を模擬した舌変形モデルの構築が今後の課題である。
(3) ヒト運動機能の神経筋骨格系シミュレータの開発（東京大学

山村直人，高木周，大阪大学

野村

泰伸）
A) 緒言
本研究テーマでは、ヒトの運動制御メカニズム解明および神経疾患に起因する運動障害の発現メカニズムの解明を目指して、大規模計算ニューロメカニクス研究の基盤プラットフォームと76

なる神経・筋・骨格系の統合シミュレータの開発を行っている。昨年度までに、筋線維への神経入力信号のもと、筋線維の収縮単位であるサルコメアを構成するミクロな収縮タンパクの確率的な挙動から、マクロな 3 次元有限要素モデルの変形を再現するマルチスケール・マルチフィジックス骨格筋シミュレータを開発した[1]。
本年度は、静止立位時の重心動揺の再現を目指して、医用画像から再現されたヒト下腿の 3次元筋骨格モデル(図 2.3.14)に対して、収縮タンパクの確率的な挙動から骨格筋の収縮挙動および関節運動を再現するための筋骨格系シミュレータの開発を行った。そのために、①3 次元骨格筋モデル内の 1 次元筋線維集団の設定、②遅筋・速筋からなる筋線維モデルの導入、③骨格・関節運動モデルの導入、④骨格筋有限要素モデルの並列化について検討を行った。
本開発により、医用画像から再現された 3 次元筋骨格モデルにおいて、ミクロな収縮タンパ
クの確率的な挙動から、骨格筋の収縮および関節運動に至るマルチスケールシミュレーションが実施できるようになり、ヒトの運動制御メカニズム解明のための基盤となるシミュレータが完成した。（「筋骨格」は筋と骨、「骨格筋」は筋のみ、を意味する。）

図 2.3.14: 3D musculo-skeletal model of lower limb for stationary standingB) 3 次元骨格筋モデルの 1 次元筋線維集団の設定
骨格筋は細長い筋線維の集合体で、骨格筋筋線維は心筋とは異なり、骨との付着部位（起始から停止）を 1 本で結ぶ細長い細胞である。骨格筋筋線維は筋収縮において、活動電位の発生と伝播を担い、収縮タンパクによる収縮力発生の起点となる。本研究で開発した骨格筋収縮モデ
ルは、複数のサルコメアが直列に並んだ 1 次元筋線維モデルを 3 次元連続体モデル内に埋め込んだモデルである（図 2.3.15）[1]。

図 2.3.15: 1D muscle fiber model embedded in 3D continuum-mechanical finite element model.

77ここでは、医用画像から再現されたヒト下腿の 3 次元筋骨格モデルに対して、その起始と停止をつなぐ複数の 1 次元筋線維の走行を設定する。筋線維の走行は起始から停止まで滑らかな曲線でつながっているため、本研究では、起始と停止を入口、出口とするような筋内のポテンシャル場において、そのストリームラインを 1 次元筋線維の走行として設定した[2]。
図 2.3.16 は下腿三頭筋モデル（ヒラメ筋および腓腹筋）の起始と停止および筋内に設定された筋線維の走行である。ヒラメ筋は 3543 本、腓腹筋は 1414 本のそれぞれの起始から停止までをつなぐ 1 次元筋線維の走行である。これら筋線維を等分割（200 分割）して複数のサルコメ
アモデルを筋線維に設置し、筋収縮モデルの解析に用いた。

(a)Origin and insertion(b) 3D muscle fiber trajectories

図 2.3.16: 3D muscle fiber distributions in triceps surae muscleC) 遅筋および速筋の筋線維モデル
骨格筋筋線維は主に、その代謝特性などの違いから、持続的な収縮が可能な遅筋線維と瞬発的な収縮が可能な速筋線維に分けられる。筋によってその含有割合は異なり、抗重力筋であるヒラメ筋では遅筋線維が支配的（90%）であり、一方、静止立位においてバランスを取るのに重要な腓腹筋では、その割合は 50%程度である。遅筋線維および速筋線維は、それがつながる運動神経線維が異なるほか、活動電位発生のしきい値や活動電位の伝播速度も異なっている。そこで本研究では、次式の筋線維の活動電位モデルにおいて、新たに遅筋線維モデルと速筋線維モデルの導入を行った。
𝜎

𝐴

𝐶

𝐼

(2.3.7)𝑉

ここで、𝜎および𝑉𝑚 は筋線維の電気伝導度および膜電位である。𝐴𝑚 は筋線維の体積に対する膜表面積の割合、𝐶𝑚 は細胞膜の単位面積あたりのキャパシタンスである。イオン電流 𝐼ion は細胞膜を通過するイオン電流の総和であり、膜電位に依存する。筋線維の違いは、細胞膜のキャパシタンス、イオンチャネルの密度、イオン電流の動態の違いなどで与えられる[3]。

78図 2.3.17 に遅筋線維と速筋線維の活動電位の挙動の違いを示す。図 2.3.17 (a)は遅筋線維と速筋線維の活動電位が発生する入力電流の違いであり、速筋線維に活動電位を発生させるためには、遅筋線維と比べて大きな入力電流が必要であること、また、活動電位の立ち上がり時間の違いが再現されている。また、図 2.3.17(b)は活動電位の筋線維長手方向の伝播速度で、速筋線

5050Slow fiber Iapp= 50.0
Fast fiber Iapp= 50.0
Fast fiber Iapp=150.0

0

-50

-100
0

5
Time /ms

10(a) Generation of action potentials

Membrane potential [mV]Membrane potential /mV維の活動電位の長手方向伝播速度は遅筋線維の 2 倍程度となっている。

2.48 m/s

01.25 m/s-50
Fast fiber
Slow fiber

-100
0
1
2
3
4
5
Distance along the sarcolemma [cm](b) Propagation of action potentials

図 2.3.17: Differences of behavior of action potential in slow and fast muscle fibers
D) 骨格・関節運動モデル
骨格筋は骨と付着しており、骨格筋の収縮力は腱を通して骨に伝えられ、関節運動として身体動作が生じる。したがって、静止立位時の重心動揺の運動を再現するためには、骨格筋に加えて骨や関節などを考慮する必要がある。ここでは、骨を剛体と仮定し、関節を拘束条件として与えることで骨格系の運動を定式化した[4]。
剛体単体の運動方程式は次式で与えられる。
𝑚𝑣

𝑓

𝐽′𝜔

𝜔

(2.3.8)𝐽′𝜔

(2.3.9)𝑛

ここで、𝑚は剛体の質量、𝑱′は慣性テンソルである。𝒗 、𝝎 および𝝎 はそれぞれ、剛体重心の並進加速度、剛体重心まわりの角速度および角加速度である。 𝒇 および𝒏 は剛体重心に加えられる並進力およびモーメントである。剛体の並進方向の運動方程式について、位置、速度および加速度の関係は Newmark-beta 法を用いて記述した。また、回転の運動方程式では、剛体の姿勢𝑹をオイラーパラメータにより記述する手法を採用した[5]。オイラーパラメータは単位
ベクトル𝒏のまわりに角度𝜃回転させたベクトル𝜽であり、姿勢𝑹は次式で定義される。
𝑹𝒏⊗𝒏

cos 𝜃 𝑰𝒏⊗𝒏sin 𝜃 𝛀𝐧(2.3.10)オイラーパラメータ、角速度および角加速度の関係は Newmark-beta 法を用いて記述した。

79

関節の拘束条件は関節の接続点における相対速度が 0 であるとして設定し、Lagrange 未定乗数法により剛体の運動方程式に付加した。また、骨格筋モデルは超弾性体有限要素モデルで定式化されており、剛体と骨格筋有限要素節点の結合は、拘束条件をその接続点における相対速度が 0 であるとして、Lagrange 未定乗数法により剛体の運動方程式に付加した。最終的に、筋骨格系の運動シミュレーションでは、骨格筋の運動方程式、剛体の運動方程式、関節に関する
拘束条件および骨格筋と剛体の接続に関する拘束条件を連成して解くことになる。
図 2.3.18 は下腿三頭筋モデル内と剛体モデルの連成による静止立位のシミュレーション例である。図 2.3.18 (a)は足関節角度 90 度を初期状態とし、その後、重力のみを付加したときの計算結果である。足関節 90 度では、アキレス腱や下腿三頭筋は伸ばされた状態にあるため、下腿三頭筋が弾性回復することにより骨格は後方に傾き、その後、重力により転倒する。また、図2.3.18(b)は初期状態からまず重心を前方に 20 mm 移動して骨格を前方に傾け、その後、重力を付加したときの計算結果である。この場合、下腿三頭筋の弾性回復による後方への移動より、重力による骨格の前方への運動が勝り、そのまま前方へ傾倒した。図 2.3.17（c）は (b)と同様に重心を前方に 20 mm 移動した後、重力とともに筋収縮を与えた計算結果である。
（b）とは異なり、収縮力によって、骨格は後方に引き戻されている。
このように、骨格筋の収縮から関節運動に至るマルチスケールシミュレーションが実施できるようになり、静止立位姿勢において、ヒトの運動制御メカニズム解明のための基盤となるシミュレータが完成した。

図 2.3.18: Coupling simulation of rigid body bone model and finite element skeletal muscle
model for human posture control. Simulation results (a) of backward fall, (b) of forward falland (c) taken into account muscle contraction.

80E) 筋骨格系シミュレータの並列化
筋骨格系シミュレータは、階層的に、複数の骨格剛体モデル、複数の骨格筋有限要素モデル、複数の 1 次元筋線維モデルおよび複数のサルコメアモデルからなるため、静止立位の重心動揺のような秒スケールの問題を解くためには、高度な並列処理が必要である。
複数の 1 次元筋線維モデルの並列化はすでに対応しており、今回、下腿の静止立位モデルの計算では、最大 1600 ノードの並列計算を実施した。また、新たに追加した骨格系剛体モデルでは、全身の骨格を考慮してもせいぜい数百のモデルなため、計算は 1 ノードで十分である。ここでは、新たに複数の骨格筋有限要素モデルの並列化について検討を行った。
骨格筋モデルの並列化では、有限要素モデルを複数の領域に分割する手法を用いた。有限要素モデルの領域分割にはグラフ理論を用いたライブラリ METIS [6]を利用し、連立一次方程式の解法には前処理付きの CG 法を用いた[7]。図 2.3.19 に 20 万自由度のモデルにおいて、分割数を 8 から 512 まで増加させたときの計算速度の向上を示す。512 分割の並列計算モデルは 8 分
割モデルに比べて約 20 倍の速度向上となった。

100
Speed upfrom 8 cores sim.

Flat MPI

10

1

10

100
Number of procs

1000

図 2.3.19: Effect of parallelization of finite element modelF) 結言
本研究では、ヒトの運動制御メカニズムおよび神経疾患に起因する運動障害の発現メカニズムの解明を目指して、力学的および生理学的に精緻な詳細筋骨格系マクロ−ミクロ階層統合モデルの開発を行った。今回は、静止立位時の重心動揺の再現するために、医用画像から再現された
ヒト下腿の 3 次元筋骨格モデルに対して 1 次元筋線維モデルを設定し、新たに骨格・関節モデ
ルの導入することで、ミクロな収縮タンパクの確率的な挙動から、骨格筋の収縮および関節運動に至るマルチスケールシミュレーションが実施できるようになり、ヒトの運動制御メカニズ
ム解明のための基盤となるシミュレータが完成した。
今後は、さらに神経モデルとのカップリングを行うことで、身体運動の神経制御メカニズムの解明を目指していく予定である。

81

参考文献
[1] Yamamura et al., IEEE 16th International Conference on Bioinfomatics and Bioengineering,pp. 338-341, 2016.
[2] Choi et al., PLoS ONE, 8(10), e77576, 2013.
[3] Shorten et al., J. Muscle Res. Cell Motil., 28, pp. 293–313, 2007.
[4] Hirasawa, Master Thesis, The University of Tokyo, 2014.
[5] Rubin, Comput. Math. Appl., 53, pp. 219-231, 2007.
[6] http://www-users.cs.umn.edu/~karypis/metis/[7] Saad, Iterative methods for sparse linear systems, 2nd edition, SIAM, Philadelphia,PA, 2003.

82

3.

心臓シミュレーションと分子シミュレーションの融合による基礎医学と臨床医学の架橋（サブ課題C）
以下では目標１、目標２に大別して実施内容を説明する。

(1) 目標１（心不全の解明と治療を目指すマルチスケール心臓シミュレーション）
心臓シミュレータ UT-Heart と分子シミュレータ CafeMol（京都大学・高田教授、理研・金田研究員）を融合させることにより、ミクロ・マクロ間の相互作用により病態が進行する心不全の解明と最適治療を可能とするマルチスケール心臓シミュレーションを実現することを目指す。これまでの
UT-Heart では、サルコメアを構成するミオシン分子のダイナミクスは、アーム部がばねで表された
1 自由度のミオシンヘッドが、ヌクレオチド状態に依存して確率的に首を振るモンテカルロシミュレーションにより表現されている。このサルコメアモデルにより生理学的に妥当な心拍動、血行動態が再現されただけでなく、心拍動におけるミオシン分子の協調性の役割などにつき重要な知見が得られた[1]。しかしミオシン分子を構成するアミノ酸残基の変異に起因すると言われる肥大型心筋症のような心疾患のメカニズム解明に切り込むためには、より実際の分子構造に立ち入ったシミュレ
ーションモデルの導入が望まれる。そこで図 3.1 に示されるように、ミオシンやアクチン分子を構成するアミノ酸一つ一つの自由度とタンパクの立体構造を考慮に入れた Cα-粗視化分子モデルに基づ
く CafeMol を用いた分子動力学シミュレーションと UT-Heart による心臓シミュレーションを連成させることにより、従来と次元を異にする心疾患の解明に挑戦する。

図 3.1: 「京」で用いて来たサルコメアモデルとポスト「京」で開発を目指す新たな分子モデル
CafeMol による分子動力学シミュレーションでは、図 3.2 にその概念を示すように、原子間相互作用のエネルギーに基づきアミノ酸粒子レベルでのポテンシャルエネルギーV が定められている。
このため Atomic Interaction based Coarse Grained (AICG) model と呼ばれる。しかし構造ベースのシミュレーションである為、参照構造となるヒトのアクチン、ミオシンの立体構造情報（V の定83

義式右辺各項における r0i ,  0i , 0i , r0ij ）が必要となる。

図 3.2: CafeMol の Atomic Interaction based Coarse Grained (AICG) model の概念及びポテンシャ
ルエネルギーと各項の係数
今年度の研究の位置付けを明確にするため、以下に昨年度（平成 28 年）の成果を要約する。
現時点における PDB (Protein Data Bank) 等のデータベースにおいては、特にヒト・ミオシンタ
ンパクの各ヌクレオチド状態 (ATP, ADP*pi, ADP,Rigor) における立体構造データが不足している。
従って、ここではヒトと配列相同性の高い生物種のミオシン構造を鋳型にしたホモロジーモデリン
グによりヒトのアクトミオシンモデルを構築した。図 3.3 はアミノ酸シーケンスの一例を示す。また
図 3.4、表 3.1 に結果を纏める。

図 3.3: Rigor 状態におけるミオシン S1 重鎖に関するホモロジーモデリングの一例

84図 3.4: Modeller から得られた 300 サンプル中ベストのミオシン構造（Rigor 状態）
表 3.1: ホモロジーモデリングの結果纏め（上:アクチン、下:ミオシン）

以上のように構築したヌクレオチド 4 状態モデル構造を用い、ATP 加水分解反応を模したミオシンの首振り運動が CafeMol によるポテンシャルの切り替えによって実現可能かフィージビリティス
タディを行った。その結果、図 3.5 に示すように粗視化 MD でも首振り運動の再現は可能であることを確認した。さらに図 3.6 に示すように、首振りの際にミオシンが発生する力のオーダーも in-vitroの 1 分子計測で測定される力の大きさ[2]と整合性がある事が分かった。

85

図 3.5: ヌクレオチド状態の変化に伴うミオシン・パワーストロークのシミュレーション

図 3.6: 分子シミュレーションから計測された各方向の発生力と実験結果[2]以上の成果を踏まえ、以下のように本年度（平成 29 年度）の研究を実施した。
分子シミュレーションと心臓シミュレーションの連成法については、自由エネルギー地形から各ヌクレオチド状態を離散化した確率モデルに置き換えることなく、直接的な連成を試みることとした。その場合、計算時間の観点から分子モデルの軽量化が望まれる。そこで首振りメカニズムのエッセンスを損なうことなく自由度が削減できるかどうか検討する為に、ミオシン構造変化の機構分析を行い、キーとなる残基群を調査した。

86最も重要なパワーストローク前後の Pi-release 状態および rigor 状態間でアミノ酸残基間のボンド
角(bond angle)および二面角(dihedral angle)の変化を調べた結果を図 3.7 に示す。パワーストロー
クに伴う変化は loop 構造をとる領域では変化が大きいが、βシートやαヘリックスの領域では小さく、これらを剛体として扱うことの可能性が示唆された。

図 3.7: パワーストロークにおける各残基のボンド角（上段）、二面角（下段）の変化（縦縞の赤：
Loop、青：Sheet、黄：Helix、青緑：その他）
さらに Relay helix と呼ばれる領域では全体としては構造変化が少ない中で、図 3.8 に示すように
C-terminal に近い狭い部位に大きな変化が見られたが、これは原子レベルの粒度でのシミュレーシ
ョン研究からパワーストロークの駆動力となる歪が蓄えられると考えられている部位に一致する。

87

図 3.8: Relay helix に認められた局所的な歪
A: CafeMol による評価結果 B: X-線回折（Koppole, 2007 [3]）から示唆される分子内歪集中
また同様に重要な歪の場所として指摘されている SH1、SH2 ヘリックスの間および SH1 ヘリッ
クスの内部にも、図 3.9 に示すように局所的な変化が観察された。

図 3.9: ATP 加水分解を力学エネルギーに変換する transducer 部分（power stroke に伴う SH1-SH2間の角度変化と SH1 の変形が認められる）
そこで Pi-release（プレパワーストローク）状態のモデルを初期構造にして、SH1 へリックス、
SH2 へリックス、Relay へリックスとそれらに関わる Go コンタクトにのみ rigor 型の歪みを与えて、CafeMol のエネルギー最小化のプロセスを実施した時、どの様な構造変化が起るか調べた。その結果、図 3.10 に示すように相当 rigor 状態に近い角度までレバーアームが回転、首振りする場合があることが分かった。但し、この結果は初期構造の Go コンタクトに依存するため、今後調査を継続する必要がある。

88

図 3.10: SH1、SH2、Relay へリックスおよびその周辺の歪エネルギーによるレバーアームの回転
そこでさらに CafeMol における bond angle と dihedral angle が梁の曲げと捩じりに対応することを考慮し、剛性比を 20：1 とした梁有限要素でミオシン分子モデルの Cα 原子主鎖を置き換え構造
力学的な観点からレバーアーム回転のメカニズムを考察した（図 3.11）。但し、文献調査から主鎖以外に重要と考えられる結合（Converter と Relay helix の間、SW2 と Wedge loop の間、βstrand
間、strand4 と P loop の間等）は図 3.12 の黒線で示されるように弱い剛性のトラス要素で結合した。
ATP を繋ぐ残基は SW1、SW2、P loop に存在するので、これら三つのループは加水分解反応時に化
学エネルギーが歪エネルギーに変換される際に重要な役割を演じると考えられる。そこでこれらの 3
ループに rigor 状態を基準とし Pi-release 状態へ至る過程で生じる相対変位を、結晶構造を参照して与えると、図 3.13(a)上段のようにレバーアームが回転した。
またアクチン分子と結合する Upper50Kと Lower50K の残基間の距離が隔たり、ミオシン分子の口が開いた変形が得られた。これは従来のATP ポケットとミオシン分子の口の開閉はシーソーのように逆位相になっているという学説を裏付けるものであると言えるが、ミオシン分子の開口量が大き過ぎる傾向にある。そこで次に、SW1、
SW2、P loop には相対変位を与えず、SH1、SH2、Relay の三つのへリックスに rigor 状態を基準とし Pi-release 状態へ至る過程で生じる歪を、結晶構造を参照して与えると、同図下段のように変形が生じた。ミオシン分子の開口量が適当になり、またレバーアームも大きく回転している。

図 3.11: 構造力学的観点からレバーアーム回転を考察するための梁有限要素モデル

89

図 3.12: 主鎖以外にトラス要素を用いて考慮した結合例（黒線；左から Converter-Relay 間、
SW2-Wedge loop 間、βstrand 間及び Strand4-P loop 間）

図 3.13(a): 梁有限要素法を用いたレバーアーム回転の構造力学的分析

図 3.13(b): 梁有限要素法を用いたレバーアーム回転の構造力学的分析
90

レバーアームの回転の方向は一方向からの図では分かり難いので別の角度からのスナップショットを図 3.13(b)に示したが、回転の方向も異なっていることが分かる。勿論、梁有限要素による解析は分子モデルそのものの挙動ではないが、分子モデルではとらえ難い構造力学的観点からの分析を容易にすることが出来る。
以上の結果は、Cα-粗視化シミュレーションによってもミオシンの重要な構造変化を捉え得る可能性があることを示している。このような調査を継続して ATP 加水分解に伴う微小な構造変化が増幅されるメカニズムを明らかにし、そのキーとなる残基群を同定することで分子モデルの自由度削減
あるいはクラスタ化を図ることが望まれる。なお、その際結晶構造を安定化する Go ポテンシャルをどう位置付けるが重要になってくると考えられる。最終的には心疾患において報告されているアミノ酸の変異がミオシン分子のダイナミクスを通じて心臓機能にどのような影響を与えるかを調べるための適切なモデル化を目指す。
本年度は昨年度作成した分子モデルの改良も以下のように行った。
昨年度作成した Pi-release 状態の分子モデルでは、鋳型となる 4PFO で構造化されていなかった
レバーアーム部分をモデリングするために、Rigor 状態(2MYS)と Pi-release 状態の Converter 領域
を構造アライメントし、Rigor 状態のレバーアームで置換した。しかし、このモデルでは排除体積項を削除したり、人為的な bridge 相互作用等の項を付与しなければ、Pi-release 状態から rigor 状態
へのパワーストロークの成功率を有意に上げる事が出来なかった。そこで今年度は Llinas らの研究
(Llinas et al, Developmental Cell 33, 401–412,2015) 記 述 に 基 づ き 、 構 造 化 さ れ て い な いPi-release のレバーアームをモデリングした。即ち、Llinas らは、ADP*Pi 状態と Pi-release 状態
を Upper-50K と N 末で構造アライメントした際、SWI と Ploop は両状態間で殆ど構造変化が無く、また Converter は限られた変化しかしないことを述べていることから、本研究では図 3.14 に示すよ
うに Scallop/1QVI を鋳型に作成した ADP*Pi 状態のヒトミオシン（モデル構造）における N 末とUpper-50K ドメインを Pi-release 状態モデルにおける対応領域に構造アライメントし、ADP*Pi 状
態のレバーアーム構造を Pi-release 状態モデルのレバーアームとする方針を取った。図 3.15 に改良
前後の構造を比較して示す。

図 3.14: Pi-release 状態モデルにおけるレバーアーム部の改良法
91

図 3.15: Pi-release 状態モデルにおける改良前後のレバーアーム部の比較
本年度は試験的連成シミュレーションを行うため分子モデルの軽量化を図った。即ち、アクチン分
子を「カーテンレールモデル」へ置き換え、以下のようにミオシン分子だけによるパワーストロー
クからレバーアーム C 末に取り付けたばね(Light meromyosin 相当)を介してマクロ有限要素モデルとの連成を行うこととした。
レバーアームの C 末に取り付けるばねは曲げモーメントも伝達できるよう 7 本とし、Pi-relase 状態モデルにおける安定位置で無歪になるものと設定した。先ず ADP*Pi 前半フェーズにおいては、安定化のため各ばね定数は 0.002857 kcal/mol/(Å^2)と仮定したが、これは rigor 状態でのばね定数(0.000290 kcal/mol/(Å^2))の 10 倍に相当する。ミオシンヘッドに対しては、アクチン長軸（z 軸）に平行なカーテンレール（Lower50K:12 本、Upper50K:2 本）を付与したが、そのばね定数は 0.000001
kcal/mol/(Å^2)と微小であり、実質上自由拡散と変わらない。なお、カーテンレールの位置は
Pi-release 状態においてアクチンと結合(native-contact)するグループの重心位置とした。またこのフェーズでは rigor 状態からのリカバリーストローク時の姿勢を素早く回復させる目的でヘッドの 2残基に対して位置拘束のばねポテンシャルを印加している。そのばね定数は 0.020000 kcal/mol/(Å^2)で、残基 ID:768、 417 である。この位置拘束は ADP*Pi 後半フェーズでは外している。次にPi-release 前半フェーズにおいては、レバーアームの C 末に取り付けた 7 本のばね定数を 0.000290
kcal/mol/(Å^2)とした。
これは Pi-release 状態後半および rigor 状態におけるばね定数と同一である。
またカーテンレールのばね定数は 0.000100 kcal/mol/(Å^2)にしており、ADP*Pi 状態よりは強く、
Pi-release 状態後半および rigor 状態の 0.1 kcal/mol/(Å^2)よりも弱い。これは、このフェーズがア
クチン結合に関して弱→強の中間状態であることを模擬したものである。このフェーズの最終ステ
ップにおけるヘッドの重心位置（拘束グループ残基 ID の重心位置）を次の Pi-release 状態後半フェーズへログを通して引き渡す。なお拘束グループ全体の重心位置がある条件を満たせば z 軸方向にも拘束を印加する様に設定しているが、その条件は殆ど満たされる事はない。Pi-release 状態後半フェーズにおいては、カーテンレールのばね定数の大きさは 0. 100 kcal/mol/(Å^2)としたが、これはPi-release 状態前半フェーズよりも強く rigor 状態と同一である。またカーテンレール方向（z 軸方向）にも拘束を設けるが、その位置については、先の Pi-release 状態前半における拘束グループ全体の重心位置を基に、その全体の重心位置からの相対位置を割り出し、その位置で拘束している。

92

以上の設定はこのフェーズが既にアクチンと強結合状態になっているという認識に基づく。最後にrigor 状態においては、ヘッドに対しては z 軸方向にも拘束を設けたカーテンレールを用いる点でPi-release 状態後半と同様であるが、本数は 13 本（Lower50K:9 本、Upper50K:4 本）を付与している。それらの拘束グループの内訳及び x-y 座標はあくまで、rigor 状態においてアクチンと結合(native-contact)する残基群及びその座標から先と同様の手順でグループ分けし、重心位置の計算を行ったものである。Pi-release 状態および Rigor 状態でのカーテンレールを図 3.16 に示す。

図 3.16: Pi-release 状態および Rigor 状態でのカーテンレールモデル（黄色表示）
以上のモデルを用い、7 本のばねの付け根は固定して 3 サイクル程度のシミュレーションを行った結果を図 3.17、3.18 に示す。図 3.17 はレバーアーム重心およびヘッド重心の時間変化を示す。図 3.17
ではレバーアーム重心座標（赤線）が+10Åから-60～-85Å程度まで変動するが、ヘッド重心座標（青線）も-20～-50Å程度の範囲で変動するので、ヘッドから見た相対的なレバーアームの変位は 30~40Å程度である。

図 3.17: 首振運動シミュレーションにおけるレバーアームとヘッド重心の z 座標変動

93

図 3.18: 首振運動シミュレーションにおける各フェーズのスナップショット
次にマクロ有限要素との連成について説明する。
今回のカーテンレールを用いた簡易的な分子モデルの枠組みでは、レバーアーム C 末部に取り付けられた 7 本のばねの歪は、刻々の C 末近傍残基の座標とばねの（もう一端の）付け根の座標の差から定まり、その歪によるばね張力がマクロの有限要素に伝えられる。一方、マクロ有限要素モデ
ルの運動解析結果から各要素の歪が求まり、それが 7 本のばねの付け根の変位としてフィードバッ
クされる具体的な手法は後述する。そのようなミクロ・マクロ間のフィードバックをどのようなタイミングで行うかが一つのポイントとなる。いまマクロ有限要素解析で 1 心周期（1s）を 1000 に刻んだステップ幅ΔT=1ms で計算を行うとすると、仮にミオシンの首振りが 1 心周期に 10 回生じても、ミオシンが何れかのヌクレオチド状態に滞在する間にマクロの解析が多数回行われることになる。換言すれば、マクロの一時間増分中にミオシンのヌクレオチド状態が変化する頻度は高くはない。では各ヌクレオチド状態に滞在するミオシン分子レバーアームに取り付けたばねの付け根の変位更新をマクロ解析の時間刻み毎に行えばよいかについて考えてみると、あくまでオーダーエスティメーションとして、半サルコメア長に相当する 1μm を 1000 ステップで割ると 1nm=10Åであるから、この程度の変位を一度に更新することは分子モデルの側から無理がある。よって分子モデル側ではマクロ解析の時間ステップΔT をさらに細かく刻んだステップΔt で計算を行う必要がある。
勿論 CafeMol においては本来の時間刻みは 10fs 程度であるが、Δt をそこまで細かくすることは計
算時間上不可能であるだけでなく、仮に 10fs 刻みで CafeMol を長時間実行したとしても粗視化分子シミュレーションの枠組みでは所詮化学状態の変化をシミュレートすることが出来ないことも事実である。従ってマクロ有限要素モデルからフィードバックされるばね付け根の変位に起因する分子構造の力学的変化（ATP ポケット寸法の変化など）から化学状態の遷移を確率的に定める経験式を別途定める必要がある。そこで次のような方法が考えられる。即ち、Δt に相当するばね付け根の変位を与えながら CafeMol を実行し分子構造の力学的変化とそのゆらぎをサンプルし、経験式から確率的に（モンテカルロ法で）化学状態遷移のスイッチを入れるかどうかを決定する。スイッチが ONになれば CafeMol のポテンシャルを切り替える。ON にならなければ次のΔt に移行する。Δt は以上の計算が可能な程度の大きさとする。なおこの場合 CafeMol の時間ステップは物理的な意味は持たない。本連成解析の概念を図 3.19 に示す。

94

図 3.19: 分子モデルと有限要素モデルの連成解析の方法
今年度は試験的連成解析を行うことを目標とし図 3.20 に示すような 8 要素からなる有限要素リン
グモデルを用いた。
有限要素法の定式化には以下のように連続体の歪エネルギー関数から応力が定まる異方性超弾性体を受動的構成側として用いた。
𝑊

𝑐 𝐼

3

𝑐

exp 𝑄
2

1

ここで、 𝐼 は左Cauchy-Green変形テンソル𝑪
𝐼

det 𝑪

Tr 𝑪

𝑄 はGreen-Lagrange歪テンソル𝑬
𝑄

𝑭 𝑭から定められる第一不変量

𝑏 𝐸

𝑏 𝐸

𝑏 𝐸

𝑪

2𝑏 𝐸

𝑰 ⁄2から定められる以下の２次形式である。
2𝑏 𝐸

2𝑏 𝐸

ここでそれぞれの成分は以下のように線維𝒇、シート𝒔、シートノーマル𝒏方向での成分表示に基づくものである。
𝐸𝐸

𝑬: 𝒇 ⊗ 𝒇, 𝐸
𝑬: 𝒇 ⊗ 𝒔, 𝐸

𝑬: 𝒔 ⊗ 𝒔, 𝐸
𝑬: 𝒇 ⊗ 𝒏, 𝐸

𝑬: 𝒏 ⊗ 𝒏
𝑬: 𝒔 ⊗ 𝒏

通常の心臓シミュレーションと同様のパラメータ(𝑐
3, 𝑏

10, 𝑏

2, 𝑏

71.8Pa, 𝑐

600Pa, 𝑏

5, 𝑏

6, 𝑏

2)を採用した。リングモデルでは𝒇に周方向、𝒔に軸方向、𝒏に半径方向を割り当てた。
またミクロ・マクロ間の具体的な連成は以下のように行った。

95

半サルコメアモデルのミオシンフィラメントのすべりΔ𝑧と連続体の線維方向ストレッチΔ𝜆の間には以下の関係が成り立つ。

1
𝑆𝐿 Δ𝜆
2

Δ𝑧

そしてこのミオシンフィラメントのすべりはレバーアーム C 末に接続されたばねの歪ξを Δ𝑧だけ減少させる。このことを考慮して以下のようにマクロの時間ステップ∆𝑇ごとのアクティブ張力を以下のように決める。
各ミクロの時間ステップ∆𝑡において、連続体としてのアクティブ張力はカーテンレールに結合しているミオシン分子レバーアーム C 末に接続されたばねにかかる力の総和から以下のように定義され
る。

𝑅𝑆
𝑆𝐴

𝐹act

𝑑𝑊rod
𝜉
𝑑𝜉

𝛿,

ここで𝑅𝑆は心筋中の筋原線維の占める割合、𝑆𝐴 はアクチンフィラメント 1 本あたりの断面積である。
ミクロ-マクロカップリング計算においては、マクロは∆𝑇ごとに力学的平衡をとる。その際、次のように計算を進める。まずは、 𝑇, 𝑇

∆𝑇 間の仮のアクティブ張力を以下のようにミクロの時間ステップ∆𝑡ごとの力の平均値とする(∆𝑇

𝑛∆𝑡)。

,

𝐹act

𝑅𝑆
𝑆𝐴 ∙ 𝑛

∆

𝛿 ,

𝑑𝑊rod
𝑑𝜉

𝜉

この際、対応する剛性も計算しておく。

,

𝐾

𝑅𝑆
𝑆𝐴 ∙ 𝑛

∆

𝛿 , 𝑘

𝑘 ,

𝑑 𝑊rod𝑑𝜉𝜉

ここで𝑘 , は k 番目のステップの以前で対象とするミオシンがカーテンレールに結合した直近のタ
イムステップ数である。ただしそれが時刻𝑇以前である場合にはゼロとする。上記の剛性を用いて連続体モデルにおいてはアクティブ張力を以下のように時刻𝑇∆𝑇での線維方向ストレッチ速度で修

正する。
𝐹act

ここで

,

𝐹act

∆𝑇
𝑆𝐿 , 𝐾4∆𝜆𝜆𝜆の項は , 𝐹act から時刻 T でのストレッチ速度の影響を取り除くためのもので、はそれを時刻𝑇∆𝜆の項∆𝑇でのストレッチ速度で入れ換えることを意味する。これにより時間幅 𝑇, 𝑇

∆𝑇

96

での心筋の剛性を適切に取り入れた安定なマクロの陰解法が実現された。Newton 法での反復では仮のアクティブ張力と剛性を再利用するのでミクロのタイムステップを繰り返し実行する必要はない。

図 3.20: リングモデルによる試験的連成解析の説明図
ここではレバーアーム C 末に取り付けた（7 本合計の）ばねの特性を図 3.21 のように 2 通り仮定し、さらに歪が負になるとミオシンが解離するモデルについてもシミュレーションを行った（合計 3 ケース）。図 3.22～24 にシミュレーション結果を示す。図 3.22 は t=0、150ms、300ms における変形のスナップショットを 3 ケースについて比較したものである。着色は周方向（線維方向）歪を表す
が、線形ばねモデル、非線形ばねモデル、非線形ばね＋歪依存解離モデルの順に変形とひずみが増大していることが分かる。非線形ばね＋歪依存解離モデルではリング内側（心内膜側に対応）での歪は 15%程度となっている。これを定量的に示したものが図 3.23 である。さらに図 3.24 には各ミオシン分子の発生力の変化が示されている。非線形ばね＋歪依存解離モデルでは収縮が早く起こる
ため他の 2 ケースと較べ発生力が小さくなっているが、負の力は殆ど生じていないことが分かる。

force [pN]

20
10

線形バネ
非線形バネ

0
-10
-20
-10

0

10

strain [nm]図 3.21: レバーアーム C 末に取り付けたばね（7 本のばねの合計）の特性

97

図 3.22: リングモデル連成解析結果の比較。着色は周方向歪を表す。

線維方向ストレッチ

1

0.95

0.9

0

線形バネ
非線形バネ
非線形バネ＆ひずみ依存解離

100

200
Time [ms]図 3.23: 周方向ストレッチの時間変化の比較

98

300

図 3.24: 有限要素に埋め込まれたミオシン分子の発生力の比較
一方、現在の UT-Heart を用い、ミオシンの状態遷移と心臓の拍動性能に関する基礎的な検討を行った。図 3.25 は UT-Heart で用いられている 2 段階パワーストロークモデルの各状態を表す模式
図、並びに第 1 段階首振り後の状態 XBpostR1 からの（ミオシンレバーアームの歪の関数としての）各遷移率曲線を示す。緑は首振り前の状態 XBpreR へ戻るリバーサルストロークの遷移率曲線を表すが、破線（No Trap Model）と実線（Trap Model）の 2 種類のパターンを仮定し心臓拍動へ及ぼす効果を比較した。実線ではレバーアームの歪（ばねの伸び）が 10nm 程度を超えると ATP ポケットが閉まることによってパワーストローク前の状態 XBpreR に戻れなくなる（trap される）との仮説に基づき遷移率をゼロとしている。図 3.26 は、No Trap Model と Trap Model に基づく心拍動の結果を比較したものである。Trap Model では高い心内圧が維持される時間、即ち拍出期が長くなっており、その結果 EF(Ejection Fraction)が顕著に増加していることが分かる。これは心筋内で弛緩のタイミングが早かった箇所が周囲から引き伸ばされて trap 領域に入り、XBpostR1 の状態を維持することで弛緩タイミングのずれを補正し、周囲と同期化することが原因である。本 Trap Model は心筋片の急な引っ張りに対する張力応答(Stretch activation)、サルコメア自励振動のカルシウムイオン濃度および温度依存性など幾つかの実験事実も上手く説明できる。勿論、今後十分な検証が必要であるが、分子シミュレーションモデル作成においては、図 3.27 に示されるようにレバーアームから伝達される張力が ATP ポケットの開閉に及ぼす影響や ADP および Pi の再結合の過程などを取り込む
必要性のあることが示唆された。

99

図 3.25: 現在の UT-Heart の 2 段階パワーストロークモデル、および第 1 段階首振り後の状態 XBpostR1からのレバーアーム歪に依存した各遷移率

図 3.26: No Trap Model と Trap Model に基づく心拍動の比較

図 3.27: 分子シミュレーションモデル作成への示唆

100本年度は以下のように心不全治療と外科手術の慢性期予測手法の研究に着手した。
心臓は負荷に対して大きくその形態を変化させて適応する。大動脈弁狭窄症、高血圧など収縮期
心室圧が過度に上昇するような条件（圧負荷）では壁厚の増加（求心性肥大）が、弁閉鎖不全症など拡張期の容積が過度に増加するような条件（容量負荷）では壁が伸展されたように変化し壁厚は減少（遠心性肥大）する。これらの変化は細胞の形態変化に基づいており、圧負荷に対しては細胞の横径が、容量負荷に対しては縦径が増加するが、こうした変化はそれぞれの方向へのサルコメア数の増加を伴う。病的負荷条件およびその是正に伴う慢性期効果予測のために、脚ブロックを伴う
心不全に対する心臓再同期療法、および先天性心疾患に対する手術療法を受けた症例の検討を自治
医科大学の協力の下に開始した。
脚ブロックのモデルとして左心室自由壁からのペーシングを行った動物実験では最も収縮が遅れる心室中隔壁が肥厚するのに対し最初に収縮する自由壁側は壁の菲薄化が起こることが報告されている（Circulation 1998）。我々がこれまでに行った左脚ブロック症例に対するモデル化研究（Okada
et al. J Mol Cell Cardiol 2016）でも同様の症例が認められたが（図 3.28 左）、一方全周性に菲薄化が認められる症例も存在した（図 3.28 右）。心室壁各セグメント（中段に色で示す）における周方向の歪の時間経過を比較したところ（図 3.28 下段）、全周性の菲薄化が認められた症例ではほとん
どのセグメントで伸展が生じており、そのタイミングは心筋の弛緩期（早期に興奮する中隔側では後半、遅れて興奮する自由壁では前半）であった。一方自由壁で肥厚が見られた症例では自由壁の収縮は心室圧が上昇してから始まる傾向が認められた。

図 3.28: 左脚ブロック症例における壁厚（上段）と色で示したセグメント(中段)における周方向の歪の心周期における時間変化（下段）。左の症例では最も収縮が遅れる青の領域で壁厚が厚く、最初に収縮する白の領域で薄くなっている。一方、右の例では全周性に収縮と過伸展が起きており、壁厚は全周性に薄くなっている。最初に興奮が起きる領域をペーシングとして記した。
この結果は圧負荷、および容量負荷モデルとして捉えることができ、また心臓同期療法による各セグメントの収縮タイミングの同期が負荷条件の是正をもたらすことが期待されるが、後ろ向き研究101であったために、十分な検討が行えるだけのデータが得られなかった。今回は形態変化について十分な検討が行えるように、治療前と治療後 6 か月の時点で心臓 3D エコーを記録する計画を承認していただき 2 例について 6 か月後の記録が終了した。これに加えて 3 例で経過観察が継続中である。
先天性心疾患はさらに大きな負荷が心臓に加わり、手術による構造の変化が大きな負荷の軽減をもたらす。先天性心疾患においても長期に渡る経過中の十分なデータが得られる症例を 5 例選択し、順次データの整理を実施している。1 例については解析に着手した。
肥大につながる分子メカニズムについては活発な研究が行われているが、ヒトの病気で見られる肥大のメカニズム、特に圧負荷と容量負荷による信号をどのように峻別しているかについては解明されていない。様々な経路が提唱されているが、その中でも近年注目されているサルコメアの構成
タンパク、タイチン(titin)を介した経路の検討には、分子シミュレーションと心臓シミュレーションの連成によってミクロとマクロの力学的相互作用を明らかにすることが重要である。

図 3.29: 半サルコメアの構造とタイチン内のメカノセンサー部位の分布。Z-帯、I-帯、M-帯の各領域に心肥大のシグナル伝達に関与するメカノセンサーが存在する。
タイチンはサルコメアの主要な構成要素であり、当初は構造の安定性に寄与するものと考えられていたが、近年、メカノセンサーとしての役割が注目されるようになった（図 3.29）。
図 3.30 にサルコメアの模式図を示す。タイチンは Z 帯から M 帯までサルコメアの半分にわたって存在する巨大分子で A 帯からはミオシンフィラメントに結合している。Z 帯、I 帯、M 帯領域のそれぞれに伸展に反応するメカノセンサーと考えられている構造がある。次図に示すように I 帯領域のセンサーは容量負荷で伸展されるが圧負荷では伸展されにくいと考えられる。
一方、収縮期にはミオシンフィラメントとアクチンフィラメントの間に架橋（クロスブリッジ）が形成され、かつサルコメア長は比較的短いために、I 帯領域は伸展を免れるが、力は下に示す黄色で示した経路で伝わり、ミオシンフィラメントに歪が生ずるものと予想される。タイチンはこの部

102

分ではミオシンフィラメントに結合しているため M 帯領域のセンサーが伸展されるものと考える。
このような可能性について Arizona 大学の Granzier 教授と意見交換を行った。
今後はマクロレベルの負荷の変化に伴う、クロスブリッジ動態およびサルコメアの挙動との関係について検討を進め、治療後の長期予後予測へとつなげる予定である。

図 3.30: 圧負荷および容量負荷におけるサルコメア長（右上）とタイチン分子の変化（左上）。クロ
スブリッジが形成される時相に加わる圧負荷は I-帯のセンサーより M-帯のセンサーに歪を加える（下）。
(2) 目標２（創薬を加速する心毒性スクリーニングシステムの開発）
本テーマでは、重点課題 1（奥野恭史・重点課題責任者）の寺田 透博士との連携により、心臓シミュレータ UT-Heart と心筋細胞イオンチャネルの分子シミュレーションを融合し、創薬における候
補物質の分子構造と遺伝子情報のみから不整脈のリスク評価を可能とするシステムを開発することを目的とする。具体的には、既に開発済の図 3.31 に示される UT-Heart と細胞パッチクランプ実験と組み合わせた「ハイブリッド心毒性スクリーニングシステム」の細胞薬理実験（赤破線内）の部
分を分子シミュレーションに置き換える。本ハイブリッド心毒性スクリーニングシステムは、図 3.32に示すように、既に 12 種類の薬剤に対し検証を終え、従来のガイドラインにおける指標（例えばQT 延長）に較べ偽陽性も出ないなど優れた判定が出来ることが示されている[4]。

103図 3.31: UT-Heart と細胞パッチクランプ実験と組み合わせた「ハイブリッド心毒性スクリーニングシステム」

図 3.32: ハイブリッド心毒性スクリーニングシステムによる 12 薬剤の解析結果（Okada et al,
Science Advances, 2015 [4]）図 3.33 はカリウムイオンチャネルの立体構造を横（細胞膜切断面）および上（細胞外）から見たものである。四つの青色領域は電位センサードメイン（VSD）の立体構造、それによって囲まれる
薄茶色部はポアドメインの立体構造である。電位センサードメインを構成する 4 本のヘリックスのうちの一つ S4 の正に荷電した残基が電位センサーの役割を果たし、脱分極状態では S4 が細胞外側に移動する（Up 状態となる）ことによりチャンネルが開く。また分極状態では S4 が細胞質側に移動する（Down 状態となる）ことによりチャネルが閉まる。

104

図 3.33: カリウムイオンチャネルの立体構造。左図は横（細胞膜切断面）からの概観、右図は上（細胞外）からの概観（重点課題 1 寺田博士作成）
分子シミュレーション開発は二つのパートに大別される。第 1 は電位センサードメインの立体構
造変化に関するシミュレーション手法の開発、第 2 は不整脈を誘発する薬剤と心筋イオンチャネルとの間の相互作用を予測するためのポアドメインでのドッキングシミュレーション手法の開発である。年次計画に従って第 2 のパートは次年度の課題とし、今年度は第 1 のパートについて以下のよ
うに研究開発を進めた。
本研究開始時点ではヒト心筋の電位依存性イオンチャネルについては立体構造が決定されていなかったため、ここではラットのカリウムチャネル Kv1.2 の S3 の C 末端側半分と S4 の N 末端側半分をラットの Kv2.1 に置換した（図 33 にピンク色で表示）キメラの立体構造（PDB ID: 2R9R）を用いた。そして図 3.34 に示すように VSD の立体構造を MD シミュレーションにより水和した脂質二重層中で平衡化した後、Jensen ら[6]の長時間シミュレーションのトラジェクトリをガイドとして、
targeted MD を実施し、Up 状態から Down 状態に至るトラジェクトリから 128 の構造を選択した。

図 3.34: VSD の分子シミュレーションに於ける初期構造の作成（重点課題 1 寺田博士作成）
次にこれら 128 の構造それぞれに対して分極状態の膜電位 Vm = −0.1 V もしくは脱分極状態の Vm= 0.0 V の電場を与えて 0.1 μs の MD を実施した後、クラスタ解析により立体構造を microstateに分類し、一定時間後に microstate i から microstate j に遷移する数をカウントした。ここで立体構造空間は多数の状態（microstate）から成り、microstate 内の構造変化は、microstate 間の状態遷移に比べて速く起こるとし、また microstate i から j への遷移確率 Tij は、それ以前にどの
microstate にいたかに依存しないマルコフ状態モデルを仮定する。このような Tij からなる遷移行列を用いたモンテカルロシミュレーションを 1 秒間実施した。そして S4 ヘリックスの R296 残基の Z

105座標から、F233 残基の Z 座標を引いた値が正であれば Up 状態、負であれば Down 状態であるとみなした。モンテカルロシミュレーションの結果、Up 状態から Down 状態あるいは Down 状態からUp 状態に至る立体構造遷移の時間スケールは 10～100 μs であることが分かった。これはチャネル全体の動作の時間スケール（10 ms 程度）に比べて十分に速いので、今後は瞬時に平衡に達するとみなすことにする。
以上のシミュレーションに基づき VSD4 量体のすべてが Up になった場合チャネルは開口状態に遷移するとみなして図 3.35 のようなチャネル開閉動作の確率モデルを考案した。ここで K(0)および

q は MD から求めたマルコフ状態モデルの遷移行列に基づくモンテカルロシミュレーションにより評価する。なお図中に本モデルから計算される矩形電位波に対する開口率の時間変化特性を併記するが、これは実験結果に基づく Zagotta ら[7]の動作モデルの結果と類似している。

図 3.35: VSD の分子シミュレーションに於ける初期構造の作成（重点課題 1 寺田博士作成）
図 3.31 に示したハイブリッド心毒性スクリーニングシステムの破線部を MD シミュレーションに基づく方法に置き換えた新たな連成手法の概念を図 3.36 に示す。このシステムによれば、新たな薬剤に対してチャネル阻害率を評価するドッキングシミュレーションは毎回行う必要があるものの、膜電位に依存するチャネルの構造変化は確率モデルに置き換えているため、分子シミュレーションを再度行う必要はない。多数の候補薬に対する心毒性評価を効率的に行うことを目指した現実的連成法であると言える。

106

図 3.36: イオンチャネル分子シミュレーションと UT-Heart の新しい連成方法
本年度は UT-Heart で用いている ORd モデルの緩徐活性型 IKs チャネル部を、以上に説明したMD に基づく確率過程モデルによって置き換え心臓の興奮伝播シミュレーションを行った。但し IKs
チャネルは KCNQ1 遺伝子と KCNE1(mink)遺伝子によりコードされるが、本分子モデルは Shaker
遺伝子によりコードされる IKur と呼ばれるチャネルであるため図 3.35 に示した動特性は IKs の動特性（図 3.37）とは大きく異なる。次年度は 2017 年に Wang らにより決定された電子顕微鏡によCurrent (Arbitrary Units)

る IKr チャネルの立体構造[5]を用いたシミュレーションに着手する予定である。

図 3.37: 緩徐活性型 IKs チャネルの矩形電位波に対する電流特性107図 3.38 は心外膜側における膜電位および IKs 電流について、従来の ORd モデルに基づく心臓シ
ミュレーション結果と今回の分子モデルに基づく心臓シミュレーション結果を比較したものである（後者を ORd model+MD IKs と表示）。ORd model+MD IKs は ORd モデルに較べ応答が速いため、総てのチャネルの作用の合計としての膜電位波形も前倒し傾向にある。また IKs 電流に関しては位相が早まるだけでなく振幅も大幅に増加している。図 3.39 では 12 誘導心電図を比較したが、同様の理由により ORd model+MD IKs では QT 間隔（図中の矢印）が短くなっていることが分かる。図3.40 で は 心 臓 の 興 奮 伝 播 な ら び に 体 表 面 電 位 の 推 移 を ス ナ ッ プ シ ョ ッ ト で 比 較 し た 。 ORdmodel+MD IKs では興奮が早く醒めていることが分かるが、これは T 波が前倒しになったことに対応する。
以上のようにチャネル分子モデルの特性を端的に反映した連成解析シミュレーションが実現できた。

図 3.38: 心外膜側における膜電位および IKs 電流に関する両モデルの比較

図 3.39: 12 誘導心電図に関する両モデルの比較

108

図 3.40: 心臓の興奮伝播と体表面電位のスナップショット比較
さらに、チャネル阻害率と不整脈発生の関係を現在の UT-Heart で調べ心毒性の全体像を明らかにした。具体的には、複数のイオン電流（遅延整流性カリウム電流 IKr、IKs、ナトリウム電流 INa、
INaL、カルシウム電流 ICaL）の重畳抑制がどのように不整脈発生に影響を及ぼすかについて主な組み合わせについて調べた。興味深いことに、例えば Verapamil は dose を増しても心電図の QT 間
隔は延長するものの不整脈は発生せず、一方、Cisapride は QT 間隔の延長と共に、常用量の 30 倍
で TdP（致死性不整脈に繋がる頻脈）が発生する。これは従来のガイドラインで危険とされる QT
間隔延長が不整脈リスクの指標として必ずしも有効でないことを示唆しており、一つの薬剤により複数のイオン電流が重畳抑制されたことに起因する結果であると考えられる。なお、QT 間隔の定義
を含め、これまで提案されている各種の心電図指標を図 3.41 に示す。では、どのような複数のイオ
ン電流の重畳抑制がどのように不整脈発生に影響を及ぼすのか、その全体像を明らかにしておくことは、今後分子シミュレーションでモデル化するイオンチャネルの選択に関わってくる。
図 3.42 に IKr と ICaL を抑制した結果の一例を示す。黒色領域は不整脈発生領域、他の色は心電図指標の一つである Tpeak-Tend の変化量（ms)を示す。これより、IKr の抑制は不整脈をもたらすが、加えて ICaL の抑制が 20%を超えると IKr の効果を緩和することが分かる。さらに、Tpeak-Tend
値は不整脈発生領域に向かって一旦減少し、その後増加に転じる傾向を示していること、並びに、不整脈発生直前での Tpeak-Tend の限界値は一定ではない、ことから不整脈発生を予測する有効な心電図指標とはなり得ないことが分かる。

109

図 3.41: これまでに提案されている各種の心電図指標

図 3.42: IKr と ICaL 電流の各種抑制率の組み合わせに対する不整脈発生の有無と Tpeak-Tend 変化
量（カラーバーの単位は ms）
本研究ではさらに全体像を明確にするため、樹木構造接近法（CART）により解析を行った。樹木
構造接近法とは、各サブ集団の特徴が最も均一になるように母集団を変数により分割していく手法である。同じ変数で何度も分割できるため母集団に潜む非単調作用や変数同士の相互作用も含めて樹木形式で可視化して理解する事が可能となる。先ず、種々のイオンチャネル阻害の組み合わせに対し予め不整脈発生の有無をシミュレーションで明らかにした後、イオンチャンネルの阻害率を変数として、不整脈を起こしやすいサブ集団と不整脈を起こしにくいサブ集団に分割を行っていく。
この分析から、イオンチャンネル A の阻害率が X%以上で、イオンチャンネル B が Y%未満の場合
は不整脈が起こるという条件式が得られる。具体的には図 3.43 の樹木構造を上から下に向かって見て行くと、最初の分岐では 5 種類のイオンチャネルのうちで IKr チャネル阻害率 75[%]を境界線としてシミュレーション結果全体を不整脈発生し易い集団(左側の集団)と不整脈発生し難い集団(右側の集団)に分ける事が出来ることを示している。次の分割では、IKr チャネル阻害率 65[%]以下(最右110

の集団)では不整脈は全く発生しない一方で、85[%]以上(最左の集団)では 98.8％に不整脈が観測されたことを示している。このように心毒性の評価においては IKr チャネル阻害が大きな役割を果たしている事が分かる。さらに IKr チャネル阻害率が 65[%]以上 85[%]以下の集団における分割から IKs
チャネル阻害が不整脈の発生を促進すること、また ICa チャネル阻害が不整脈の発生を抑制する効
果のあることが分かる。以上の調査により、重畳効果も含めて不整脈発生に対して主要な役割を果たすチャネルは IKr、IKs、ICaL の三つであり、これらを中心に分子シミュレーションモデルを開発するのが良いことが分かった。なお以上の計算結果はデータベースとして公開されている[8]。

図 3.43: 樹木構造接近法(CART)による五つのイオンチャネルの不整脈発生に対する影響度解析
参考文献[1] Washio T, Okada J, Takahashi A, Yoneda K, Kadooka Y, Sugiura S, Hisada T, Multiscale heart
simulation with cooperative stochastic cross-bridge dynamics and cellular structures, SIAM
Multiscale Modeling and Simulation, 11(4),965-999(2013)[2] Sugiura S, Kobayakawa N, Fujita H, Yamashita H, Momomura S, Chaen S, Omata M, Sugi H,
Comparison of unitary displacements and forces between two cardiac myosin isoforms by the opticaltrap technique: Molecular basis for cardiac adaptation, Circ Res,82:1029-1034(1998)[3] Koppole S1, Smith JC, Fischer S,The structural coupling between ATPase activation and recovery
stroke in the myosin II motor, Structure,15(7):825-37(2007)[4] Okada J, Yoshinaga T, Kurokawa J, Washio T, Furukawa T, Sawada K, Sugiura S, Hisada T,
Screening system for drug-induced arrhythmogenic risk combining a patch clamp and heart111

simulator , Science Advances, 1（4）, e1400142(2015)[5] Wang W, MacKinnon R, Cryo-EM Structure of the Open Human Ether-à-go-go-Related K+
Channel hERG, Cell,169,3, 422–430.e10(2017)[6] Jensen MØ, Jogini V, Borhani DW, Leffler AE, Dror RO, Shaw DE. Mechanism of voltage gatingin potassium channels. Science, 336, 229–233, 2012
[7] Zagotta WN, Hoshi T, Aldrich RW. Shaker, potassium channel gating. III: Evaluation of kinetic
models for activation. J. Gen. Physiol. 103, 321-362. 1994[8] Okada J, Yoshinaga T, Kurokawa J, Washio T, Furukawa T, Sawada K, Sugiura S, Hisada T,
Arrhythmic hazard map for a 3D whole-ventricles model under multiple ion channel block British
Journal of Pharmacology (2018）DOI: 10.1111/bph.14357

112

4.

プロジェクトの総合的推進

プロジェクト全体で密な連携を取りながら運営が円滑に進むよう、業務計画を策定・審議し、運営委
員会や研究会を主催し、また成果を積極的に公開した。即ち、前段階からの成果を進化させるべく目標と体制を最適化して計画を立案・策定し、その計画に沿ってプロジェクト運営のためにテレビ会議、電話会議による運営委員会を 10 回開催し、全体の進捗状況を確認しつつサブ課題参画各機関による連携を推進した。得られた成果は学会で発表するほか、重点課題 2 主催のシンポジウムやワークショップ各 1
回において紹介した。計画や成果は外部有識者による諮問委員会を 4 回開催し重要事項を審議に掛けた。
また、全重点課題を含むシンポジウム等においてはポスター掲示などで成果をアピールし、重点課題広報責任者ワーキンググループにおいては計算科学の発展に資する議論を深め、活動により実践した。さらにホームページ運営では、成果を随時掲載することで、コンテンツ拡充を推進した。またニュースレター（電子版）6 編を発刊し、ホームページに掲載した。
HPCI 計画推進委員会、HPCI コンソーシアム通常総会、今後の HPCI コンソーシアムの活動に関する意見交換会、システム・重点課題アプリケーション連携推進委員会等、ポスト「京」開発に関係のある会議に参加し現状把握のための情報収集を行った。ソフトウェア普及においては、重点課題 2 のター
ゲットアプリケーションである Genomon2 の講習会などを 2 回開催した。
(1) 運営委員会
プロジェクト運営を円滑かつ効率的に行うために、重点課題 2 責任者、サブ課題責任者などで構成した運営委員会を実施している。研究開発及びプロジェクトの総合的推進の活動進捗の確認や、計算資源の配分、課題についての調整等を、サブ課題参画各機関と協調して行った。
平成 29 年度は以下の日程で開催した。
第1回

4 月 24 日(月)東京大学医科学研究所にて開催

第2回

5 月 19 日(金)東京大学医科学研究所にて開催

第3回

6 月 16 日(金)東京大学医科学研究所にて開催

第4回

7 月 21 日(金)東京大学医科学研究所にて開催

第5回

8 月 18 日(金)東京大学医科学研究所にて開催

第6回

9 月 25 日(月)東京大学医科学研究所にて開催

第7回

10 月 20 日(金)東京大学医科学研究所にて開催

第8回

11 月 20 日(月)東京大学医科学研究所にて開催

第9回

12 月 19 日(火)国立科学博物館講堂にて開催

第 10 回2 月 23 日(金)東京大学医科学研究所にて開催(2) 諮問委員会
プロジェクトの推進に資するため外部有識者からなる諮問委員会を開催した。研究開発の成果やプロジェクト体制、今後の方向性などについて、諮問委員から幅広く意見を伺い平成 29 年度の研究開発活動に反映させた。

113A) 諮問委員（◎委員長）
礒田治夫（名古屋大学脳とこころの研究センター

梶谷文彦（川崎医療福祉大学

高井義美（神戸大学大学院医学研究科

◎樋口

知之

教授）

名誉教授)

（情報・システム研究機構

教授）
統計数理研究所 所長）

B) 開催状況2017 年 5 月 24 日(水)成果報告書内容の審議および助言

2017 年 7 月 28 日(金)中間評価ヒアリング前の助言依頼

2017 年 9 月 5・6 日(火・水)中間評価ヒアリング後の報告および今後の方針への助言依頼
2017 年 12 月 4 日(月)中間評価結果の説明および実施計画への意見交換(3) シンポジウムおよびワークショップの開催
他機関と連携したシンポジウムやワークショップ等を開催した。
平成 29 年度に行った活動概要は以下の通り（共催、協力、協賛、後援を含む）。
A) シンポジウム(a) 名称：その予防・医療、時代遅れです－ビッグデータ×
シミュレーション×ポスト「京」＝DNA＊♰∞－
開催日：2017 年 11 月 13 日(月)場所：秋葉原 UDX ギャラリー
参加者数：92 名
主催：ポスト「京」重点課題 2
協賛：「システム癌新次元」文部科学省科学研究費補助
金 新学術領域研究 (複合領域 4701) がんシステムの新次元俯瞰と攻略

図 4.1: シンポジウムの様子

内容：医療応用をめざし最先端の生命科学研究を行っている研究者とともに、計算科学と生命
科学の視点から、ポスト「京」コンピュータを用いた医療技術・予測医療の可能性を考えることを目的としてシンポジウムを開催した。電気通信大学脳科学ライフサポート研究センター正本和人教授を基調講演者として招聘し、社会的対策が急務とされる認知症の発症予防法の確立に大きく寄与する可能性がある脳血液循環と脳機能との関わりについて最新の研究講演を行った。また各サブ課題の研究開発の最新状況については、今年度は参加研究者 3 名が登壇した。
アンケート回答では、「非常に充実した内容」「分かりやすく大変面白かった」「ポスト「京」の必要性を理解することができ大変有意義であった」等、満足度の高い結果を得た。
(b) 名称：第 4 回 大型実験施設とスーパーコンピュータとの連携利用シンポジウム
開催日：2017 年 9 月 8 日(金)場所：東京・秋葉原 UDX

6 階カンファレンス

主催：高輝度光科学研究センター、総合科学研究機構、高度情報科学技術研究機構

114

協賛：理化学研究所、ポスト「京」重点課題 2～9、他 6 機関
協力：ポスト「京」重点課題１
後援：文部科学省、兵庫県、茨城県(c) 名称：第 4 回「京」を中核とする HPCI システム利用研究課題 成果報告会
開催日：2017 年 11 月 2 日(金)場所：東京・コクヨホール
主催：高度情報科学技術研究機構
協賛：HPCI コンソーシアム、スーパーコンピューティング技術産業応用協議会
協力：理化学研究所計算科学研究機構、国立情報学研究所、他 11 機関
後援：計算科学振興財団、ポスト「京」重点課題実施機関(d) 名称：2nd International Symposium on Research and Education of Computational Science
開催日：2017 年 12 月 18 日(月)～12 月 19 日(火)場所：Koshiba Hall, School of Science, The University of Tokyo
主催：The Computational Science Alliance of the University of Tokyo
協 賛 ： Joint Usage/Research Center for Interdisciplinary Large-scale Information
Infrastructures, HPCI Consortium, Post-K Computer Priority Issue 1-9, RIKEN Advanced
Institute for Computational Science(e) 名称：The 8th AICS International Symposium
Roadmap on computer science and computational science in the future AICS and HPC
communities
開催日：2018 年 2 月 7 日（水）～8 日(木)場所：RIKEN Advanced Institute for Computational Science
主催：RIKEN Advanced Institute for Computational Science
後 援 ： Priority issues (9 issues) on Post-K computer, High Performance Computing
Infrastructure Consortium, Research Organization for Information Science and Technology(f)

名称：見える化シンポジウム 2018 「シミュレーション可視化の未来～計算科学ミュージアム
実現に向けて～」開催日：2018 年 2 月 24 日(土)場所：東京・日本橋ライフサイエンスハブ
主催：理化学研究所計算科学研究機構
共催：ポスト「京」重点課題 1、7、9
協力：ポスト「京」重点課題実施機関、高度情報科学技術研究機構

115B) ワークショップ(a) 名称：ポスト「京」重点課題 2 ワークショップ
開催日：2017 年 12 月 4 日(月)場所：東京大学医科学研究所
参加者数：35 名
主催：ポスト「京」重点課題 2
研究開発の進捗確認および報告の場、またプロジェクト
参加機関内での実質的な相互連携の推進を目的として、
プロジェクト参加研究者を対象にワークショップを開催した。参加研究者をはじめ諮問委員からの質疑応答では

図 4.2: ワークショップの様子

活発な議論が交わされた。
(b) 名称：KOBE HPC Summer School 2017
開催日：2017 年 9 月 4 日(月)～9 月 8 日（金）場所：理化学研究所計算科学研究機構

6 階講堂

共催：神戸大学計算科学教育センター、理化学研究所計算科学研究機構、他 3 機関
後援：ポスト「京」重点課題 9 課題、高度情報科学技術研究機構(c) 名称：2017 年度 理化学研究所神戸キャンパス一般公開
開催日：2017 年 10 月 14 日(土)場所：理化学研究所計算科学研究機構

6 階講堂

一般来場者数：3,509 名
内容：理研計算科学研究機構 6F 講堂に設けられた重点
課題の研究紹介や工作体験などを行う「神戸スパコンシ
ミュレーション」にて、重点課題 2 の研究内容を紹介したポスター展示と人体の臓器パズルの組み立て体験を実施した。生命科学への関心を促すため臓器について解図 4.3: ブースの様子

説した資料も提供した。
重点課題 2 ブースへの来場者数は約 300 名を記録し大変好評であった。
(d) 名称：RIKEN AICS HPC Youth Workshop FY2017
開催日：2018 年 2 月 5 日(月)～2 月 6 日(火)場所：理化学研究所計算科学研究機構

6 階講堂

主催：理化学研究所計算科学研究機構
後援：ポスト「京」重点課題 9 課題、計算科学振興財団、HPCI コンソーシアム、高度情報科学
技術研究機構(e) 名称：KOBE HPC Spring School 2018
開催日：2018 年 3 月 14 日(水)～ 3 月 16 日(金)116

場所：神戸大学計算科学教育センター 208 セミナー室
共催：神戸大学計算科学教育センター、理化学研究所計算科学研究機構、東京大学情報基盤センター、兵庫県立大学大学院シミュレーション学研究科
後援：ポスト「京」重点課題 9 課題、高度情報科学技術研究機構、計算科学振興財団(4) ポスト「京」重点課題広報責任者ワーキンググループ
ポスト「京」の研究成果やシミュレーションの重要性・必要性を、重点課題全体として効果的に普及を行うこと、また広報に関する重点課題間の連携推進等を目的として、広報責任者ワーキンググループが開催された。理化学研究所計算科学研究機構・重点課題 9 課題・高度情報科学技術研究機構の広報責任者および文部科学省で構成され、重点課題全体で実践する成果の普及活動等、計算科学の発展に資する議論・活動を行った。
アウトリーチ・広報活動を効果的に行うために、前年度より引き続いてプロジェクト全体で共通
する統一的な広報ツールを制作した。個々の重点課題の活動においてもプロジェクトの全体像や他
重点課題の取り組みが俯瞰的かつ直感的に理解されることを目的とする。具体的には、統一ロゴやシミュレーションの意味を広く深く掘り下げたポスター”シミュレーション図”等を制作し様々な媒体で情報発信を行った。
組織戦略の一環として、ポスト「京」ユーザ、周辺研究者関係者や一般等を対象としたポスト「京」とポスト「京」で行われるシミュレーション科学のブランディング確立に向けての活動を本格的に開始した。効率的に進めるために、ポスト「京」ブランドコンセプトタスクフォースを設置し協議を重ねワーキンググループで意見を取り纏めた。
平成 29 年度に開催したワーキンググループは以下の通り。
第7回

7 月 4 日(火)理研東京連絡事務所会議室、理研 AICS（TV 会議）にて開催

第8回

10 月 3 日(火)理研東京連絡事務所会議室、理研 AICS（TV 会議）にて開催

第9回

12 月 14 日(木) 理研東京連絡事務所会議室、理研 AICS（TV 会議）にて開催(5) ホームページ（日本語・英語）
昨年度に引き続きサイトコンテンツの充実および更なる拡充を行った。今年度は専門知識のない閲覧者を対象にした情報発信として、専門用語ページや研究開発内容を比較的容易な言葉を用いた解説ページを新しく設けた。また、全重点課題の統一ロゴや他重点課題や研究機関との連携活動を掲載することにより、プロジェクト全体の認知度向上のための横断的取り組みを紹介した。

117

図 4.4: 日本語版
図 4.5: 英語版(6) ニュースレター、パンフレット発行（電子版）
一般および企業等を対象として、研究開発の進捗や成果などを分りやすく説明することで研究活動の理解を深めるとともに、社会における重要性や価値の理解促進を図るためニュースレターを発行した。また、研究体制の変更に伴い、研究活動の概要を紹介したパンフレットを改訂した。

118

図 4.7: Vol.5, 6（9, 11 月）

図 4.6: Vol.3, 4（5, 7 月）

図 4.8: Vol.7, 8（1, 3 月）

(7) セミナー、講習会等の開催
A) 講習会SiGN-BN による遺伝子ネットワーク解析実習では、理論的な仕組みや遺伝子ネットワークのサイズに応じた三つのネットワーク構造推定アルゴリズムについて講義および実習を行った。
Genomon2 Tutorial では、主にがんのゲノム・トランスクリプトーム解析に興味のある実験研
究者、医学研究者を対象に、がんの後天的変異・構造変異・融合遺伝子を検出する一般的原理、
ゲノム解析に必要な実践的なテクニックを紹介する講習会を行った。ハンズオンセッションでは、Genomon2 の体験実習を通じてゲノム解析の流れについて理解を深めることを目指した。
開催場所：東京大学医科学研究所
表 4.1: 平成 29 年度講習会の詳細
7 月 27 日

遺伝子ネットワーク解析実習

1:30-5:30pm
10 月 30 日
1:00-4:30pm

遺伝子ネットワーク推定ソフトウェア SiGN-BN の講義と実習

Genomon2 Tutorial

がんゲノムシーケンス解析の原理と Genomon2 の紹介
Genomon2 ハンズオンセッション119

４－３．活動（研究会・受賞・書籍等）
1.

研究会

発表した成果（発表題目、口頭・発表者氏名発表した場所

発表し国内・ポスター発表の別）（所属機関）（学会等名）た時期

外の別心筋アクトミオシンの粗視化分子

金

第３回数理物理・統計力学

2017年

モデル

（UT-Heart）

セミナー

12月

心臓有限要素解析の医療・創薬へ

岡 田の応用（招待講演）（UT-Heart）

連続体と分子モデルを繋いだ心拍鷲動シミュレーション（招待講演）（UT-Heart）

心臓有限要素解析の医療・創薬へ

岡 田の応用（招待講演）（UT-Heart）

拍動する心臓のマルチスケールシ

鷲ミュレーション（招待講演）（UT-Heart）

2.

田

亮
純 一

尾

巧
純 一

尾

巧

10月

第2回CDMSI(ポスト「京」重2017年

点課題（7）)研究会

7月
2017年

第31期非線形CAE勉強会

6月

第２回数理物理・統計力学

2017年

セミナー

6月

国内
国内
国内
国内

受賞等
受賞者氏名（所属機名称関）

第30回バイオエンジニアリン

大谷智仁、足利洋志,

グ講演会 Outstanding Poster
Presentation
第30回バイオエンジニアリン
グ講演会 Outstanding Poster

授賞機関(学会名等)
日本機械学会バイオエ

受賞した

国内・国

時期

際の別

2017 年 12和田成生（大阪大学） ンジニアリング部門

月15日

武石直樹、伊井仁志,

2017 年 12日本機械学会バイオエ

和田成生（大阪大学） ンジニアリング部門

Presentation

月15日

日本機械学会バイオエ2017 年 12

ンジニアリング部門

月15日

第26回業績賞

和田成生（大阪大学）

2017 年武田医学賞

小川誠司（京都大学） 武田科学振興財団

平成 29 年度 科学技術分野の文部科学大臣表彰 若手科学者 片岡圭亮（京都大学） 文部科学省
賞
3.

2017年

第32期 非線形CAE勉強会

国内

2017 年 11
月 13 日2017 年 4月 19 日

国内

国内

国内
国内

国内

メディアへの情報発信、ウェブサイト等での情報公開

名称コンピュータシミュレー
ションによる心臓研究日付
2017年4月

説明

備考

著者所属

京都府立医科大学

その他[大学杉浦 清了(東大大学院特別講義

院特別講義]/UT-heart 研究所)

120

4.

書籍

タイトル

執筆者（所属機関）

書籍名

発行年月

数理モデル研究による腫瘍内不均

新井田厚司、宮野悟（東実験医学 増刊2018 年 2

一性と治療抵抗性への挑戦 新京大学）

Vol.36-No.2 2018月

実験医学 増刊2018 年 2

Vol.36-No.2 2018月

骨髄異形成症候群の病態とクローン進化

小川誠司（京都大学）

４－４．実施体制
1.

業務主任者（受託者（委託先））
役職・氏名

国立大学法人東京大学

医科学研究所教授

宮野悟

E-メールアドレス： miyano@ims.u-tokyo.ac.jp
TEL：03-5449-5615（再委託先）
役職・氏名

FAX: 03-5449-5442

国立大学法人大阪大学

大学院基礎工学研究科

教授

和田

成生

E-メールアドレス：shigeo@me.es.osaka-u.ac.jp
TEL&FAX：06-6850-6170（再委託先)
役職・氏名

株式会社 UT-Heart 研究所

代表取締役会長

久田俊明

E-メールアドレス：hisada@ut-heart.com
TEL：03-3410-0216（再委託先）
役職・氏名

FAX：03-3410-0216

国立大学法人京都大学

医学研究科

教授

小川

誠司

E-メールアドレス：sogawa-tky@umin.ac.jp
TEL： 075-753-9283（再委託先）
役職・氏名

FAX：075-753-9282

学校法人自治医科大学

学長

永井

良三

E-メールアドレス：rnagai@jichi.ac.jp
TEL：
2.

0285-58-7006

FAX：0285-44-5019

業務項目別実施区分
業務項目

担当機関

担当責任者

（１）

東京都港区白金台4-6-1

ヒトゲノム解析センター長

大量シーケンスによるがんの国立大学法人東京大学

宮野悟

個性と時間的・空間的多様性・ 医科学研究所ヒトゲノム解析セ
起源の解明（サブ課題 A）

ンター
京都府京都市左京区吉田近衛町

121

大学院医学研究科教授

国立大学法人京都大学

小川誠司

大学院医学研究科（２）

大阪府豊中市待兼山町1-3

大学院基礎工学研究科

データ同化生体シミュレーシ

国立大学法人大阪大学

和田ョンによる個別化医療支援（サ大学院基礎工学研究科

教授

成生

ブ課題B）
東京都文京区本郷 7-3-1

大学院工学系研究科

国立大学法人東京大学

高木

教授

周

大学院工学系研究科
大阪府吹田市山田丘1-1

大学院情報科学研究科 教授

国立大学法人大阪大学

松田秀雄

大学院情報科学研究科
大阪府豊中市待兼山町1-3

大学院基礎工学研究科

国立大学法人大阪大学

野村

教授

泰伸

大学院基礎工学研究科
大阪府吹田市山田丘1-8

歯学部附属病院

国立大学法人大阪大学

玉川

准教授

裕夫

歯学部附属病院（３）

東京都世田谷区野沢3-25-8

代表取締役会長

心臓シミュレーションと分子株式会社UT-Heart研究所

久田俊明

る基礎医学と臨床医学の架橋

栃木県下野市薬師寺3311-1

学長（サブ課題C）

学校法人自治医科大学（４）

東京都港区白金台4-6-1

ヒトゲノム解析センター長

プロジェクトの総合的推進

国立大学法人東京大学

宮野

シミュレーションの融合によ永井悟

医科学研究所ヒトゲノム解析セ
ンター
3.

経理担当者
（受託者（委託先））
役職・氏名

国立大学法人東京大学

医科学研究所

外部資金戦略チーム

齋藤

主任

利弥

E-メールアドレス：t-gshikin@ims.u-tokyo.ac.jp
TEL：03-5449-5138

FAX： 03-6409-2017（再委託先）

122

事務部

研究支援課

良三

役職・氏名国立大学法人大阪大学
研究協力係

係長

大学院基礎工学研究科
武中 典子

E-メールアドレス：ki-kenkyukyoryoku@office.osaka-u.ac.jp
TEL：06-6850-6142

FAX：06-6850-6145（再委託先）
役職・氏名

株式会社 UT-Heart 研究所
代表取締役会長

久田

研究開発部

俊明

E-メールアドレス：kida@ut-heart.com
TEL：03-3410-0216

FAX：03-3410-0216（再委託先）
役職・氏名

国立大学法人京都大学
経理・研究協力課

医学系研究科外部資金掛長医学・病院構内共通事務部

矢野哲也

E-メールアドレス：a40gaishi@mail2.adm.kyoto-u.ac.jp
TEL：075-753-9483

FAX：075-753-4347（再委託先）
役職・氏名

学校法人自治医科大学
経理第３係

係長

総務部総務経理課髙巣嘉子

E-メールアドレス：uetakeh@jichi.ac.jp
TEL：0285-58-7022

FAX：0285-40-8014

4.

知的財産権の帰属「知的財産権は乙に帰属することを希望する。」

5.

委託契約書の定めにより甲に提出することとされている著作物以外で委託業務により作成し、甲に納入する著作物の有無
（

有

・

無

）

123