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

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

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

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

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

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

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

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

４－３.活動（研究会等） ....................................................................................................................... 85
４－４.実施体制...................................................................................................................................... 88
別添 学会等発表実績
別添 実施計画i１．委託業務の題目「個別化・予防医療を支援する統合計算生命科学」
２．実施機関（代表機関）
機関名

国立大学法人東京大学

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

ふ り が な
課題
代

責任者

表

氏

名

所属部署名
連絡先

みやの さとる生年

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

宮野 悟

月日

※2016 年 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

役職

主任

Fax. 03-6409-2017

E-mail: t-gshikin@ims.u-tokyo.ac.jp

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

大学法人大阪大学、株式会社 UT-Heart 研究所、学校法人自治医科大学、国立大学法人岡山大学と密接に連携し、研究開発を実施する。
用語集
３．

ゲノム

オミクス

生物種にとって必要な遺伝子セット。ヒトの場合は 23 対 46 本の染色体のＤＮＡに含まれる全ての遺伝情報を指す。
生物の体の中にある分子全体、特に、ゲノム、RNA、タンパク質、代謝物質などを網羅的に調べる研究分野。

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

大量シーケンスによるがんの個性と時間的・空間的多様性・起源の解明（サブ課題A）
がんはゲノムの変異によって生ずる疾患で、個体別、空間的（腫瘍の場所）、時間的多様性を持っている。がんの起源とその多様性の理解は、がん治療戦略、がんの予防法と超早期発見にイノベーションを起こし、副作用に優しく個人ごとに効き目のよい薬を創出するための戦略上大変重要である。そのため、次の目標を設定し研究を実施する。
(1) 世界の追随を許さない大規模がんオミクスデータ解析システムの開発
平成 28 年度は、大規模がんオミクスデータ解析システム（Genomon 及び SiGN）を発展させ多機能化を図り、スーパーコンピュータ「京」（以下、「京」）でそのシステムの構成を確認する。Genomonについては、I/O に関する詳細設計コデザインに寄与する。また、システムの開発の中で、がん 500
試料を使って「京」で各機能を確認する。
(2) がんの起源と集団内多様性を理解するためのシステムの開発
平成 28 年度は、システム開発のための基礎データの獲得と基礎となるモデルを構築することを目標にし、複数多部位エクソーム等の解析を実施し、数理モデルを検討する。
(3) 疾患に関連する胚細胞バリアントの効率的な解析システムの開発
平 成 28年 度 は 、 胚 細 胞 バ リ ア ン ト の 効 率 的 な 解 析 シ ス テ ム の 設 計 を 行 う 。 そ し て 、 ゲ ノ
ムシーケンスデータに基づいた高感度な胚細胞バリアント検出法の開発を開始する。平
成 29年 度 ま で こ の 検 出 法 の 開 発 を 続 け 、 シ ス テ ム を 完 成 さ せ る 。
2.

データ同化生体シミュレーションによる個別化医療支援（サブ課題B）
ポスト「京」を活用した大規模生体物理シミュレーションと生体計測データとを様々なレベルで同化・融合させることにより、実測データを重視する医療に受け入れられる計算機シミュレータを開発し、シ2

ミュレーションで得られる物理情報を活用した個別化医療支援の実現を目指す。そのため、次の目標を設定し研究を実施する。
(1) 生体シミュレータの開発および高度化
平成 28 年度は、個別化医療支援に供する生体シミュレータの計算モデルを構築し、ポスト「京」での実行に向けたプログラムを開発する。具体的には、生体機能を生み出す循環系、筋骨格系、口腔系の動態を再現する物理の原理に基づく計算モデルを構築し、それらを数値的に解析するための「京」による大規模計算コードの移植と改良を実施する。また、これらを統合する脳神経系を理解する上で重要となる全脳循環シミュレータの開発に向けたモデル構築を行う。
(2) データ同化手法の開発と個別化医療支援への応用
生体シミュレータと計測データを融合・同化させるための手法を開発する。平成 28 年度は、そのための計測データや医用画像データを収集するとともに、「京」による大規模ライブセルイメージ
ング解析技術を開発する。
3.

心臓シミュレーションと分子シミュレーションの融合による基礎医学と臨床医学の架橋（サブ課題C）
タンパク機能（遺伝子）の変化が最終的に臓器レベルにどのような変化を及ぼすかを物理的に解明・予測することができれば、その応用は無限に広がる。ポスト「京」のパワーを活かし、本サブ課題では、世界で初めて心臓シミュレーションと分子シミュレーションを融合させた真のマルチスケールシミュレーションを実現する。これにより計算科学の歴史に新たなマイルストーンを築くとともに、現実の医療・創薬にブレークスルーをもたらす。
本サブ課題では具体的に二つの目標を掲げている。
目標 1. 心臓病、特に「心不全」の根本的解明と最適治療を可能とするマルチスケール心臓シミュレータの開発を目的とする。モデル化は心不全研究において得られた遺伝子レベルの変化、発現調
節、細胞内情報ネットワークに関する膨大な情報、臨床情報を効果的に取り入れて推進する。
平成 28 年度は、サルコメア分子モデルの開発と臓器モデルの高度化を行う。具体的には基礎医学研究に基づいてサルコメアモデルを開発する。連成法の「京」による検討も行う。また、心不全治療と外科手術の急性効果を検討する。
目標 2. 重点課題 1 において進められるイオンチャンネルの分子シミュレーションと UT-Heart を融合することで、候補物質の分子構造と遺伝子情報のみから不整脈のリスク評価を可能とするシス
テムの完成を目標とする。
平成 28 年度は、イオンチャネル分子モデルおよび連成手法を開発する。具体的には、膜電位に依存し開閉するイオンチャネル分子モデル、並びに UT-Heart との連成法を開発する。
4.

プロジェクトの総合的推進
プロジェクト全体の連携を密としつつ円滑に運営していくため、運営委員会や研究会の開催等、参画各機関の連携・調整にあたる。
特に、プロジェクト全体の進捗状況を確認しつつ計画の合理化を検討し、必要に応じて調査あるいは外部有識者からなる外部諮問委員会を開催して意見を聞くなど、プロジェクトの推進に資する。ニュースレター（電子版）を発行する。ホームページを運営する。シンポジウムを開催し研究の進捗と連携を3

推進する。プロジェクトで得られた成果については、積極的に公表し、今後の展開に資するが、UT-Heartのコードについては知財の観点から公開は行わない。また、ソフトウェアの普及のための教育セミナーを開催するとともに、およびそれに伴った人材育成活動を行う。このために広報を担当する人材とリサーチアドミニストレータを雇用する。

4

４－２．実施内容（成果）1. 大量シーケンスによるがんの個性と時間的・空間的多様性・起源の解明（サブ課題A）(1) 世界の追随を許さない大規模がんオミクスデータ解析システムの開発
Genomon2（改定版の番号を付している）は、偽陽性の変異を効果的にフィルターするプログラム、GenomonMutationFilter (https://github.com/ Genomon-Project/ GenomonMutationFilter)、多数のコ
ントロール検体を効率的に利用して偽陽性を除くプログラム、EBFilter(https://github.com/Genomon-Project/EBFilter)などを備え、特に複雑な配列パターンを有する全ゲノム領域においても高度に正確な変異検出を可能とした。また、 GenomonSV(https://github.com/Genomon-Project/GenomonSV) により、大規模な構造変異だけではなく、数十塩基から数百塩基の中程度のサイズの構造変異の検出が可能となった。さらに、ruffus, drmaaなど最新のジョブ管理ツールを効率的に利用したGenomonPipeline(https://github.com/Genomon-Project/GenomonPipeline)により多数検体の効率的な解析を容易にした。
このGenomonの移植に関して、I/Oに関する詳細設計コデザインの一つとして新たに「京」向けのバー
チャル・グリッドエンジンを開発し、「京」でそのシステムの構成を確認した。これにより多数検体解析
手法の実装が進展した。
もう一つの大規模がんオミクスデータ解析システムである SiGN に関して、「京」互換機上で高度なネットワーク推定機能を実証した。
本解析システムを使用した成果として、平成 27 年度にがんの個別化医療・予防の実現のために大量シーケンスデータ解析が必要であることを示したことに続いて、平成 28 年度は、大規模がんオミクスデー
タ解析システムの一つである Genomon およびその利用環境 VGE の開発を継続して行った。その開発の中で、2 つの重要な成果を得た（詳細は D）の項目を参照）。平成 27 年度に行った成人 T 細胞白血病リンパ腫(ATL)の全ゲノムシーケンス解析（49 検体）およびオミクス解析（400 以上の検体で、ATL 研究
としては世界最大）により ATL の遺伝子変異の全容が解明されたが(Kataoka K, et al. Nature Genetics.
47(11):1304-1315, 2015)、その中の一つに免疫チェックポイントパスウェイの変異があった（PD-L1 遺伝子の構造異常など）。それを手がかりに米国データベース TCGA の 10,210 検体のデータを解析し、またゲノム編集技術によるバリデーションを行い、がん細胞が免疫細胞を回避するメカニズムを解明した(Kataoka K, et al. Nature. 534(7607):402-406, 2016)。そこでは、これまでショートリードシーケンス
解析の盲点であったミッドレンジ(20-500bp)の構造異常検出を、独自の方法で Genomon が高精度に行いブレークスルーした。これにより免疫チェックポイント阻害剤の有効性をゲノム異常の検出でスクリーニングする可能性が示された。第二の成果は、時系列のゲノムシーケンス解析により、骨髄異形成症候
群 (MDS) に お け る ク ロ ー ン 進 化 を 解 明 し た こ と で あ る (Makishima H, et al. Nature Genetics.
49(2):204-212, 2017)。急性白血病(AML)を起こす 2 ステップの遺伝子異常のパターンが明かになった（2,250 検体データを解析）。即ち、MDS から AML へ進行する際に、2 つのタイプのゲノム異常が細胞に蓄積し異常クローンが増大すること、そして、低リスク MDS に「タイプ 2 ゲノム異常」が起こると白血病リスクの高い状態（高リスク）に移行し、「タイプ 1 ゲノム異常」が起こると急性白血病が起こることが判明した。さらに一方で、SF3B1 変異陽性症例は、タイプ 1、タイプ 2 異常を持たない場合、ほとんど白血病を発症しないことも観察された。
このような世界最高レベルの成果が世界の強豪グループを一気に追い抜く形で得られたのは、普通に5は扱えない規模の大規模シーケンスデータを対象とし、スーパーコンピュータを駆使して大規模データ
解析を行える世界の追随を許さないシステムが構築されつつあることによると考える。
A) Virtual Grid Engine の開発(a) 概要重点課題 2 のターゲットアプリケーションである Genomon2 の「京」/ポスト「京」上での動作のために、Virtual Grid Engine (VGE)の開発を実施した。VGE は非グリッドエンジン環境に置いて、MPI のプロセス空間を疑似的なグリッドエンジン環境としてユーザアプリケーションに提供するミ
ドルウェアである。
ゲノムシーケンスや RNA シーケンスデータの解析ではパイプラインが一般的であり、本体および内部利用の外部ソフトウェア等の更新が頻繁（~10 回/年）なバイオインフォマティクス分野において、ソフトウェアの個別 MPI 化は発見に至るまでの時間的競争という生命科学的優位性の観点から現実的な対応ではない。グリッドエンジンのポスト「京」への搭載については検討中（後述）であるが、「京」を含めた HPCI 計算機群や他のペタフロップス級スーパーコンピュータ上での柔軟な運用のためには、必要不可欠なソフトウェアである。平成 27 年度に実施した VGE の設計（図 1.1）に基づき、本年度は VGE の実装および Genomon2 の移植・動作検証を実施した。
VGE は、(i) python による実装、(ii) mpi4py モジュールを通した MPI プログラム、(iii) マスター
ワーカーモデルによるグリッドエンジンライクなアプリケーションタスク実行、の 3 つの特徴を備えている。
実装言語に選択した python は、現在、バイオインフォマティクス分野で利用される主要プログラ
ミング言語の一つである。その他の利用言語としては R、perl やシェルスクリプトが上げられる。
一般に HPC 分野で利用される Fortran や C などの高級言語はこの分野での利用はごく一部に限られている。Python で実装することにより、VGE のようなミドルウェアの同分野における普及に大きな効果が期待できる。

6

図 1.1: バーチャル・グリッドエンジンを用いた Genomon2 の多数検体解析実施概念図
mpi4py モジュールは python から MPI ライブラリを利用可能にするパッケージであり、「京」にもインストールされている。同モジュールを利用することにより、C/Fortran ラッパーを用いずに
python から直接 MPI 並列プログラムを作成・実行することができる。複数言語の混成は環境依存性を高め、移植性を著しく下げる傾向があるため本ライブラリ利用効果は高い。
グリッドエンジンの並列処理は、システムのリソース（CPU や計算ノード）に空きができ次第、順次割り当て・実行が行われる。これは、並列処理が指定されたタスク内のプロセス同士にプロセス間通信が必要ない、独立プロセスの集合だからである。MPI プログラムでこのような処理を実現するのに最も適したアルゴリズムがマスターワーカーモデルである。VGE では最初に使用するリソース全体に対して一つの MPI プロセスを起動する。このプロセスのうち特定ランク（通常はランクゼロ）をマスタープロセスとし、残りのリソースをワーカープロセスに割り当てる。アプリケーシ
ョン側からは、VGE の提供する vge_task 関数を qsub/pjsub の代わりに使用し、バッチジョブシス
テムではなく VGE へジョブを投入させる仕組みとなっている。
(b) 実装図 1.2 に VGE のフローチャートを示す。VGE は主として 3 種のプロセスで構成される。VGE は
起動時に指定されたプロセス数の MPI プログラムとして動作を開始する（VGE 開始）。初期設定を読み込んだのち、マスターランクとそれ以外で動作が区分される。マスターランクの場合、pythonの multiprocessing モジュールを利用して 2 つの子プロセスが起動する。
7

VGE開始フロー3

前処理/MPI処理
workerランク

通信処理

フロー1

MPIプロセス並列処理

フロー2

masterランク
multiprocessing並列

MPIジョブコントロール

Genomonジョブコントロール

前処理
masterへジョブリクエスト
前ジョブ結果送信

前処理

前処理
共有オブジェクト

MPI通信

NO

データ参照

ジョブ情報

masterからデータ受信

MPI終了？

新ジョブ確認/共有データ更新

新ジョブ・情報作成
YES

YES

ソケット通信要求内容処理ジョブデータNO

ジョブある？

受信データ処理

MPI通信

Pythonデータ

Genomonからデータ受信
内部プロセス間
データ共有

workerからデータ受信

NO

Genomonと接続OK？

Pythonmultiprocessing
manager

GenomonPipeline

masterへデータ送信

MPI
終了フラグ

YES

データ参照

新ジョブ？
実行ジョブ終了確認？
VGE終了？

宅配ボックスライク

workerへデータ送信

共有データ更新/送信データ作成

ジョブ実行
終了ジョブ確認/データ更新
ジョブステータス処理NO

Genomonへデータ送信
データ参照

ソケット通信YESMPI終了？

VGE終了？

YES
worker終了

Pythonデータ

NO
待機MPIコントロール終了
VGE終了

図 1.2: バーチャル・グリッドエンジンのフローチャート（Genomon2 実行の場合）
一つは MPI ジョブコントロールプロセス、もう一つはアプリケーションジョブコントロールプロ
セス（Genomon ジョブコントロール）である。
VGE を利用するアプリケーションは、vge_task 関数を通して VGE へジョブを投入する。投入されたジョブを受け付けるプロセスがアプリケーションジョブコントロールプロセスである。このプ
ロセスはマスターランクに常駐し、アプリケーションとソケット通信によってジョブ情報を送受信する。すなわち、アプリからのジョブ実行リクエストを受信しジョブプールへデータを追加、ジョブ終了時に終了ステータス等の情報をアプリケーションへ送信する。MPI ジョブコントロールプロ
セスも同様にマスターランクに常駐する。このプロセスはアプリケーションジョブコントロールプ
ロセスとワーカープロセスの 2 種のプロセスと通信を行う。アプリケーションジョブコントロール
プロセスとの通信には python multiprocessing モジュールの提供する共有オブジェクトを利用している。共有オブジェクト内に作成されたジョブプールからキューイングおよび投入済みジョブの状
態更新を行う。
マスターランク以外のランクがワーカーである（図 1.2 青線左部分）。ワーカーは MPI ジョブコン
トロールプロセスと MPI 通信によって情報の送受信を行う。VGE 起動後、各ワーカーは待機状態で起動され、MPI ジョブコントロールプロセスからのジョブ実行要求を待つ。ジョブ情報を受信す
ると、ワーカーは指定されたジョブをワーカープロセス内にサブプロセスとして起動する。その後、サブプロセスの終了まで待機し、ジョブステータスを MPI ジョブコントロールプロセスへと送信し8

再び待機する。
アプリケーションの実行が全て終了した後、VGE の終了処理が実行され全プロセスへと終了シグ
ナルが送信される。
(c) 実行テスト
VGE の動作確認のため、Genomon2 を、サンプルデータを用いて「京」上で実行した。Genomon2の VGE 環境への対応は stage_task.py ファイルの修正のみである。詳細に関しては Genomon2 の章で述べる。データサイズはおよそ 61MB、mutation call および SV（Structural Variation: 構造変異）の検出に対応したサンプルデータである。MPI プロセス数を 2～81（マスター1+ワーカー数
なので最少プロセス数は 2 となる）と変化させた場合の解析時間を表 1 に示す。データサイズが小さいため加速率の伸びが低いが、ワーカー数 5 程度までは 8 割弱の加速率を達成している。また、シリアル版とワーカー数 1 の差が 12 秒程度であり、VGE によるオーバーヘッドも小さなものであることが確認できた。
表 1.1:「京」上での Genomon2 実行時間
ワーカー数

MPI
プロセス数

実行時間(秒）

加速率

並列化効率(%)

1

21130.54

1

100

2

3

596.17

1.90

95.0

5

6

296.43

3.81

76.2

10

11

203.25

5.56

55.6

20

21

169.04

6.69

33.5

30

31

139.51

8.10

27.0

40

41

130.05

8.69

21.7

50

51

131.85

8.57

17.1

60

61

126.498.94

14.9

129.65

8.72

10.9

80

81
シリアル版

1142.30

B) Genomon2 の移植作業およびコデザイン(a) 概要重点課題 2 のターゲットアプリケーションである Genomon2 の移植作業を実施した。MPI 環境への対応は VGE 側で実施しているため、グリッドエンジン環境部分を VGE に変更することが作業内容となる。
コデザインに関しては、Genomon2 が移植作業中ということから計算資源の 8 割以上を使用する
BWA（Burrows-Wheeler Aligner[1]）のホットスポット解析、およびボトルネック部分の SIMD 化
等による高速化検証と必須演算に関する回路・命令セットの開発要求を設計チームへ報告した。また、Genomon2 は他ターゲットアプリに比べ I/O 負荷の高いアプリケーションである。パイプライ
ン中のファイル I/O に関しても調査結果を富士通に報告済みである。また、グリッドエンジンの機能に関して、ポスト「京」のバッチジョブスケジューラへの組み込み可否についてスケジューラ開9

発チームと検討中である。
(b) Genomon2 の移植
前節で述べたように、本移植に関しては Genomon2 自身の大幅な改変を伴うことなく、「京」をはじめとした MPI 環境での動作を目標に設定している。その意味で、Genomon2 自身の移植作業は
SGE/UGE 等のグリッドエンジン環境から VGE 環境への対応部分のみとなっている。Genomon2
はワークフローマネジメントシステム（WMS）として ruffus を採用しており、パイプライン中の各タスクからはジョブ投入のための専用関数 stage_task を通してバッチジョブシステムへとジョブをサブミットするシステムとなっている。よって、ここでは stage_task で qsub/pjsub コマンドを発行する部分を vge_task（VGE の提供するジョブ投入専用 python 関数。VGE の章を参照。）を使用す
るように数行書き換えるだけで作業が完了した。
(c) コデザイン
まず初めに BWA のチューニングに関して報告する。BWA は Genomon2 の中で 8 割以上の計算機
リソースを使用する主要演算ソフトウェアである。処理内容は、次世代シーケンサから出力された
ショートリードデータ（FASTQ 形式）を標準ゲノム配列（リファレンスデータ）と照合し、位置を決定すること（マッピング）である。BWA は全ゲノム（Whole Genome Sequence: WGS）のマッ
ピングツールのデファクトスタンダードである。BWA の高速化で注意すべき点は、他のアプリケー
ションで良く見られるホットスポットに相当する関数が無く、3～4 か所の演算に 10%～30%程度の負荷が分散している[2]ことである。我々はまず、小規模テストデータを用いて「京」と Intel における BWA の実行時間および簡易プロファイリングを実施した。図 1.3 上にあるように、Intel マシ
ンでは数分で終わる計算が「京」で数十倍の時間が掛かっていることが判明した。
「京」で極端に実行時間が遅くなっている部分は worker2 関数であった。worker2 は Smith-Warterman アルゴリズ
ムを用いてリードの詳細なマッピングを実施する部分であり、内部で ksw_u8 関数を使用する。本関数には Intel x86_64 アーキテクチャ用 SSE2 イントリンシック命令が使用されており、BWA の移植の際に我々が開発した代替命令（mysse2）を使用している個所である。

10

図 1.3: BWA のプロファイリング
（上:Intel および「京」での実行時間比較、下:Intel Xeon 上での実行時間内訳）
mysse2 は BWA で使用される SSE 命令を、非 Intel 環境で動作可能にすること目的に開発したライブラリであり、高速化等の処理は施されていない。調査の結果、mysse2 が多大な時間を消費していることが判明したが、「京」上には ksw_u8 に相当する演算回路および命令が存在しないことが判明した。そこで我々は、SW を実装している ksw_align2 関数を書き換え、ksw_u8 を使用しない SW
アルゴリズムのナイーブな実装として ksw_i32 関数を開発した。本関数を用いて実行時間を測定したところ、FX100 上で mysse2 を使用した場合に比べて 2.6 倍程度の高速化を達成した（図 1.4）。
一方、Haswell でのオリジナル実行時間との比較からわかるように、ksw_i32 は ksw_u8 を使用す
る演算に比べまだかなり性能面での差が大きい。
SW アルゴリズムで要求される演算は 8 ビット整数
11

の SIMD である。Sparc64VIIIfx および富士通コンパイラには現在のところ該当する演算回路および命令が無く、これ以上のソフトウェア的な高速化には限界があることが判明した。よって、8 ビッ
ト整数に対する演算回路および命令セット（最低限、max 演算）に関しての開発を要望として開発チームへ報告した。

図 1.4: ksw_i32 を使用した BWA の実行時間
次に、I/O の調査結果について報告する。Genomon2 ではパイプライン中のタスク間での情報のやり取りは基本的にファイル渡しとなっている。これはバイオインフォマティクス分野でのパイプラインで一般的な方法である。Genomon2 では多数検体の解析を目標にしている。そのため、プロセス間で再利用可能なサイズの大きいファイル等はパイプライン実行後もシステムに残す方が I/O シ
ステムの負荷軽減およびパイプライン自体の実行効率向上に繋がる。このような観点から、
Genomon2 における主要入出力ファイルとその利用タイミング、使用後のファイル消去可能性等を分析した（表 1.2）。
最後に、ジョブスケジューラへの要望である。Genomon を含むパイプラインはグリッドエンジンを念頭に作成されることがほとんどである。今回、VGE の開発を通してグリッドエンジンに必要な機能を整理した。基本的にはプロセス間通信を必要としない多数プロセスの一括実行機能と複数ジ
ョブ（パイプライン中の多段のタスクは、スケジューラからは独立のジョブに見える）の依存関係
処理、および計算ノードからのジョブ投入機能が基本機能である。これに加えて、「京」等の大規模計算環境ではパイプライン中のジョブはできるだけ連続に処理できるようにすることが実用上の観点から必要となる。これらの要望・機能はジョブ充填率等の面で既存のジョブスケジューリング方法ではそのまま実装することは難しい。現在、スケジューラチームとの検討を実施中である。

12

表 1.2:Genomon2 ファイル I/O 一覧
ファイ
1個あたりの
ファイル名
ファイル数
前処理 Task1A Task3A Task4A Task5A Task1B Task6A Task7A
備考
ル番号
ファイルサイズ$INDIR/testdata/normalreference*.bam
ファイル23と同程
1#list（シンボリックリンク先）
度$INDIR/testdata/normalreference*.bam.b
ファイル25と同程
2 ai
#list
度
（シンボリックリンク先）
3 GRCh37.fa.amb
1
100KB
◇
これら5ファイルはTask4Aのリファレンス用入力ファイ
4 GRCh37.fa.ann
110KB◇
ルである。Task4Aでは1708プロセス動作し、これらの5 GRCh37.fa.bwt
1
3GB◇
全プロセスから参照される。また、検体が異なる実行
6 GRCh37.fa.pac
1
750MB
◇
においても共通に利用される。（Genomon2の実行に
7 GRCh37.fa.sa
1
1.5GB
◇
おいてユーザーIDは同じでジョブIDが異なる状況でア
8 $INDIR/sample.csv
1KB◇◇ ・・・入力ファイル
9 control_panel.txt
1
KB◆◆ ・・・Task7A 終了後も保存される出力ファイル
10 $INDIR/testdata/${TYPE}.bam
2
160GB
◇
★ ・・・ファイルの生成と削除を行うタスクが同じ場合の一時作業用出力ファイル
11 sample_${TYPE}/*_1.fastq
2
250GB
◆
◇
▲ ・・・ファイルの生成と削除を行うタスクが異なる場合の一時作業用出力ファイル
12 sample_${TYPE}/*_2.fastq
2
250GB
◆
◇
× ・・・一時作業ファイルが削除されるタイミング
13 single_end_output.txt
2
KB◆
S ・・・シンボリックリンクを作成（カッコ内はシンボリックリンク先を表す）
14 unmatched_first_output.txt
2
KB◆
TYPE・・・normal または tumor
15 unmatched_second_output.txt
2
KB◆* ・・・0-9 までのいずれかの数値、*の個数は桁数を意味する
1708
293MB
#list・・・$INDIR/sample.csv の[controlpanel]に含まれるリストの数
16 sample_${TYPE}/1_****.fastq_split▲◇×(=854×2) (=250GB/854)
1708
293MB
17 sample_${TYPE}/2_****.fastq_split▲◇×(=854×2) (=250GB/854)
18 fastq_line_num.txt
1
KB
◆
◇
1708
585MB
19 sample_${TYPE}_****_bwa.sam★×(=854×2) (=500GB/854)
1708
187MB
20 sample_${TYPE}_****.sorted.bam▲◇×(=854×2) (=160GB/854)1708
21 sample_${TYPE}_****.sorted.bam.bai
10MB～
▲
◇×
(=854×2)
1708
187MB
22 sample_${TYPE}_****.sorted.bam.tmp★×(=854×2) (=160GB/854)
23 sample_${TYPE}_markdup.bam
2
160GB
◆
◇
24 sample_${TYPE}_markdup.metrics
2
KB◆◇
25 sample_${TYPE}_markdup.bam.bai
2
10MB～
◆
◇
26 sample_${TYPE}_markdup.bam.md5
2
KB
◆
S(ファイ
27 reference*/reference*.markdup.bam
#list
◇
ル1)
S(ファイ
28 reference*/reference*.markdup.bam.bai
#list
◇
ル2)
$database/GRCh37/GRCh37_noScaffold
29
1
KB
◇
_noDecoy.interval_list
30 sample_tumor.fisher_mutations.**.txt

30

KB

★×

sample_tumor.realignment_mutations.**.t
31
xt

30

KB

★×

32 sample_tumor.indel_mutations.**.txt

30

KB

★×

33 sample_tumor.breakpoint_mutations.**.txt

30

KB

★×

sample_tumor.simplerepeat_mutations.**.
txt

30

KB

★×

30

KB

★×

30

KB

▲

1

KB

34

35 sample_tumor.ebfilter_mutations.**.txt
sample_tumor_mutations_candidate.**.hg1
369_multianno.txt
sample_tumor_genomon_mutations.result.t
37
xt

◇×
◆

[1] Li H. and Durbin R. Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics,
25:1754-1760, 2009.
[2] Ahmed N., et al. Heterogeneous Hardware/Software Acceleration of the BWA-MEM DNA Alignment
Algorithm. IEEE/ACM International Conference on Computer-Aided Design (ICCAD), pp.240-246, 2015.
C) ネットワーク解析ツール(SiGN)開発の進捗
遺伝子ネットワーク推定ソフトウェアである SiGN は非線形回帰によるベイジアンネットワーク推定を用いた SiGN-BN、時系列データを対象として状態空間モデルを用いた SiGN-SSM、L1 正則化法を用いた SiGN-L1 から構成されるが、そのうちの SiGN-BN は適用するデータのサイズによって使用するアルゴリズムを使い分ける。SiGN-BN に実装されているアルゴリズムのうち、Para-OS と呼ばれるアルゴリズムは、適用可能な遺伝子数が 30 程度と限られるが非線形回帰に基づくスコア上で数13

学的な意味での厳密な最適解を求めることができる高精度な遺伝子ネットワーク推定アルゴリズムである。このアルゴリズムの改良を「京」互換機である Fujitsu FX-10 上で行い、離散スコアを用いて36 遺伝子までの遺伝子ネットワーク推定が可能なことを実証した (Honda et al., 2016)。これまでの報告されている最大サイズは 33 変数（遺伝子）である。この実装は、ベイジアンネットワークの構
造探索アルゴリズムの研究に目的を絞ったものであるため、遺伝子やノンコーディング RNA のマイクロアレイや RNA-seq などの計測データを用いた遺伝子ネットワークの推定には直接使用できない。
そのため「京」および「京」互換機である FX-100 を用いて、非線形回帰モデルスコアの実装を行った。また同時にアルゴリズムの改良を試みた。新アルゴリズムを FX-100 上で実行したところ、864
ノード 27648 コアを用いて 34 遺伝子の遺伝子ネットワークを 21 時間半で計算可能なことを実証した。
また必要なメモリーも大幅に減少しており（3.5GB→0.7GB、28 変数 324 計算ノード利用時）、計算
ノードあたりのメモリー量が増えない状況においてもより大規模な最適遺伝子ネットワーク推定が可能になった。当該アルゴリズムは「京」互換機およびヒトゲノム解析センターのスパコン向けの実装
が SiGN-BN のウェブページにて公開中である。
D) Genomon による大量シーケンスデータの解析によるがんの個性と時間的多様性の解明(a) がん細胞が免疫から逃れるメカニズムの解明
i. 研究の背景
がんは今や国民の半数が罹患し、死因の 1/3 を占める深刻な疾患となっており、超高齢社会を迎えて現在、増加の一途を辿っている。一方、その治療成績は、集学的な治療の努力にも関わらず伸び悩んでおり、近年開発が進んでいるいわゆる「分子標的薬剤」も、一時的にはしばしば著効を示すものの、これらの薬剤に耐性を有するがん細胞の出現によって、患者さんの多く
はやがて再発することから、がんの治療成績への寄与は、なお、不十分と言わざるを得ない。
一方、近年我が国の研究者による発見を契機として開発された抗 PD-1 抗体および抗 PD-L1
抗体が、悪性リンパ腫や、肺がん、悪性黒色腫をはじめとする多くのがん種において、顕著な臨床効果を示し、末期のがん患者においてもしばしば治癒をもたらすことが報告されている。
一般に生体には、本来がん細胞を除去するための免疫監視のメカニズムが備わっていて、がんの発症が抑制されていると考えられている。こうした、いわゆる「がん免疫」は、様々な免疫担当細胞によって担われているが、このうち、直接的にがん細胞の殺傷に関わる細胞障害性 T細胞（CTL）は、がん細胞に特異的に発現する多様な分子を通じてがん細胞を認識し、これを殺傷することにより、がん免疫を担っている。PD-L1 は、がん細胞に発現する「免疫チェック
ポイント」を構成する一群の分子の一つで、がん細胞を攻撃する CTL が発現する PD-1 分子に作用して、CTL を抑制することにより、こうした CTL からの攻撃を逃れることができると考えられる（図 1.5）
。抗 PD-1 抗体および抗 PD-L1 抗体は、がん細胞によるこの免疫回避の仕組みを無効にすることによって、免疫担当細胞が再びがん細胞を認識し、これを攻撃することを可能にすることによって、抗腫瘍効果を発揮すると考えられる。このように抗 PD-1 抗体および抗PD-L1 抗体をはじめとする免疫チェックポイント阻害薬では、生体にそもそも備わっているがん免疫を再度活性化することによって、抗腫瘍効果を得る、極めてユニークな治療法といえる。
また、免疫システムは、患者の体内にある極めて多様ながん細胞を認識して殺傷できることが、

14

通常の分子標的薬と異なった大きな特色である。

図 1.5: がん細胞の免疫回避の仕組みと免疫チェックポイント阻害剤の標的
しかし問題点として、免疫チェックポイント阻害薬は特定のがん種では極めて高い奏功率が得られるものの、他のがん種では、効果の得られる患者の割合が大変低く、同治療の費用が極めて高額におよぶことを考えると、こうしたがん種については全ての患者に投与することが困難である点が挙げられる。従って、治療効果を予測する優れたマーカーを開発することにより、
例え頻度は低くても、治療効果が期待できる患者さんを正確に同定して、高い治療効果を得る
ことが大変に重要であるが、現在までに、治療効果を高精度に予測できるマーカーは知られていなかった。また、こうしたマーカーを見出すためには、そもそもがん細胞がどのようにして
PD-L1 分子を発現して免疫を回避するのか、また、なぜ治療効果を示すがんと示さないがんがあるのか、に関する理解が不可避であるが、そうしたメカニズムについては多くが不明のままであった。
ii. 研究方法と結果(ア) がんにおける PD-L1 遺伝子のゲノム異常の発見
今回の発見の契機となったのは、同研究グループによって昨年報告された、成人 T 細胞白血病リンパ腫（adult T-cell leukemia lymphoma: ATL）のゲノム解析研究である。ATL はヒト
T 細胞白血病ウイルス 1 型（human T-cell leukemia virus type-1: HTLV-1）の感染によって生じる血液がんの一つで、日本、特に西南日本を中心に好発する、非常に悪性度の高い白血病である。既存の抗がん剤では十分な治療効果を得ることが難しく、同種造血幹細胞移植を除いて根治的な治療手段は知られていない。同研究グループは、昨年の研究に続いて、49 例の ATLの全ゲノム解析を行った結果、13 例（27%）において PD-L1 遺伝子の後半部分を共通に侵す
ゲノムの構造異常を同定した（図 1.6）。

15

図 1.6: ATL における PD-L1 のゲノム異常
これらの構造異常には、染色体の欠失や逆位、転座、重複など様々なタイプの異常が含まれるが、全ての例で、蛋白質に翻訳されない「3′非翻訳領域」と呼ばれる領域の欠損が見出された。さらに、興味深いことに、これらの PD-L1 遺伝子のゲノム異常を持つ症例の全例で、PD-L1
遺伝子の発現の顕著な上昇が認められた（図 1.7）。

図 1.7: ATL における PD-L1 ゲノム異常と発現の関係

PD-L1 遺伝子の 3′非翻訳領域の異常を認めた症例のうち、約半数では PD-L1 蛋白の構造が保たれていたが、残りの例では PD-L1 蛋白が途中で切断されていた（図 1.8）。

16

図 1.8: PD-L1 ゲノム異常を持つ ATL 患者における PD-L1 蛋白の構造
実際、病理診断でもしばしば用いられている免疫染色という方法を用いて、がん細胞のPD-L1 の発現を検討したところ、構造異常のない PD-L1 蛋白を発現するがん細胞は、PD-L1の前半部分（N 末端）、後半部分（C 末端）を認識する抗体のいずれでも検出されるが、後半部分を欠失した PD-L1 蛋白を発現する例では、前半部分を認識する抗体では認識されるものの、後半部分を認識する抗体では認識できなかった（図 1.9）。一方、後半部分が切断された
PD-L1 蛋白も、PD-L1 蛋白の機能に重要な部分（細胞外および膜貫通ドメイン）は保たれており、正常の PD-L1 と同様に機能することが確認されたことから、後半部分を欠失した PD-L1蛋白を発現するがん細胞も免疫の監視から逃れることができると結論した。

図 1.9: PD-L1 ゲノム異常を持つ ATL 患者における免疫染色(イ) 様々な悪性腫瘍における PD-L1 遺伝子異常の解析
こうした一連の ATL に関する解析から、PD-L1 遺伝子の 3′非翻訳領域の異常を介した遺伝
子発現上昇は、ATL のみならず、様々ながん種で重要な免疫回避のメカニズムとなっている可能性が予測できるので、米国のがんゲノムアトラス（The Cancer Genome Atlas（TCGA））に登録されている、33 種類の悪性腫瘍からなる 10,210 例のがん試料の遺伝子解析データを用
いて PD-L1 遺伝子の 3′非翻訳領域の異常を探索した。その結果、肺がん、胃がん、食道がん、大腸がん、腎がん、膀胱がん、子宮頸がん、子宮体がん、頭頸部がん、悪性黒色腫、B 細胞リ
17

ンパ腫を含む、12 種類のがん種を含む、計 32 症例において、PD-L1 遺伝子の 3′非翻訳領域の異常が同定された（図 1.10）。本異常は、特に B 細胞悪性リンパ腫（8%）と胃がん（2%）
で高頻度に認められ、さらに興味深い点として、TCGA 症例の一部ではウイルスゲノムのPD-L1 遺伝子領域への組み込みによって、PD-L1 遺伝子の 3′非翻訳領域の異常が生じていた。
ウイルス感染は、子宮頸がん、頭頸部がん、肝がん、胃がんなど、様々ながんの原因になることが知られているが、今回の解析は、一部のがんでは、ゲノムへのウイルスの組み込みによる3′非翻訳領域の異常が PD-L1 遺伝子の発現上昇の原因になっていることを明らかにした。

図 1.10: 様々な悪性腫瘍における PD-L1 ゲノム異常(ウ) ゲノム編集技術を用いた 3′非翻訳領域の機能の検討
一方、ウイルスの組み込みの有無に関わらず特筆すべき点として、これらの 32 例の全例でPD-L1 遺伝子の高度な発現上昇が認められたことが挙げられる。これらのうち、12 例については各々のがん種の中で最も高い発現を示していた（図 1.10）。このことは、ATL と同様に、これらのがん種においても、PD-L1 遺伝子の 3′非翻訳領域の異常が遺伝子発現の調節に大変重要であることを示唆している。そこで、これを検証するために、CRISPR-Cas9 システムという最先端のゲノム編集技術を用いて、ヒトおよびマウス由来の様々な細胞株について、

PD-L1 遺伝子の 3′非翻訳領域の異常（欠失ないし逆位）を導入して、それが PD-L1 遺伝子発現に及ぼす影響を調査した。その結果、PD-L1 遺伝子の 3′非翻訳領域の異常を導入した細胞では、いずれの細胞においても顕著な PD-L1 発現上昇が認められた。このことから、3′非翻訳領域の異常が、一連のがんで PD-L1 遺伝子の発現上昇の原因であることが確認された（図1.11）。実際、3′非翻訳領域を欠いた PD-L1 の転写産物は、正常の PD-L1 転写産物と比較して、細胞内で安定に存在することが確認され、このことは、正常の細胞では、その 3′非翻訳領域を介した PD-L1 の転写産物の安定性の制御を通じて遺伝子発現レベルが調節されていること、
また、3′非翻訳領域の異常を有するがん細胞では、それが破綻を来す結果、転写産物が安定化することによって、PD-L1 の遺伝子発現の上昇が生じていることを示唆している。

18

図 1.11: ゲノム編集を用いた PD-L1 ゲノム異常の導入により PD-L1 高発現が誘導される(エ) 3′非翻訳領域の欠失による PD-L1 遺伝子の発現上昇は、がん細胞の免疫回避と増殖を促進
する
最後に、3′非翻訳領域の欠失により、生体内でがん細胞にどのような結果がもたらされるかについて、マウスの移植モデルを用いて検証した。すなわち、マウスのリンパ腫由来の細胞株である EG7-OVA をマウスの皮下に移植した後、免疫を促進させる薬剤（アジュバント）を投与してがん免疫を誘導することにより、腫瘍の縮小効果を通常の EG7-OVA 細胞と 3′非翻訳領域を欠失させた EG7-OVA 細胞について検討した。通常の EG7-OVA 細胞を移植したマウスでは、がん組織内に T リンパ球（CTL）が浸潤し、がん細胞は縮小した。しかし、上記の方法
で 3′非翻訳領域を欠失させて PD-L1 遺伝子の発現上昇を誘導した EG7-OVA がん細胞では、CTL の浸潤やがん細胞の縮小はほとんど認められない（図 1.12）。
19図 1.12: Pd-L1 ゲノム異常によりがん細胞の増殖と免疫回避が促進される
この結果は、PD-L1 ゲノム異常を持つがん細胞は、PD-L1 遺伝子の発現上昇を介して、CTL
によって担われるがん免疫を回避できることを示唆している。さらに、3′非翻訳領域を欠失させた EG7-OVA 細胞を移植したマウスにアジュバントと併用して抗 PD-L1 抗体を投与すると、通常の EG7-OVA 細胞を導入した場合と同様に、CTL の浸潤やがん細胞の縮小が認められた（図 1.13）。これらの結果は、免疫チェックポイント阻害剤を用いることにより、PD-L1 遺伝
子の 3′非翻訳領域の異常を有するがんの治療が可能であることを示すと考えられる。

20

図 1.13: 抗 Pd-L1 抗体の投与により Pd-L1 ゲノム異常による腫瘍増殖・免疫回避は阻害される
iii. まとめ
今回の研究によって、がん細胞が免疫を回避する新たなメカニズムが解明されるとともに、

PD-L1 遺伝子の発現調節において、3′非翻訳領域が重要な役割を担っていることが明らかになった。ほとんどのがん種において、本ゲノム異常が存在する症例は限られているが、がん全体
でみれば、最も普遍的かつ高頻度に認められるゲノムの構造異常の一つであると考えられる。
さらに、本研究で示された最も重要なことは、PD-L1 遺伝子の 3′非翻訳領域の異常は、これを有するがん細胞が積極的にがん免疫を回避していることを示しており、従って、抗 PD-1/PD-L1抗体を用いた免疫チェックポイント阻害に対する顕著な臨床効果が期待される症例を同定するための有用な遺伝子マーカーになり得るという可能性である。また、現在抗 PD-1/PD-L1 抗体の効果予測に用いられている免疫染色において、PD-L1 の後半部分（C 末端）を認識する抗体が使用されていることの問題点も指摘している。特に、PD-L1 の遺伝子異常が高頻度に認められる一部のがん種（ホジキンリンパ腫など）において免疫チェックポイント阻害剤が極めて有効で、しばしば治癒も期待できることに鑑みると、この可能性は早急に検証されるべきであると考える。
iv. 今後の展望と期待される波及効果
今回の研究成果を踏まえた、今後の最も重要な課題は、本異常を有するがんに対して、抗PD-1/PD-L1 抗体を用いた免疫チェックポイント阻害が、期待されるような顕著な臨床効果を示
すかどうかを臨床試験によって検証することである。実際に、顕著な臨床効果が確認された場
合には、本ゲノム異常がヒトのがんで最も普遍的かつ高頻度に認められるゲノムの構造異常の一つであることを考えると、その波及効果は極めて大きい。また、本異常をバイオマーカーとして用いた分子標的治療の実現する可能性がある。実際、本異常のバイオマーカーとしての意義を検証するために、「再発または難治性の成人 T 細胞白血病・リンパ腫に対するニボルマブの第Ⅱ相医師主導治験【UMIN000020601】」が現在鹿児島大学を中心として進行中である。
21

(b) 骨髄異形成症候群におけるクローン進化の解明－がんの時間的多様性・起源の理解
i. 背景
急性白血病の原因は、血液のもとになる骨髄細胞におきるゲノム異常と考えられているが、具体的にどの異常が原因になっているのか、またそれらが起こるタイミングなどはこれまで不明であった。従来研究では、1）解析症例数が少ない、2）解析遺伝子の数が不足、3）時系列で複数回解析される例が少ない、などが原因で、十分な成果が得られていなかった。今回、急性白血病の前がん病変である骨髄異形成症候群と、そこから発生する二次性の急性白血病を、1）
これまでになく多数例において、2）網羅的な遺伝子解析を、3）時系列で複数回行うことによって、急性白血病の原因となる遺伝子異常を明らかにした。
ii. 研究手法・成果
骨髄異形成症候群および二次性急性白血病 699 例において、ゲノム上の全ての遺伝子領域（およそ 2 万遺伝子）あるいは 60 あまりの重要遺伝子を対象に次世代シーケンサにより異常を発見、
し、そのうち、122 例では時系列で複数の検体を採取し解析した。さらに既に報告されたデータ
ベースより変異情報を入手し、合計2,250 例を解析した。

図 1.14: 骨髄異形成症候群から急性白血病振興への 2 ステップ
① 一症例で平均 10 個の遺伝子変異が認められ、変異の数・異常クローンのサイズ・変異遺伝子の種類は、病期が進行するにつれ、経時的に増加する傾向がある。
② 病期進行に伴いどの遺伝子変異が起こるのかを調べたところ、高リスク骨髄異形成症候群症例から二次性白血病に進展する際に高頻度に認められるタイプ 1 異常（FLT3, PTPN11,

WT1, IDH1, NPM1, IDH2, NRAS 遺伝子変異）と、低リスクから高リスクに骨髄異形成症候群が進行する際に起こるタイプ 2 異常（TP53, GATA2, KRAS, RUNX1, STAG2,

ASXL1, ZRSR2, TET2 遺伝子変異）が判明した（図 1.14）。また、この解析により、非常に興味深いことに、逆に白血病進展を起こしにくい遺伝子変異も見い出された（SF3B1 遺伝子）。
③ 時系列での解析の結果、タイプ 1 異常は経過中に新たに獲得されるか、もともとある異常クローンのサイズが増加する傾向が著明であることも解った。
④ 生存期間の解析によると、タイプ 1 異常をもつ症例は最も短期間で白血病進展を起こし、早22

期死亡することが判明した。一方、タイプ 2 異常は、タイプ 1 異常に比較すると白血病進展までの期間は長めだが、タイプ 1 あるいはタイプ 2 異常のない症例に比較すると白血病進展は有意に短期間で起こる ことも解った（図 1.15）。一方、SF3B1 変異陽性症例は、タイプ 1、タイプ 2 異常を持たない場合、ほとんど白血病を発症しない。

図 1.15: タイプ 1 とタイプ 2 の遺伝子異常の臨床的意義
今回の研究により、世界で初めて、2,000 例以上の骨髄異形成症候群および二次性急性白血病
において網羅的遺伝子解析を行い、白血病発症には遺伝子変異の蓄積が必須であることを明白に証明し、白血病進展に重要な2 つのタイプの遺伝子群を明らかにした。
iii. 波及効果、今後の予定
これらの遺伝子異常は、骨髄異形成症候群の低リスクから高リスク症例への進行、および二次性急性白血病進展を予測するマーカーとしてスクリーニングに用いられることが期待される。具体的には、タイプ 2 の遺伝子をシーケンスし骨髄異形成症候群の高リスクへの進行を予測し、タイプ 1 の遺伝子を確認し急性白血病への移行を予測する。
白血病進展を予測するマーカーが明らかとなれば、臨床応用により早期に白血病を発見し治療を開始し、治療効果を高められる可能性がある。
そのためには、今回の研究の結果を基に臨床研究を推進する必要がある。
(2) がんの起源と集団内多様性を理解するためのシステムの開発
A) 早期大腸がんの腫瘍内不均一性の解析－がんの空間的多様性の理解
がんの分岐進化により生み出される腫瘍内不均一性はがんの治療困難性の原因の一つと考えられ注目を集めている。次世代シーケンサの登場により、一つの腫瘍の複数の領域から DNA サンプルを取得し、シーケンスする多領域分割シーケンスにより、様々ながん種での腫瘍内不均一性の存在が明ら

23

かになってきた。いくつかのがん種については既知のドライバー変異がサブクローナルに存在し、腫瘍内不均一性の少なくとも一部はダーウィン進化により形成されていることが明らかになっていた。
我々はこれまでに、業務協力者の九州大学別府病院三森功士教授と進行大腸がん 9 例の外科手術サンプルの多領域分割シーケンスを行い、広範な腫瘍内不均一性を見出す一方で、大腸がんにおいては著明なサブクローナルなドライバー変異の蓄積は見られなかった。しかしながら「京」を用いたがんの進化シミュレーションによって突然変異率が十分高い時に中立進化によって実験データと同様の広範な腫瘍内不均一性が再現できることを見出し、進行大腸がんにおいては中立進化が主な腫瘍内不均一性の生成原理であるという結論に至った。
その一方で、早期大腸がんの腫瘍内不均一性についてはほとんど研究がなされていなかった。そこ
で H28 年度は早期大腸がん 10 例の内視鏡手術サンプルの多領域分割シーケンスを行い、進行大腸がんデータとの比較解析を行った。進行、早期がんそれぞれの多領域変異プロファイルから進化系統樹を作成しクラスタリング解析を行ったところ、進行大腸がんと比べて早期大腸がんにおいては幹が短く枝が短い傾向があり（すなわちクローナル変異の割合が少なく、サブクローナル変異の割合が多い）
異なる形状を示すことが明らかになった（図 1.16）。また、既知のドライバー変異の分布を評価した
ところ、進行大腸がんに比べて、ドライバー変異が有意にサブクローナル変異として存在する傾向がある、すなわち早期大腸がんの腫瘍内不均一性は進行大腸がんに比べてダーウィン進化の寄与が大きいことが明らかになった（図 1.17）。これにより固形がんの時間的・空間的多様性の理解が進んだ。

図 1.16: 進化系統樹のクラスタリング24

図 1.17: ドライバー変異の分布
以上の結果より、大腸がん進化の初期では腫瘍内不均一性の生成原理は主にダーウィン進化により生成され、がん化が進行すると中立進化に切り替わるという新たな大腸がんの進化モデル、ダーウィ
ン中立遷移モデル（Darwinian-neutral transition model）を構築した（図 1.18）。これにより、来年度は腫瘍内不均一性の生成原理の切り替えがどのような機構で起こっているかを、がんの進化シミュ
レーションにより明らかにするという方向が見えてきた。

図 1.18: ダーウィン中立遷移モデル

25B) 超並列シミュレーションによるがんの進化原理探索に向けたフレームワークの整備
これまでに我々はがんの進化シミュレーションモデル BEP モデルを構築している(Uchi et al.
PLoS Genet. 12(2):e1005778, 2016)。BEP モデルはセルオートマトンモデルの一種で、パラメータセ
ットに依存して各細胞がゲノムに異なる組み合わせ変異を蓄積していきながら、細胞集団が不均一性を生成しつつ成長していく様子を再現できる。H28 年度はこの BEP モデル用いた超並列シミュレー
ションによりがんの進化原理を探索するためのフレームワークを整備した。具体的には 1)シミュレー
ションモデル構築・改良、2)様々なパラメータセットで超並列シミュレーション、3)各シミュレーシ
ョンより得られた変異プロファイルの要約統計量を計算、4)インタラクティブな可視化ツールを用いてパラメータ依存性を解析、の 4 ステップを東大医科研ヒトゲノム解析センターのスーパーコンピュ
ータ上でシームレスに行えるパイプラインを構築した。特に 4)の可視化のステップについては 4 次元のパラメータ空間全体をウェブブラウザ上でインタラクティブに探索できる、並列シミュレーション
結果の可視化手法一般としても新規性の高いツール BEP エクスプローラーの開発に成功した。以上の H28 年度に開発したリソースを用いて来年度以降がんの進化原理の探索を行っていくための基盤ができた。

図 1.19: BEP エクスプローラーのスナップショット(3) 疾患に関連する胚細胞バリアントの効率的な解析システムの開発
胚細胞バリアントを検出するためのデータ解析手法を検討した。複数非適合対照検体（患者検体の対照ではない検体）を利用した低頻度体細胞変異（細胞集団中に少ない割合の細胞が有している体細胞変異）をできるだけ高感度に検出する手法を検討することにより、胚細胞バリアントを高感度に検出することを考えた。一般に、アレル比率が低くなるにつれて、「真の変異」と、「シーケンスエラーなどの様々なエラーやアーチファクト」との区別が難しくなるために、既存の変異検出法ではアレル比率 1〜5％以下の変異の検出は難しいと考えられている。そこで、H28 年度は、低頻度体細胞変異を高感度に検出す
る方法について研究開発を行った。
まず、H28 年度の研究の基礎となっている技術について述べる。これらは、サブ課題 A の業務参加者

26

らが、本重点課題が始まる前に開発したものである。非適合対照検体の利用の利用に関して、偽陽性を除くために、これまでの VarScan2 [1] などで用いられてきた Fisher 正確検定の P 値に加えて、非適合対照検体の変異アレル比率について以下の知見を得ていた（図 1.20）。
•

真の変異については、非適合対照検体のほぼ全てにおいて、その場所で変異アレルが観測されない。

•

シーケンスエラーについては、非適合対照検体のいくつかで、その場所において変異アレルが複数観測される。
これにより、シーケンスエラーが起こりやすい領域があり、また、シーケンスエラーの起こりやすさ

はある程度検体間で共通しているということがわかる。

図 1.20：左の変異候補については、非対照検体の同一箇所におけるミスマッチ率のうち多くは 0 になった。それに対して、右の変異候補については、複数の検体で高いミスマッチ率が観測された。この結果
から、右の候補箇所は、シーケンスエラーが生じやすい箇所において、偶然に腫瘍検体に特異的にシーケンスエラーが生じた偽陽性と考えられる。
この観察結果に基づき、複数の非適合対照検体を利用して体細胞変異の検出精度を向上することが可能となる。例えば、候補となる体細胞変異から偽陽性を除くための簡易的なフィルターとして、次のよ
うな手続きが考えられる。
1. 体細胞変異の候補について、複数の非適合対照検体（例えば 10 検体）における、変異アレルの比率を調べる。
2. 事前に決めた条件を上回っていたらその候補を除去する（実際のルールとしては例えば、「変異アレ
ル比率が 2％以上の検体が 3 検体以上観測される」など）。
ここで、何検体の非適合対照検体を用いれば良いのか、また、具体的なフィルタリングルールをどのように設定するか、などの問題が起こってくる。これらを解決するために、経験ベイズ理論に基づいて、理論的に妥当な基準で、複数の非適合対照検体を利用してフィルターする方法論を創り、アプリケーシ
ョンとして実現したものが EBCall [2]である。これにより、通常のエクソーム解析において効率的に偽陽性を除去しつつ、感度の高い変異検出が可能になることを実証した。世界最先端の成果が数十の論文

27

として発表されている。そして、この EBCall による偽陽性を除去は、低頻度体細胞変異にも有効である。
実際、2%～10%の低アレル頻度の変異についても、80%～90%の正答率が得られ、さらに真の変異と確認された変異の変異アレル比率のヒストグラムから主要なクローンの集団の他に、低アレル頻度のサブクローンの集団を観察することができた。
更なる高精度化のためのアイディアは SNP との関連を見ることである。体細胞変異の候補となる場所の周辺を、IGV（http://software.broadinstitute.org/software/igv/）などのブラウザーで観察していると、しばしば、体細胞変異の候補となる場所の周辺に、ヘテロ SNP が観察される。さらに詳細な観察をすると、「体細胞変異の変異アレルと、ヘテロ SNP のアレルの共起関係」について以下のような興味深い事
実が観測された（図 1.21）。
•

真の変異の場合には、体細胞変異の変異アレルは周辺のヘテロ SNP のどちらか一方アレルとのみ共起する。これは、体細胞変異が通常、片方のハプロタイプにのみ生じるためである。

•

シーケンスエラーの場合は、体細胞変異の変異アレルはヘテロ SNP の両方のアレルと共起する。これは、シーケンスエラーは両方のハプロタイプにランダムに生じることを示唆する。

図 1.21：真の体細胞変異、シーケンスエラーにおいて、それぞれ近隣のヘテロ SNP との共起関係の概念
図。左の図においては N_a=4、N_A=0 となり、N_a/(N_A+N_a)=1 となる。右の図においては N_a=2、N_A=2となり、N_a/(N_A+N_a)=0.5 となり、0 または 1 から大きく隔たりがあるため、偽陽性と考えられる。
この観察結果に基づき、周辺のヘテロ SNP を利用することで、体細胞変異の検出精度を向上することができる。例えば、候補となる体細胞変異から偽陽性を除くための簡易的なフィルターとして、次のよ
うな手続きが考えられる。
1.

体細胞変異の候補について、周辺のヘテロ SNP（ここではメジャーアレルが A、マイナーアレルが aと表す）との共起関係を調べる。ここでヘテロ SNP は dbSNP などに登録のある信頼できるものに制限するのが良いであろう。

2.

ヘテロ SNP の 2 つのアレル（A、a)について、体細胞変異の候補と共起している（同一のシーケンス
リードに生じている）回数をカウントする（N_A、N_a とする）。

3.

共起しているアレルの比率(N_a/(N_A+N_a))が 0 または 1 から大きく隔たりがあれば（例えば 0.1 ~ 0.9の中に入っているなど）、該当する体細胞変異を信頼できない変異としてフィルターする。
さらに、こうした手続きを、ベイズ統計学により理論的に妥当な形で定式化した変異検出手法が28

HapMuC [3]である。通常のエクソーム解析において周辺の SNP を考慮することにより性能が向上することを観察している。このテクニックは、候補の体細胞変異とヘテロ SNP が十分に近くにあり、複数本のシーケンスリードがこれらをカバーできる必要がある。そのために、特に 100bp ほどのリード長では適用範囲が限られてしまうという弱さがあるが、今後リード長が長くなることで、こうしたテクニックの有用性が増してくると考える。
以上が、H28 年度の開発の基盤となっている研究である。本サブ課題で用いている次世代シーケンサ
は Illumina 社製のもので、多くの場合リードをペアとしてシーケンスする。その際に、インサートサイズ
が小さい場合、二重に読まれる箇所が生じる。実はこうしたオーバーラップの問題は、シーケンス解析のコミュニティでは長らく問題になっていた。一つの主要な問題点は、変異アレルのダブルカウントにより、実際よりも統計量が高くなってしまい、偽陽性に繋がってしまうことである。このために、二重にシーケンスされた箇所については、BamUtils (http://genome.sph.umich.edu/wiki/BamUtil)などのソフトウ
ェアを用いて、重複して読まれた領域のうち一本のリードをマスクすることで、ダブルカウントを防ぐというようなことがしばしばなされてきた。また、最新の samtools (version 1.3.0: http://www.htslib.org)では、mpileup コマンドにより塩基の集計を行う際に、overlap している箇所については、一方の base quality
を 0 にして、二重カウントを防ぐような仕様になっている。
本年度の研究として、我々は、このオーバーラップの現象に着目して、真の体細胞変異とシーケンス
エラーにより、以下のような現象が見られることを発見した（図 1.22）。
•

真の変異の場合は、オーバーラップしている場所に置いて、両方のリードペアのコールが一致したものになる。

•

シーケンスエラーの場合には、しばしばオーバーラップしている場所のコールに不一致性が生じる。
すなわち、シーケンスエラーは、両方のリードペアである程度の独立性をもって生じている。

図 1.22：対になっているリードはリードペアを表す。真の体細胞変異、シーケンスエラーにおいて、それぞれオーバーラップ箇所における一致性の概念図。左の図においては、RR=2、RA=0、AA=3 となり、オーバーラップにおける不一致の割合 RA/(RA+AA)が 0 となる。一方で、右の図においては、RR=1、RA=3、
AA=1 となり、RA/(RA+AA)=0.75 となり、偽陽性と考えられる。
この観察結果に基づき、このオーバーラップの情報を利用して、体細胞変異の検出精度向上することができる。例えば、候補となる体細胞変異から偽陽性を除くための簡易的なフィルターとして、次のよ29

うな手続きが考えられる。
1.

体細胞変異の候補について、その場所を両方のリードペアが覆っている各々のリードペアについて、それらのベースコールを以下の 3 つに分類して集計する

2.

(RR) 両方ともリファレンスアレル

(RA) 一方がリファレンスアレル、もう一方が変異アレル

(AA) 両方が変異アレル

少なくとも一方が変異アレルを有するリードペアのうち、オーバーラップにおいてベースコールが不一致の割合 (RA/(RA+AA)) を調べ、これがある閾値以上であれば（例えば 0.2 など）
、候補を除
去する。
さらに、我々は、こうした方法論をベイズ統計学により定式化し、新たな変異検出手法として OVarCall[4] (Moriyama et al. 2016: https://github.com/takumorizo/OVarCall）を開発した。実際、この手法によりオー
バーラップを考慮すると、エキソームシーケンスデータからの変異検出の性能が向上することをシミュ
レーションデータ、および実データを用いて確認した。
このテクニックを低頻度体細胞変異検出（〜1%）に応用する際には、いくつもの問題点が残っている。
第一に、シーケンス直前のライブラリの段階で既に生じているエラー（元々の生体内ではなかったが、サンプル調整の過程で生じたと考えられるものなど）については、オーバーラップにおいて、両方のリードペアで一致した変異が観測されるはずであり、このテクニックでは除去できない。また、このテクニックが利用できるのは、インサートサイズが実際に読まれるリード長と比較して短い場合に限定されてしまうという問題もある。そのために、オーバーラップにおける変異の一致性については、「これが満たされていないならば偽陽性と言えるが、満たされているからといって真の変異とは言えない」、という
必要条件的なものになるであろう。
[1] Koboldt et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome
sequencing. Genome Research. 22(3): 568-576 (2012).
[2] Shiraishi Y, Sato Y, Chiba K, Okuno Y, Nagata Y, Yoshida K, Shiba N, Hayashi Y, Kume H, Homma Y, Sanada
M, Ogawa S, Miyano S. An empirical Bayesian framework for somatic mutation detection from cancer
genome sequencing data. Nucleic Acids Research. 41(7): e89 (2013).
[3] Usuyama N, Shiraishi Y, Sato Y, Kume H, Homma Y, Ogawa S, Miyano S, Imoto S. HapMuC: somatic
mutation calling using heterozygous germ line variants near candidate mutations. Bioinformatics 30(23):3302-3309 (2014).
[4] Moriyama T, Shiraishi Y, Chiba K, Yamaguchi R, Imoto S, Miyano S. OVarCall: Bayesian mutation calling
method utilizing overlapping paired-end reads. IEEE Trans Nanobioscience. 2017 Mar 1. doi:
10.1109/TNB.2017.2670601.

30

(4) 用語集
DNA の場合、DNA を構成するヌクレオチド(A, T, C, G)の結合
順序（塩基配列）を決定することをシーケンスするという。
シーケンス
4.

本研究で用いているイルミナ社の次世代シーケンサでは、100文字程度の長さに断片化した DNA 試料の配列を読み、コンピ

4-1.

ュータに出力する。RNA も同様にシーケンスできる。

(1)

DNA や RNA のシーケンスデータからマッピングやデータ解
Genomon析などを行い、遺伝子変異や融合遺伝子を見つけ出すプログ
ラム群の総称。
DNA チップなどで得られる遺伝子発現データやタンパク発現

SiGN (gene network

データ、miRNA を含む noncodingRNA 発現データなどから細estimation software )胞内分子の発現制御システムのモデル化・予測を行うスーパ
ーコンピュータ用ソフトウェア群。
ヒトゲノム中の全エクソン領域（直接的にタンパク質の構造エクソームを決定している部分）の集合体のこと。エクソンを分断している領域はイントロンと呼ぶ。
生殖細胞（卵や精子）における遺伝子変異（バリアント）で、胚細胞バリアント子孫の全ての細胞に組み込まれるもの。生殖細胞内のバリア
ントは親から子孫へと受け継がれる。生れてから何らかの原因でゲノムに入る変異の対語。

4.

4-1.

データ同化

(2)
4.

4-1.

サルコメア

(3)
4.
1.

4-2.
コデザイン

実測値を活用して、シミュレーションモデルから現実に近い解が得られるようモデルを改良すること。
筋肉が収縮するための最小構造単位。2 種類のタンパク質から成る。
ハードウエアとソフトウェアの設計を、開発の初期段階から協調して行うこと。開発期間短縮やシステム最適化を狙える。

(1)
グリッドエンジン
バーチャル・グリッドエンジン(VGE)

ネットワークバイオインフォマティク
ス

1 台のマシン上で実行していたプログラムを、複数の計算機上で円滑に逐次実行することを可能にする管理ソフト。
グリッドエンジンの搭載されていないスーパーコンピュータ
上で、MPI を利用してグリッドエンジン機能を提供するために開発されたミドルウェア。
遺伝子の発現調節が相互に影響しあって調節されている様子
を、遺伝子の相互作用として表した図、あるいはその手法。
生命について研究する際に、生命が持っている「情報」と言えるもの、例えば遺伝子やタンパク質の構造などを分析する
分野

31

MPI (Message-Passing
Interface)

MPI Forum が制定している、並列計算器における複数プロセス間での通信規格の名称。MPICH や OpenMPI がその代表的実装例。
インタープリター言語（1 行ごとに逐次的に処理する言語）のpython実装

一つ。
あらかじめ機能（/関数）だけが決められている箇所を、具体的なプログラムとして作成している部分や、その作業
OS とアプリの間に入るシステムのようなもののうち、データ
ベース管理システムやトランザクションモニターのような、ミドルウェアオペレーティングシステムの機能の拡張、あるいはアプリケーションソフトウェアの汎用的（共通的）な機能を集めたものを指す
ソフトウェアを構成する単位であって、ある一連の独立したモジュール機能を実現するプログラム部分。
その部分だけで独立してコンパイル可能な単位であり、かつ再利用可能な単位。

ライブラリ

汎用性の高い複数のプログラムを再利用可能な形でひとまとまりにしたものだが、実行ファイルでない場合もある。
例えば日本語入力システムやウイルス対策ソフトのように、常駐システム起動時にプログラムをメモリーに読み込み、そのま
まメモリー上に保持すること。

ジョブ

コンピュータにさせる仕事の単位。一連のプログラム群
送るデータをキューと呼ばれるデータ領域に保持し、データ

キューイングを受ける側の処理が完了するのを待たずに次の処理へ移る方式

ロード

記憶媒体にセーブされたデータ（セーブデータ）から作業内
容や各種データを復元すること。
変数やサブルーチンの有効範囲である名前空間。

パッケージ

モジュールと同様に、変数名やサブルーチン名の重複を防ぐために用意されている機構。
プログラミングやソフトウェア開発の分野で、ソフトウェア

ラッパー

やプログラム部品などが提供するクラスや関数、データ型などを本来とは異なる環境や方法で利用できるようにしたもの
ある操作の結果が次の操作の入力になるようにハードやソフ

パイプライン

トを並べておく技術。操作が並行して行われているように見せることができる。

遺伝子発現

遺伝子の働き、つまり RNA（またはタンパク）が合成される

32

こと
バッチジョブソケット通信処理手順をファイルに記述しておき、一定量のデータを一括処理する方法。実行に際し臨機応変なユーザ対話は行わない。
TCP/IP プロトコルを用いたコンピュータ間の電話概念的通信。IP アドレスと ポート番号が電話番号に相当プログラムがコンパイルされ、機械語に翻訳された状態をオ
プジェクトという。全てのオブジェクトとシステムライブラ

共有オブジェクト

リが結合（リンク）されて初めて実行可能状態になる。オブ
ジェクトの中には、複数プログラムから利用されるものがあり、そのようなものが共有オブジェクト（shared object）として特別な処理を受ける。

SV（Structural Variation:

染色体中に数百塩基から数百万塩基にわたる大規模な変化を

構造変異）

伴う変異。

BWA（Burrows-Wheeler
Aligner）

DNA のショートリード配列の参照ゲノムに対する位置同定（マッピング）を行うソフトウェア。現在のデファクトスタ
ンダードソフトウェアである。

SIMD (single instruction

複数のデータ列に、1 つの命令を施す操作。コンピュータ内でmultiple data)

並列処理されている形態。

I/O 負荷頻繁にデータを出し入れすることにより、入出力関連のハード ウェアやネットワークにかかる負荷のこと。
複数ユーザによる共同利用環境において、ユーザごとのプロ

スケジューラグラム実施要求（ジョブ）を割り当て可能な資源内に置いて効率的に順位付けし、実行制御する管理機構のこと。

サブミット

プログラムの実行要求をスケジューラに渡すこと。

ショートリードデータ

Wet な実験と一次配列解析による、ショートリードの配列と（FASTQ 形式）

各塩基のクオリティを示したファイルの一種
ヒトゲノムは核ゲノムとミトコンドリアゲノムから成る。ヒ
トゲノム計画は 2003 年に解読終了が宣言された。ただし複数

標準ゲノム配列

国複数人の情報から多数決的に決められたのであって、例えば日本人によく見られる構造多型の情報は、この国際参照配列には含まれないのではないかとも考えられている。

プロファイリングプログラムの高速化等のために、各処理における実行時間、消費資源量等の分析データを取得すること。
マルチメディア拡張命令セット。複数の浮動小数点数演算をSSE (Streaming SIMD
Extensions) 命令
Haswell

同時に実行する機能や、マルチメディア処理に頻繁に使う命令が追加され、3D グラフィックスの描画などに適する。
インテルによる第 4 世代 Core アーキテクチャとなる CPU。ノ33

ート PC（モバイル PC）市場がターゲット。
人間がある言語で書いたプログラムをコンピュータが理解でコンパイラ

きる機械語に一括翻訳すること。その後、プログラムを実行することになる。
確率変数の定性的因果関係を矢印で結んだ図であり、変数間の定量的依存関係は条件付き確率で表現してある。本件ではベイジアンネットワーク

遺伝子とその発現量が確率変数とその値に相当する。静的モ
デルでは依存関係はサイクルを持たないが動的モデルでは可能である。
遺伝子間の制御の因果関係を推定するために、遺伝子発現デ

ベイジアンネットワーク

ータなどを利用し、非線形回帰モデルなどと組み合わせベイ

推定

ジアンネットワークの構造を推定する手法。遺伝子ネットワ
ーク推定に用いられる。

L1 正則化法
非線形回帰

離散スコア

過学習を防ぐために加える正則化項の指数が 1 である回帰法。
出力が入力の一次結合で表現できない場合を扱う統計学上の解析手法。
度数表で表されるような離散データにさらに重みをつけて評価したもの。
タンパク質をコードしない、つまりタンパク質へと翻訳されノンコーディング RNA

る伝令 RNA（messenger RNA, mRNA）を除く全ての RNA の総称であって、一部は最近発生、疾患との関わりが明らかになりつつある。

マイクロアレイ

検査・実験対象物（例えば DNA やタンパク質）を多数固定して、相同性や反応を、検出したり定量したりする方法。
細胞の中の mRNA や miRNA の配列を解読して、発現量の定RNA-seq

量、新規転写配列を解析する手法。細胞内の RNA の絶対量を数えることができるようになる可能性がある。
シーケンサの出力するデータは DNA/RNA の断片情報である。

マッピング

この断片（リードという）を何らかの参照配列と比較することで、そのリードの位置を決定する作業をマッピングという。

PD-1

腫瘍を攻撃する活性化 T 細胞の表面に現れる受容体タンパ
ク質（別名：CD279）
抗原を提示している細胞、つまり腫瘍の表面上に発現していPD-L1るタンパク質。PD-L1 が PD-1 に結合すると活性化 T 細胞の活動が抑制される。

免疫チェックポイント

免疫過剰な攻撃を続けたり、正常な細胞を攻撃したりしない
ように抑制する体のチェック機構。

34

3′非翻訳領域
転写産物
EG7-OVA 細胞

メッセンジャーRNA においてタンパク質コーディング領域よりも下流にあって、タンパク質を生成させない領域
RNA ポリメラーゼによりゲノム DNA を鋳型として、DNA に相対する塩基で合成された全ての RNA。
卵白アルブミン（OVA）、マウスリンパ腫細胞（EG7）
例えば、がんに特徴的な物質が産生されるときに、そのようマーカー

な物質のうち体液中（主として血液中）で測定可能なものが、いわゆる「腫瘍マーカー」として臨床検査の場で使われる。
特定の障害の保有または発生リスク、あるいは特定の障害に関する遺伝子の保有リスクが高い個人を同定するために計画される遺伝学的検査。 ゲノム上にランダムに遺伝的な変異をスクリーニング

導入させた後、それらの処理によってランダムに変異を持った生物集団に対して特定の表現型を示す個体を単離し、その個体がどのような遺伝子に変異を持っているかを同定する。
それによりその表現型にかかわる遺伝子と特定できる。

4.
1.

分類対象の集合を、内的結合 (internal cohesion) と外的分離

4-2.
クラスタリング解析(external isolation) が達成されるような部分集合に分割すること。

(2)
クローナル変異既存の有機体の一部から作られた同一の複製

BEP (Branching

細胞集団が不均一性を生成しつつ成長していく過程を解明す

Evolutionary Process)るためのモデル。

35

2. データ同化生体シミュレーションによる個別化医療支援（サブ課題 B）(1) 全脳循環シミュレータの開発（阪大 和田Ｇ）
A) ボクセル型流体解析の基礎的手法構築
ポスト「京」の計算性能を活かした超大規模並列計算を行うため、直交格子を用いた解析アプロー
チを構築した。この際、適用した従来手法に内在する問題点を検討し、これを克服する新たな定式化を提案した。また、解析で必要となるプリ・ポスト処理に関して機能の充実化を図った。
(a) 方法論の基礎的検討
直交格子を用いたシミュレーションは領域分割による並列化と相性が良いが、血管の任意形状を表現するための手法が必要となる。Weymouth(J Comput Phys, 230, 2011, 6233-6247)により提案されている Boundary data immersion(BDI)法では、媒質相（流体相、固体相など）を記述する滑らかな指標関数を導入することで直交格子において任意形状を取り扱っており、数値アルゴリズムの観点から非常に優れたアプローチとなっている。しかしながら Weymouth の BDI 法では、非定常非圧縮性 Navier-Stokes 方程式への適用の際、時間刻みが定常解に含まれるような定式化がなされているため、定常解を議論する際にその取扱いを注意する必要がある。そこで本研究では、これらの問題点を回避する新たな定式化を提案した。
流体相および固体相での支配方程式は以下のように与えられる。

ì Ñ× v = 0
in Ω fïïí r (¶t v + v × Ñv) = -Ñp + hÑ2 v in Ω fïin Ω sïî v = 0(1.1)

ここで、流体領域 Wf では非定常非圧縮性 Navier-Stokes 方程式が、固体領域 Ws では静止剛体条
件が与えられている。BDI 法では、有限幅の相間の変化領域 G を導入し、以下ように与えられる滑らかな相関数を定義する。
ì 1
in Ω f \ Gïï
c (x) = í Î (0,1) on G
ï0in Ω s \ G
ï
î

(1.2)

その後、相関数を用い各相の支配方程式を混合化することで、解析領域 Wf =Wf Ws での混合化方
程式を得る。解法にフラクショナルステップ法(Dukowicz & Dvinsky, J Comput Phys, 102, 1992,
336-347)を適用すると、Prediction 相はì æ *öv - vnhï+ Adv(v n ,v n-1 ) ÷ = -Ñp n + Ñ 2 (v n + v* ) in Ω f
ï rç2í è Dtøï *in Ω sïî v =0(1.3)と与えられる。ここで、上添え字は時間ステップを意味する。オリジナルの BDI 法では、Prediction
相における流体の支配方程式の右辺を時間積分した形で与え混合化を行っているため、定常状態となる時間無限大の極限もしくは時間刻み dt 無限大の極限においても時間刻み dt の影響が入ってしまう。そこで我々は右辺の時間離散化を踏まえた上で定式化を行った。この結果、Prediction 相に36

おける混合化方程式は、
2ææ r h 2 ö -1 *r h ö÷ n
n
n-1
n
ç
Ñ
(
c
v
)
=
r
Adv(v
,
v
)
Ñp
+
+
v
ç
÷
ç Dt 2 ÷è Dt 2 ø
è
ø

in Ω

(1.4)

と与えられ、時間刻み dt 無限大の極限では、上式は（圧力項を除いて）定常支配方程式の離散化表
記に帰着する。
(b) プリ・ポスト処理の充実
理化学研究所が提供する符号付距離関数作成ライブラリ(V-SDFlib, VCAD, RIKEN)を利用し、画像処理により抽出された血管表面ポリゴンをボクセル格子上の体積率関数に変換した。また、領域
分割並列によりバラバラに出力されたデータに対し、それらを結合するシステムを構築した。さらに二次関数を用いた移動最小二乗法を利用することで、ボクセル格子での体積率等値面における壁面せん断応力の評価手法を新たに確立した。以下に例として、ボクセル格子において構築された体
積率関数の等値面(=0.5)、本手法により評価された壁面せん断応力を示す。

図 2.1: ボクセル格子と構築された体積率関数

図 2.2: ボクセル格子での体積率等値面(=0.5)におの等値面（=0.5）

ける壁面せん断応力(WSS)の評価(c) アプリ開発状況
本アプリ(Voxel based Brain Blood simulation -Flow: VBB-Flow)は他のアプリの骨格となるため、単体性能、並列化性能などの面を常に改良・更新しながら開発を行っている。並列化技法は直交領域分割による MPI＋OpenMP のハイブリッド並列である。本離散化においては、近隣のステンシルが必要となる。そこで、各担当ノードに自身の領域に加え境界領域を設定し変数を保持し、必要に応じて隣接ノード間で各変数を非同期通信した。また、ポアソン方程式の求解には本モデルでは対角スケーリング付き Bi-CGSTAB 法を用いた。この際、各担当ノード内でベクトル・行列積演算を行い、必要に応じて隣接ノード間で非同期通信、また全領域での和が必要とされるスカラー値を集団通信により処理した。現在、864 プロセス(MPI)×16 コア(OpenMP)の並列解析まで達成しており、今後、単体性能および並列性能の調査をしていく予定である。また、ボトルネックとなるポアソン
方程式の求解に幾何的マルチグリッド法の導入を検討している段階であり、次年度以降本格的に実装する予定である。

37B) 血球流動を伴う脳微小循環解析に向けた取り組み
脳循環機能は、微小循環スケールにおけるニューロン、グリア細胞などの神経系と血管系の連成による血行動態によって変化するとされており、脳循環機能の解明に向けた微小循環における血行動態のメカニズムの解明は重要な課題の一つである。微小循環においては、血液中に含まれる固体成分の一つである赤血球と血管のサイズが同程度になるため、三次元的な変形を伴う血球流動を考慮しつつ、数マイクロメートルから数十マイクロメートル径の複雑に分岐する血管網を取り扱う必要がある。また、メカニズムの本質を明らかにするためには、実際の血管形態を再現および定量化し、適切な条件のもと解析を行い結果の検討に移る必要がある。
これを達成するため本年度は、(a)計測画像を反映した脳微小血管モデル構築手法の改良、(b)血球流動を考慮した解析アプローチの構築を行った。
(a) 計測画像を反映した脳微小血管モデル構築手法の改良
前年度では、血管新生のアナロジーを利用した血管成長モデルを提案し、二光子顕微鏡観察で取得した脳微小血管画像と併せることで、マウス大脳皮質内の微小血管形態を構築した。しかしながら同手法では、血管径を反映する仕組みがなかったため、得られた形態モデルは全て一様な血管径となっていた。本年度はこの問題点を克服する取り組みを行った。
二光子顕微鏡により得られたマウス大脳皮質内微小血管群のスライス画像（1024*1024ピクセル、
250スライス）に対し、画像処理ソフトAmira(Visage Imaging, Berlin, Germany)を用いて細線化を行い、中心点群と関連する内接球半径を抽出した。
血管新生のアナロジーを利用した血管成長モデルでは、血管新生因子濃度cが以下の反応拡散方程式により輸送される。

¶t c = DÑ2c - acf + S(X;R)(1.5)

左辺第一項は時間変化項、右辺第一項は拡散項、右辺第二項は血管相における血管新生因子の消費項、Sは低酸素細胞からの発生項を表す。発生項では、抽出された中心点群を低酸素細胞と見なし
血管新生因子が放出される。この際、発生するソース領域は内接球半径Rを持つ球とすることで、血管径をモデルに反映させた。その他、数理モデルの詳細などは前年度に記載した手法と同様である。
画像処理により抽出された中心点群と内接球半径および初期血管相、また、異なる角度から見た本手法により構築された微小血管形態を以下に示す。結果より、大脳皮質の深さ方向に延びる太い血管が見られ、またそこから横方向に枝分かれした微小血管によりそれらが架橋している様子が見てとれた。これらは、穿通動脈および(準)毛細血管を表わしており、動脈から静脈への循環系がモデ
ルにより再現された。
現在の問題点として、構築されたモデルが内接球による凹凸を引きずってしまっているため表面平滑化の仕組みを盛り込むこと、また画像輝度の低さに由来する末端血管の存在が挙げられる。今後はこれらの改良を行っていくと同時に、構築された形態モデルを血流解析に適用していく。

38

図 2.3: 画像処理による抽出され鳥瞰図

画像平面視点た中心点群と内接球半径 R のカラーマップおよび初期血管相の等値

図 2.4: 本手法により構築された大脳皮質内微小血管形態モデル

面(=0.5)(b) 血球流動を考慮した解析アプローチの構築
直交格子を用いた任意形状内の流れ解析を行うため、開発中のVBB-FlowにImmersed boundary法(Peskin, 1972)による流体と赤血球の力学連成機能を新たに追加した。ここで、分岐を多数含む血管内を長時間かつ長距離流れる赤血球を扱う際、膜応力評価の安定性の向上、および血管分岐部における数値的不安定を緩和する必要がある。そこで本研究ではそれぞれ以下の工夫を行った。膜応
力評価の安定性向上に関して、曲げ応力を評価する新たな離散化手法を提案することで、従来用いられていた平均化操作を排除し数値振動を抑制した。具体的には、血球膜を表わす三角形要素の頂点において曲率Kを定義し重み付き残差法によりその値を求めた。血管分岐部における数値的不安定は血球表面を構成するメッシュが血管に埋まることが関係している。そこで、血管形状を表現する
符号付距離関数をボクセル格子において定義し、血球表面格子を構成する頂点に、壁面からの距離に応じた反発力を加えることで、表面格子の血管形状への接触を抑制した。
本アプローチの妥当性を検証するため、血流のみかけ粘度に関して、先行研究(Pries et al、 1992)と比較した。円管流路内に血球を複数配置し周期境界条件のもと圧力駆動させることで準定常解を求めた。この際、血球が存在しない系での流量（理想的にはハーゲン・ポアズイユ流れ）と血球を含む系での定常状態における流量比を評価することで相対みかけ粘度mrelを算出した。以下に円管直径Dと相対みかけ粘度mrelの関係を示す。なお、血球体積率（ヘマトクリット）が30%で一定になる
ように系のサイズは適宜変更した。以下に流れが十分発達した状態での赤血球の様子およびシミュ
レーションから得られた相対見かけ粘度と対応する実験結果のフィッティングカーブ(Pries et al,1992)を示す。赤血球の変形挙動はそれぞれ過去に報告されているSingle-file、 Multiple-fileを再現した。また、径が細い血管においては、赤血球は中心軸に対して対称な形状をとっているが、径が大きくなるにつれて、横を向いたように分布することが分かった。さらに、Fahraeus-Lindqvist効
果として知られているよう、円管直径の変化に伴い相対見かけ粘度が変化し、先行研究と同様の傾向を示すことを確認した。これらより、相対見かけ粘度の血管径依存性を定性的に再現できたと考39

えられ、本アプローチの妥当性を示すことができた。

図 2.5: 異なる円管径における準図 2.6: 円管直径 D と相対みかけ粘度 mrel の関係。Empirical fit

定常状態での血球流動は実験値のフィッティングカーブを意味する (Pries et al,

の様子

1992)

次に多数の分岐を含む血管系に対し本アプローチを適用し解析を行った。その一例を以下に示す。
結果より、複雑に分岐する微小血管系においても、流体を介した血球同士の相互作用を捉えることができ、その時系列変化も捉えることが可能であることが分かった。本シミュレーションでは、3sec
ほどまで計算が安定して走ることを確認しており、結果に関する統計的な議論が十分可能である。

t=0.25 (sec)t=0.5 (sec)t=2 (sec)

図 2.7: 分岐血管を流れる赤血球シミュレーションの時系列結果(c) アプリ開発状況
開発中のボクセル型流体・固体解析プログラム（VBB-Flow）を拡張することで派生アプリ
VBB-MicroVasc および VBB-Capsule を開発した。領域分割による MPI 並列において、各領域に含まれる計測点を分割コアに割り当て通信を行った。それぞれ、最大で 12000 ノード(MPI)×8 コア

40

(OpenMP)、864 プロセス(MPI)×16 コア(OpenMP)のハイブリッド並列を達成している。
C) 全脳循環血流解析に向けた画像ベース数理モデル血管形態モデルの構築
脳循環は、内頚動脈および椎骨動脈からウィリス動脈輪、脳表層動脈、脳実質微小血管網と巡り、同程度のサイズの静脈系経路をたどり戻ってくる。脳血管構造は階層的なツリー構造を持つだけではなく、各階層においてネットワーク構造（サイクル構造）を持つことが知られており、これらが脳循
環全体の調整に関係していると考えられている。これらを明らかにするためには、全脳循環シミュレータによる本質的な理解が必要となる。脳血管モデルは主に医用画像（血管造影画像）を用いて構築されるが、医用計測の空間分解能に限界があるため、細動静脈以下の細い血管は抽出されない。そこ
で、数理モデルを用い医用画像より抽出された血管形状と組み合わせることで、全脳血管形態モデルの構築を試みた。なお、ここではポスト「京」を用いることで、数拍動の脳循環の血流動態を再現できる規模の計算モデルの構築を想定している。
(a) 時系列血管造影画像を用いた画像処理によるヒト血管形態モデルの構築
はじめに医用画像を用いた形態構築に関して説明する。所定の倫理手続きを経て大阪大学医学部
より提供されたヒト脳の 4D-CTA 画像を用い血管形状の構築を行った。使用した画像は、約 1 心拍分を 22 フレームに分割して撮影されており、空間分解能は 0.468 mm×0.468 mm、解像度は 512×512、スライス幅は 0.5 mm、スライス枚数は 280 枚であった。4D-CTA では、静脈から注射された造影剤が脳血管内を循環する様子が撮影されるため、段階的に輝度値が高い塊の領域が画像として得られる。この特性を利用し、時相間で画像を差し引きすることにより、動脈系および静脈系の分離を行った。フレームごとの造影領域の抽出形状および本アルゴリズムの概念図を以下に示す。

図 2.9: 時系列造影剤画像を用いた動静脈の分離アルゴリズムの模式図

図 2.8: 造影領域の抽出形状の時系列変化

画像処理ソフト Amira(Visage Imaging, Berlin, Germany)を用い二値化処理を行い、本分離アル
ゴリズムにより抽出された血管形態を以下に示す。抽出結果より、動脈および静脈がある程度の血管径まで抽出され、明確に分離されていることが確認できた。

41

図 2.10: 分離アルゴリズムを用いて抽出された動脈(赤色)と静脈(青色)の血管形状管径まで抽出され、明確に分離されていることが確認できた。

(b) 離散的血管生成アルゴリズムによる脳表層血管形態の構築
時系列血管造影画像から画像処理により構築した血管形態モデルを基に、脳表層を巡る細い血管を数理モデルにより構築することを試みた。
本研究では、同サンプルの計測画像は CT 計測によるものしかないため、脳実質を直接抽出するのは困難である。そこで、頭蓋の内側領域を抽出し縮小することで擬似的な脳実質として表面形状モデルを構築した。これらは画像処理ソフト Amira(Visage Imaging, Berlin, Germany)を用いて行われ、最終的に疑似的な脳実質表面が三角形メッシュにより構築された。
構築された脳実質表面メッシュを利用し、以下の血管生成アルゴリズムを提案した。
 (a)で得られた脳血管メッシュを用い、脳実質表面の各三角形メッシュにおいて血管生成の不可に関するマスク処理を行う。脳実質表面メッシュと脳血管メッシュ間の重心座標の最短距離があるしきい値を下回った場合、その表面メッシュは画像処理により構築された大きな血管系が存在するとし、数理モデルによる血管生成の対象からはずす。
 (a)で得られた脳血管形態に対し細線化処理を行い、各末端血管の中心点位置座標を求めておく。
この位置座標と最も距離が近い脳表面メッシュを選び出し、数理モデルにおける血管生成を開始するメッシュとする。これらは動静脈でそれぞれ複数存在する。
 設定した開始点（開始メッシュの重心）から血管の分岐生成を行う。数理モデルの全ての開始
点を同一世代の親点とし、それぞれの親点は三角形要素の辺を共有する隣の三角形要素の格子点が空いていれば、1 つもしくは 2 つの隣り合う重心点へ子世代の点（子点）を選び、血管を分岐生成する。今度は子点が同様のアルゴリズムで次の世代の点を選び、これを繰り返し行うことで、重心点を埋め血管生成を行う。分岐生成アルゴリズムの概要を以下に示す。赤、青それ
ぞれ動脈の末端、静脈の末端から発生した数理血管を表している。赤血管の末端点 p、q がそれ
ぞれ隣り合うメッシュ点はどちらも 2 つある。このとき、q 点は 2 方向自由に選べるが、p 点は子点を選択できないため分岐生成プロセスはストップする。世代を進めて分岐生成を繰り返し、全ての血管の分岐生成プロセスがストップした時点でアルゴリズムを終了する。
42

図 2.11: 離散的血管生成アルゴリズムの概要
以下に、本アルゴリズムを用いて構築された数理血管、医用画像から構築された血管、脳実質表面を重ねた図を示す。赤色が動脈系、青色が静脈系を示している。脳実質表面に沿って全体を覆うように数理モデル血管が生成されており、主要な動脈・静脈から脳実質を張りめぐる細かな血管までモデル化された全脳動静脈モデルを構築することができた。このモデルは脳全体の血流動態を解明するための血流解析シミュレーションに今後用いることができる。
現在の問題点として、本アルゴリズムは血管同士の重なりを許していないため、脳実質表面メッ
シュが荒い場合、動脈系と静脈系の並走部分（動脈と静脈で入り組んだ構造となっている部分）が極端に少なくなることが挙げられる。そのため、表面メッシュをより細かくするとともに、生成された末端の血管の消去、再分岐を行うリモデリングアルゴリズムを加えることで、動脈と静脈がより入り組んだ数理モデルを構築できると考えられる。また、現段階では仮想的に血管を配置しているが、後述する脳血管網の計測データから得られる動静脈血管の配置や血管密度、血管径の分布などを反映させたモデリングが可能となっている。

図 2.12: 離散的血管生成アルゴリズムにより構築された数理血管、医用画像から構築された血管、
脳実質表面を重ねた図（赤色が動脈系、青色が静脈系を表わす）

43D) マウス大脳における脳微小血管構造の 3 次元計測と再構成法（電気通信大学

正本）(a) 目的全脳循環代謝シミュレーションを行うにあたり、脳微小領域における血管構造の詳細を把握するための微小血管構造の撮像および計算空間上へのリアリスティックな血管のネットワーク構造を再現するための 3 次元再構成法の検討。
(b) 方法
実験には成体マウス 2 例を用いた。脳血管を造影するために蛍光色素スルホローダミン 101 (10mM in saline、 シグマアルドリッチ)を画像撮像前に、イソフルラン麻酔下でマウスの腹腔に投与（8µL/g body weight）し、生体二光子顕微鏡を用いて左頭頂部の脳微小血管構造を撮像した。マウス大脳への光の透過性を確保するためにマウスの頭蓋の一部（直径 3 mm）を除去しカバーグラスで置き換える頭窓をあらかじめ作製した。
励起波長は 920 nm に設定し、25 倍の水浸対物レンズを通して、1024 × 1024 ピクセル（解像度 0.5 µm/pixel の画像サイズで 4 か所モザイク状に約 1 mm 四方の領域において二光子顕微鏡撮像を行った。各撮像領域において脳表から深さ方向へ約 1.2 mm までを
1.0～4.0 μm のスライス間隔で撮像を行った。得られた画像は Matlab 上で再構成し、撮像深さ間でのコントラストの調整およびノイズ除去を行い、4 か所の撮像画像を 1 つの画像に結合した。本研究で行った動物実験のプロトコルおよび実験計画に関しては、電気通信大学動物実験倫理委員会による承認を得た上で行った。
(c) 結果
生体二光子顕微鏡法を用いて撮像したマウス大脳における脳微小血管構造に関する 3 次元再構成
画像の一例を図 2.13 に示す。脳表から皮質に連絡する細動脈および細静脈を明瞭に観察し、皮質内で連続的に接続する毛細血管のネットワーク構造画像を得た。深さ間での画像コントラストの調整
およびノイズ除去後に得た信号背景輝度の比から、脳表から最深部約 1 mm の深さにおける血管構造を描出できていることを確認した（図 2.14）。
(d) 考察および今後の展望
本研究結果により大脳 1 ミリ立法の微細領域における血管構造画像を得ることが可能になった。
また撮像画像を計算空間上に 3 次元再構成することで、脳循環シミュレーション用のリアリスティ
ックな血管構造画像を提供することが可能である。ただし、本研究で得られた血管領域は脳全体の一部であり、今後は深部脳のイメージング技術の開発および領域間での血管構造の差異について検討を進める必要がある。さらにこのような脳の血管構造は同一領域においても成長や加齢とともに変化する可能性がある。特に脳血流の低下は認知機能の低下に先立つことが知られていることから、脳血管の微細なネットワーク構造の変化と脳病態の発症や脳機能の可塑的な変化との関連についてさらに検討を進める必要がある。同時に、血管周囲のグリアあるいは神経細胞の機能活動を画像化することが可能になれば、時間空間的な事象の因果関係を明らかにすることが期待される。

44

図 2.13: マウス大脳における脳微小血管ネットワークに関する 3 次元再構成画像の 1 例（画像の上向きに皮質表面方向、下向きに皮質深部方向を示し、大脳皮質の約1 mm 立方の領域における脳微小血管ネットワーク構造を生きたまま撮像することが可能になった）

図 2.14: 脳表からの深さ方向における血管領域の平均ピクセル輝度(signal)と背景輝度(noise)(2) データ同化血流解析（阪大和田Ｇ）

A) 物理的整合性を考慮したフィードバック制御によるデータ同化
患者個別の計測速度データを反映した血流解析の実現に向け、フィードバック制御に基づくデータ
同化血流解析の基礎的方法論を確立しその妥当性を検証するとともに、実際の患者の脳動脈形状および PC-MRI 計測データを用いた解析を行いアプローチの実現可能性を検討した。脳循環の場合、主要な動脈においても PC-MRI による血流計測には 20%程度の誤差が含まれるため、こうしたデータ同化
手法は、脳動脈瘤など様々な脳血管障害に関与する血行力学的因子を解明していく上で極めて有効な手段となる。

45

(a) 方法論の基礎的検討と確立
数値流体解析(Computational Fluid Dynamics: CFD)および医用計測により得られる速度の差を偏差とし、これに由来する外力を入力としフィードバック制御を行った。フィードバック制御による先行アプローチとして、速度差に由来する外力を運動方程式に直接作用させる方法が提案されているが(Funamoto et al., Ann Biomed Eng, 33, 2005, 415-428)、実際の流れには存在しない外力を考慮することで、非物理的な流動加速や渦度生成を引き起こす可能性があり、これらは誤差レベルが高い計測速度を用いる場合において重大な問題となる。そこで本研究では、フィードバック制御
における入力による影響を圧力境界条件に組み込むことで、物理的整合性を満足するアプローチを新たに提案する。
支配方程式として非圧縮性 Navier-Stokes 方程式を適用した。

 v  0(2.1)

r (¶t v + v × Ñv) = -Ñp + hÑ2 v(2.2)

血流は出入口を伴う閉空間の流れであるため、境界条件として壁面で速度固定、出入口面で圧力境界条件を適用した。

on Γ w
v  0  pn   n v   Pin n on Γin  pn   v   P n on Γ
n
out
out


(2.3)

ここで、v は速度ベクトル、p は圧力、Pin、Pout は入口・出口に課せられた境界圧力（応力ベク
トルの法線成分）である。
フィードバック制御として、CFD の速度ベクトル v および計測の速度ベクトル U を用い偏差ベク
トルを e=U-v と定義し、入力ベクトルを f=Kp*e と定義した、Kp は比例ゲインである。本研究では、二つの仮定を用い、入力ベクトル f による寄与を境界圧力 Pin、Pout に担わせた。
仮定 1：圧力境界により駆動される流れ場において、運動方程式の外力としての入力ベクトル f はス
カラーポテンシャル q のみによる保存力として表現される。すなわち、f=-q となる。
仮定 2：圧力境界により駆動される流れ場は、仮定 1 で定義されるスカラーポテンシャルを境界圧力に用いることで表現される。すなわち、それぞれ出入口境界において、Pin=q、Pout=q となる。
ベクトル恒等式よりスカラーポテンシャル場は以下の式を解くことで求めた。
 2 q    f(2.4)

なお、境界条件は Neumann 条件とした。
並列化の容易さを考慮し、直交格子を用いた有限差分法により離散化を行った。この際、任意形状を Boundary data immersion 法(Weymouth, J Comput Phys, 230, 2011, 6233-6247)により扱う。
解法にはフラクショナルステップ法(Dukowicz & Dvinsky, J Comput Phys, 102, 1992, 336-347)を適用した。解析には開発中のボクセル型流体・固体解析プログラム（VBB-Flow）を用いた。
提案アプローチの妥当性を検証するために、ハーゲンポアズイユ流れに対する数値テストを行った。ここで、計測データを模擬するために任意の空間解像度を指定し等間隔に疑似データを配置し、46

その点におけるポアズイユ流れの解析解に正規乱数に従うノイズを加えることで速度ベクトルを得た。ノイズレベルとして SNR(signal-to-noise ratio)を正規分布の標準偏差の逆数と定義し、これを指定することで任意のノイズレベルを持つ疑似計測速度を作り出した。以下に SNR=20 とした疑似
計測速度データおよび推定された速度場の結果を示す。ポアズイユ流れの特徴である二次分布が良く再現されていることから、データノイズに対するロバスト性が示された。数値シミュレーション
では、計測の空間解像度より細かい解像度を採用し解を得ることが可能であり、壁面せん断応力など局所速度場に依存する物理量を評価することが可能となる。推定速度場より評価された壁面せん
断応力(Wall Shear Stress: WSS)と厳密解の正規化誤差 Ewss とノイズレベルとの関係を調べた結果、
Ewss=0.0609（ノイズなし）、0.0617（SNR=100）、0.0836（SNR=20）、0.101（SNR=10）となり、
ノイズレベルの上昇に従い誤差が大きくなることが分かった。しかしながら、SNR=10 においてもその誤差は 10％ほどであり、また誤差が大きくなる箇所は出入口付近であったため、実用的な解析が可能であることが示された。

図 2.15: ノイズを含む疑似計測

図 2.16: データ同化解析による全速度ベクトルの様子

格子点での速度ベクトル

図 2.17: 主流方向速度に関する
厳密解との比較(SNR=20)(b) 医用計測データを用いた解析
内頸動脈を対象とし、血管形状を CT 装置により取得した DSA(Digital Subtraction Angiography)
画像から画像処理ソフト AMIRA5.4.2(Visage Imaging, Berlin, Germany)で構築し、計測速度ベク
トルを MRI 装置により取得した Phase Contrast-MRI(PC-MRI)画像から構築した。この際、CT 計
測と MRI 計測の手技の違いにより位置ずれが生じるため、剛体の並進・回転に伴う 6 自由度成分を離散的に調査し、形状を評価するパラメータを最小にする値を採用することで位置ずれを修正した。
DSA 画像の解像度は 0.1 mm、PC-MRI 画像の解像度はおよそ 0.5 mm であった。PC-MRI 計測では一心拍あたり 15 フレームの速度情報が得られるが、ここでは、最大収縮期をとる時相を採用した。
データ同化シミュレーションにおいて、ボクセル格子幅をおよそ 0.2 mm に設定した。流れ場が収束するまで時間発展を繰り返して定常解とした。異なる患者および部位における脳動脈瘤 10 例に対し解析を行った。なお、医用計測は大阪大学が定める所定の倫理手続きを経て大阪大学医学部により行われたのち、匿名化を行った画像が提供されている。
二例の解析に関して、データ同化解析に用いた PC-MRI 速度ベクトル、データ同化により求められた全 CFD 格子点での推定速度ベクトル、これを用いて評価された壁面せん断応力を以下に示す。
結果より、少ない計測データを用いた場合でも、CFD で採用した高解像度場において、非物理的な

47

振動なく速度場が推定されている様子が確認でき、それにより評価される壁面せん断応力も生理学的に妥当な範囲で得られていることが確認できた。本二例において、計測点における推定速度と計
測速度の入口平均速度で正規化された残差はそれぞれ 20%と 15%となった。この残差の違いを生み出す要因を調べるため、
全 10 例での入口断面の等価直径と残差の関係をプロットし確認したところ、等価直径の上昇に伴い残差が減少するといった負の相関関係が得られた。よって、計測データ数あるいは管形状に対する相対的な PC-MRI 解像度が推定精度に影響を与えることが示唆された。

PC-MRI 速度ベクトルの様子

推定速度ベクトル

壁面せん断応力

PC-MRI 速度ベクトルの様子

推定速度ベクトル

壁面せん断応力

図 2.18: 医用計測データを用いたフィードバック制御に基づくデータ同化血流解析の結果(c) アプリ開発状況
開発中のボクセル型流体・固体解析プログラム（VBB-Flow）を拡張することで派生アプリ
VBB-PFCDA を開発した。領域分割による MPI 並列において、各領域に含まれる計測点を分割コアに割り当て通信を行った。8 プロセス(MPI)×16 コア(OpenMP)のハイブリッド並列環境において、本研究で示した定常速度場を得るのに 20 分程度であり、実用問題に適用した場合でも良好な並列性能を達成している。
B) 三次元変分法によるデータ同化(a) 方法論の基礎的検討
患者個別の脳動脈における血流場および関連する壁面せん断応力などの血行力学因子の解明は瘤破裂と関連し、患者個別速度場を反映することが可能なデータ同化血流解析は非常に強力なツールとなる。本研究では、脳動脈瘤付近の血管形状および計測速度場を利用することで、無駄な演算を省いたデータ同化解析アプローチの可能性を模索した。

48

フィードバック制御に基づくアプローチでは推定速度と厳密解の誤差が最小となる保証がない。
一方、変分法によるデータ同化解析では、適当な分布仮定のもつ各データ群間の分散を最小にする
ことで得られる新たな解は厳密解に最も近くなることが数学的に保証されている。脳動脈瘤付近の一部のデータを用いて解析する際、この最小解の保証が推定精度に影響を与えると考えられる。そこで、三次元変分法をベースにアプローチを確立した。本手法では、非定常非圧縮性 Navier-Stokes
方程式を拘束条件（付帯条件）とし、モデル速度 v と計測速度 U の差が最小値をとる出入口境界圧
力 g および h を求める最適化問題を設定した。これらは以下のように記述される。

arg min P ,P (v, g,h) =
in

out

2
2
2
1
e
a
Q(v) - U + Ñ 2s g + Ñ 2s h + ¶t v
2
2
2
2
2
2

ì
ï Ñ×v = 0
ï r (¶t v + v × Ñv) = -Ñp + hÑ 2 v
ï
í v=0
ï
ï - pn + h¶n v = -gn
ï - pn + h¶ v = -hn (h = 0)
n
î

in Ω
in Ω
on Gw
on Ginon Gout(2.5)

上式の評価関数において、第一項は誤差項、第二項はノイズによる悪条件性を緩和する Tikhonovの正則化項、第三項は慣性による解の振動を抑えるペナルティ項である。計測速度は定常速度であるため、上式を時間発展的に解き、血流場が定常状態になるまで行った。
(b) 数値実験による検討
本解析では、デジタルサブトラクョン血管造影(DSA)画像から脳動脈瘤形状を抽出した。脳動脈瘤の抽出には、汎用ソフト Amira5.4(Maxnet)を用いた。抽出方法として、脳動脈瘤付近での輝度値の闘値を設定し、瘤とそれ以外の生体組織で二値化画像を作製し、この画像を基に三角形要素で表現した表面形状を作成した。データ同化解析では、計算格子の格子幅は 0.2 mm、時間刻みは 0.5 ms、また重み関数αは 10 とし、正則化パラメータεは L-カーブ法を用いて決定した。
適当に g、h を定め求解した結果を準厳密解とし、これを基にノイズが含まれない疑似計測データ(Noise-free)および正規乱数に従うノイズが含まれる疑似計測データを作成し、同化解析に利用した。
Noise-free データおよびこれを用いた計測位置での同化解の速度ベクトルの分布を以下に示す。
ここで、ベクトルの長さを全てそろえて表示している。同化解の速度ベクトルは Noise-free データと定性的に一致した。このとき、 疑似計測データとの残差 Em および準厳密解との誤差 Eref はそれぞれ 9.8E-3 および 1.1E-2 であり、同化解の速度場は基準解と定量的に一致した。データ同化解
析の解像度は疑似計測データの解像度の 3 分の 1 であるので、提案手法は計測データよりも高解像度で脳動脈瘤内の血流場を取得できる。次に計測のノイズレベルを変化させた際の Em および Erefの変化を示す。疑似計測データとの残差 Em は計測データに計測ノイズが含まれるので、誤差レベ
ルの増加に伴いその分だけ誤差が大きくなっていった。その一方、計測データの元となった基準解と同化解との全計算点での誤差 Eref は十分に小さかった。SNR=20 の条件では、Eref は 5.0E-2 以下であり、Em と比較しても 4 分の 1 以下であった。これらより、提案手法は脳動脈瘤形状のような複雑な血流でも計測ノイズに対するロバスト性を有することが示された。

49

図 2.19: 疑似計測速度と推定速度の様子(Noise-free)

図 2.20: 疑似計測データとの残差 Em と準厳密解との誤差 Eref に関するデータノイズの影響

脳動脈瘤の壁面せん断応力(WSS)の分布を以下に示す。(i)は準厳密解、(ii)は Noise-free データでのデータ同化解析の結果、(iii)は Noise-free データに SNR = 20 の計測ノイズと位置ずれの誤差 δX =0.4 を加えた計測データでのデータ同化解析の結果である。(ii)と(iii)での WSS の分布は準厳密解と定性的に一致した。このことから、計測誤差を含む計測データを用いた場合でも正確に WSS を推定できることが確認できた。なお、一般的に、PC-MRI の計測では血管壁近傍の速度データの誤差が大きくなる。本手法を実際の計測データに提供する場合には、採用する計測データに重みをつけるなどの改善が必要であろう。

図 2.21: 動脈瘤内の壁面せん断応力分布。準厳密解(i)と Noise-free データを用いた推定結果(ii)およびノイズデータ(SNR=20)を用いた推定結果(iii)。
(c) アプリ開発状況
開発中のボクセル型流体・固体解析プログラム（VBB-Flow）を拡張することで派生アプリ
VBB-VADA を開発した。全ボクセル格子において流体部のみを配列に再定義することで if 文による負荷を軽減し、これらに OpenMP による共有メモリー並列化を行った。16 コア並列環境において、本研究で示した定常速度場を得るのに 20 分程度であり、瘤付近を対象とする実用問題に適用した場
合でも良好な並列性能を達成している。

50

(3) 脳動脈瘤コイル塞栓治療の治療支援に向けた計算力学フレームワーク構築（阪大 和田 G）
コイル塞栓術は脳動脈瘤の破裂を予防する血管内治療法の一つであり、カテーテルを通して金属製のコイルを瘤内へ留置し、血流を滞留させ瘤全体を血栓化する。これにより瘤へ流入する血流を塞栓し、その破裂を未然に防ぐ。本研究では、コイル留置による脳動脈瘤塞栓の機序を明らかにし、患者個別に適切なコイルや留置手順を選択する治療支援の創成を目指す。今年度において、まず、コイル塞栓時の瘤内血流を安定かつ高精度に解析するための数値流体計算のフレームワークを構築した。本フレームワークを用いて、コイルの分布形状の違いが瘤内血流の滞留度合に与える影響について検討した。さらに、
コイル塞栓時の瘤内における血栓形成を現象論的に表現する数理モデルを構築し、瘤内の血栓化による塞栓効果をコンピュータ上で表現することに成功した。以下にこれらの概要を述べる。
A) 瘤内のコイル分布形状が血流動態に与える影響
コイルが留置された脳動脈瘤の表現にあたり、CT 画像から構築した患者個別の瘤形状を用い、さらに先行研究（Otani et al., Med. Biol. Eng. Compt., 2016）にて構築したコイル挿入シミュレーショ
ンを実施することでリアリスティックなコイル形状を構築した。コイルが留置された瘤内の血流動態を安定かつ高精度に解析するため、コイル留置後の瘤内の流体領域を Volume fraction (VOF)法により直交格子中に投影し、有限差分法による血流の数値流体計算を行った。壁面境界条件には BDI 法(Weymouth and Yue, J. Compt. Phys., 2011)を用い、Fractional step 法に基づく分離解法により非圧
縮性 Navier-Stokes 方程式と連続の式を連立して解いた。入口境界に生理学的に妥当な流量を設定し、出口境界の圧力を一定とした。ここで計算の簡略化のため、入口境界の速度は一定とした。計算結果の一例を図 2.22 に示す。解析結果が先行研究の解析結果と良好に一致したことを確認した。従来の数値流体計算は主に非構造格子を用いた有限要素法・体積法による解析が主であり、コイルが密に留置された複雑な流体領域に対し非構造格子の構築が必須なため、前処理のコストが膨大になることが問題になっていた。本手法は直交格子法を用いたことで数値流体計算の前処理を大幅に自動化することに成功し、これにより、多数の症例や様々なコイル分布形状を用いての大規模な血流解析の実施が現実的に可能となった。
構築した手法を用いて、瘤内部におけるコイルの分布形状が瘤内血流の滞留度合に与える影響を検討した。コイルの挿入条件を変えて挿入シミュレーションを行い、複数の異なるコイル分布形状（n=102）を得た。ここで挿入するコイルは全 102 例において同じものとし、コイルの充填率が臨床において適切な値の範囲 (>25%) である 27%となるよう設定した。102 例についてそれぞれ血流解析を行い、コイル留置前後における瘤内血流の運動エネルギー・せん断率の体積平均の変化率をそれぞれ取得した。臨床において予後の病態に影響すると指摘されているコイル分布の指標として、瘤のネック領域における局所的なコイルの充填率に着目し、瘤と母血管との境界面から 0.5 mm の区間の瘤
領域体積におけるコイルの充填率をネック領域のコイル充填率として定義した。このネック領域のコ
イル充填率に対して、瘤内全体における血流の運動エネルギー・せん断率との相関を調べた。それぞれの結果を図 2.23 に示す。瘤内全体のコイル充填率を 27%とした場合でも、ネック領域のコイル充
填率は 5%から 35%程度までばらつきが生じており、このネック領域のコイル充填率と、血流の運動
エネルギー・せん断率には負の相関が生じた。特にせん断率について、値は 20%から 70%程度まで大きくばらつきが生じた。この結果は、コイルが密に充填されている場合でも、ネック領域におけるコ51

イルの留置が不十分であれば、運動エネルギー・せん断率ともに十分に減少せず、瘤内血流の塞栓が達成されない可能性を示唆している。以後は他の症例を用いて同様の解析を行い、コイル分布が瘤内
血流の滞留度合に与える影響について症例ごとの違いを検討する。

図 2.22: コイル塞栓時の瘤内血流の流線図（コイル充填率 27%）
。矢印は流れの方向を示す。

(a)R2=0.603
n=102

30

Flow shear rate ratio (%)

Kinetic energy ratio (%)

40

20

10

70

(b)

60

R2=0.673
n=102

50

40
30
20

0

5
10
15
20
25
30
Coil packing density in neck region (%)

5
10
15
20
25
30
Coil packing density in neck region (%)

図 2.23: コイル塞栓前後における瘤内の運動エネルギー（a）、せん断率(b)の変化に対するネック部の局所的なコイル充填率の寄与。
B) コイル塞栓時の脳動脈瘤内における血栓形成の現象論的モデル構築
前節にて構築したフレームワークを応用し、コイル塞栓術後の脳動脈瘤内における血栓の形成を現
象論的に表現する数理モデルを開発した。従来の血栓形成に関する生体外実験の知見を参考に、血流中のせん断率の低下に応じた血流中の血栓率   0,1 をシグモイド関数で表わすとともに、血栓率の増大により血流に生じる流体抵抗を、均質な多孔質体と等価な流体抵抗としてモデル化した。実際の数値計算においては、計算手順を二段階に分け、まず前節にて構築したコイル塞栓術後の瘤（充填率 27%）について、血栓による流体抵抗を考慮せず数値流体計算を行い、血流場の定常解を取得した。その後、得られた血流場に対して血栓の数理モデルを導入し、血栓形成とそれに伴う血流場の変化が平衡状態となるまで計算を行った。
コイル留置後の瘤内部断面について、血栓形成の数理モデルを導入した後の初期状態と平衡状態に52

おける血栓率の分布を図 2.24 に示す。初期状態において、コイルの留置により瘤内部の血流が滞留し、瘤の中心から奥にかけて血栓率が最初から大きくなった一方で、瘤のネックから中心にかけて高いせん断率が残存しており、血栓率はほぼゼロとなった。この初期状態に対して、血栓化が進行した平衡
状態では主に瘤の中心部の血栓率が大きくなった。しかし、ネック領域の血栓率は局所的に小さくなり、初期状態よりも血栓率が小さくなった領域が見られた（図 2.24 右赤枠）
。これは瘤内の血栓形成を考慮することで初めて観察された現象であり、血栓形成と血流動態との相互作用により生じたと考えられる。今後は数値流体計算による血流解析だけでなく、生理学的知見を取り入れた血栓形成モデ
ルを組み合わせた本解析手法の精緻化を進め、コイル留置から瘤塞栓に至る全体の機序解明を目指す。

図 2.24: 初期状態（左）と平衡状態（右）における瘤内部断面の血栓化率。
(4) 大規模ライブセルイメージング解析技術の開発（阪大 松田 G）
血管中の細胞の動態を非侵襲で長時間観察できる蛍光顕微鏡で得られる膨大な時系列画像データを高精度かつ高速に処理するための、大規模ライブセルイメージング解析技術を「京」を使って開発している。
和田グループの開発している全脳循環代謝シミュレータでは、脳血管障害や脳機能障害の解析プラッ
トフォームの構築を目指しているが、主要な脳血管障害の一つである脳血栓症の発症メカニズムの解明
では、慢性的に進行する動脈硬化とそれに伴う血栓の形成の過程のシミュレーションが重要と考えられる。
そこで、本研究では、生体内部を精密かつ経時的に観察できる二光子励起顕微鏡を使って、血管内のライブイメージングデータを取得し、動脈硬化の進行に伴う血管内部の変化をイメージング解析により解明することで、脳血栓のシミュレーションのデータ同化に必要な計測データを取得することを目指し
ている。

53

図 2.25: 炎症刺激の前（左）と刺激後（右）の生体の血管内の白血球の観察画像（提供 石井優教授（大阪大学生命機能研究科））

図 2.27: 白血球を認識する隠れマルコフモデル（p は状
態遷移確率、y は特定の輝度の出力確率）

図 2.26: 開発したライブセルイメージング解析技術の概要

近年、動脈硬化は炎症性病態であることが明らかとなりつつある（P. Duewell et al., Nature, 2010 など）ことから、動脈硬化の進行を炎症を誘発する刺激剤（LPS、リポ多糖）の投与（以下、炎症刺激と略す）による血管内での炎症進行の経過でモデル化している。炎症刺激後の血管内皮の状態変化を、緑色蛍光で標識した白血球（炎症時に血管内皮の細胞に接着・浸潤することが知られている）の動きの変化として観察している。
平成 28 年度に取得したライブイメージングデータは、マウス生体の微小血管中を二光子励起顕微鏡で撮影した経時観察動画像であり、毎秒最大 30 フレーム、各フレーム 512×512 ピクセル、24 ビットカ
ラーであり、1 回の観察ごとに 2000 フレーム（毎秒 30 フレームだと約 67 秒）の画像を取得した。血管の部位や炎症刺激後の経過時間を変えて 200 回以上撮像しており、総データ量は約 40 万フレーム（約137G バイト）となっている。図 2.25 に刺激剤投与前（左）と投与による炎症刺激後 9 分時点（右）での同じ血管内の観察画像を示す。図 2.25 で赤く染まっている領域がマウス生体内の血管であり、画像中54

で血管とその周辺の領域のみ表示している。緑色で表示されている領域が白血球であり、炎症刺激後は白い丸で囲った箇所のように特定の領域に白血球が凝集している現象が見られ、血管内皮への接着が生じていることが示唆される。
図 2.26 に今回開発した大規模ライブセルイメージング解析技術の概要を示す。取得した経時観察動画
像は RGB の 3 チャンネルからなるカラー画像であるが、緑色蛍光で標識された白血球は G（緑）チャン
ネルでしか現れないため、まず緑チャンネルの画像データのみを抽出する。次に、その画像データから白血球の領域を検出する。これは原理的には輝度の大きい領域を検出すればよいはずであるが、白血球によって輝度値にばらつきがあり、またノイズがあるため、特定の閾値で切ると誤検出が生じる。そこで本研究では、画像中の各ピクセル（画素）が白血球か背景かを判別する隠れマルコフモデルを構築した（図 2.27）。これによりフレームごとの白血球の通過による各ピクセルの輝度値の変化をモデル化する
ことで、あるフレームでそのピクセルが白血球の領域か背景かを判別することができる。
白血球の移動量（速度ベクトル）の計算には、ディープマッチング（P. Weinzaefel, Intl. Conf. Comp.
Vision, 2013）を使用した。ディープマッチングは、2 枚の画像間で類似した輝度を持つ領域を照合する
アルゴリズムで、領域間の対応を 8×8 ピクセル程度の小領域から始めて段階的に画像全体にまで行う方法である。領域間の対応付けを畳み込み計算により行うことから、対応する領域間で形の変形があっても照合できる特徴を持っている。白血球は剛体ではなく変形しながら血管内を移動するため、オプティ
カルフローなどの固定された領域どうしを照合する手法ではうまくできないことがあり、ディープマッ
チングを使うことで高精度に照合することができる。しかし、ディープマッチングは、畳み込み計算を階層的に繰り返すことから計算量が非常に大きく、1 回（隣接する 2 フレーム間）の移動量計算に約 27秒かかるため、30 回分（6 万フレーム、観察時間は合計約 33 分）のライブイメージングデータでさえ、計算に約 450 時間もかかることになる。本研究では、ディープマッチングについて、「京」で OpenMP（スレッド並列）と MPI（プロセス並列）を組み合わせたハイブリッド並列化により高速化を実現した。
その結果、「京」の 300 ノード（2400 コア）で 6 万フレームの移動量計算時間を約 24 分と、観察時間以下にまで短縮することができ、逐次計算と比べて 1000 倍を超える速度向上を達成できた。前述のように、現時点で保有しているイメージングデータのみでも 40 万フレーム以上あるため単純に計算すると約3000 時間かかることになるが、「京」を使うことで約 2.7 時間まで短縮できることになる。
ディープマッチングによる移動量計算では、画像全体で領域の照合を行うため、前述の白血球の領域
認識と組み合わせることで、白血球の領域のみの移動量を求めることができる。計算結果の例を図 2.28に示す。図 2.28 では、求められた白血球（緑色の領域）の移動の方向と大きさを白い矢印で表示している。

図 2.28: ディープマッチングによる白血球の移動量計算の例

図 2.29: 白血球の移動量の頻度分布55

次に、白血球の移動量の分布を、炎症刺激の前と、刺激後 4 分および 9 分のそれぞれの時点で集計したものを図 2.29 に示す。
図に示すように炎症刺激により白血球の血管内皮への接着が見られることから、炎症刺激の前と比べて、炎症刺激後の時間経過とともに移動量が減少していくことが想定されたが、図2.29 ではむしろ炎症刺激後の方が、移動量が大きいものが増えている。この原因としては、炎症刺激によりその部位に白血球が数多く集まるため、血管内皮に接着するものだけでなく、血管内を速く遊走す
るものが多く現れたことが考えられる。これを判定するために、次年度では、時間的に隣接するフレー
ム間での移動量の計算に加えて、より長い時間区間での白血球の追跡（トラッキング）を行うことで、接着と遊走の白血球を分けて移動量を解析する手法の開発に取り組む予定である。
(5) 血栓症シミュレータの開発と血小板粘着現象解析（東大 高木 G）
本研究は、心疾患・脳血管疾患の一因となっている血栓症を対象としたものである。血栓症を再現可
能な数値シミュレーションツールを開発し、原因解明、疾患予測・予防、治療支援、新薬開発などに貢献することを最終的な目的としている。血栓形成は多くの因子が複雑に絡み合い連鎖的に起こる事象であるが、本年度は特に血栓形成の初期段階に相当する血小板の粘着に着目した研究を行った。血小板粘着は、血小板表面上の膜糖タンパクの一種である GPIb alpha (Glycoprotein Ib alpha)と血中に存在する凝固系因子の一つである VWF (von Willebrand factor)との結合を介して実現される。VWF は血管壁損
傷部位に露出したコラーゲンと結合し、血小板と血管壁との粘着に介在する。
開発アプリの概要を述べる。解析手法はこれまで開発されてきた流体構造連成解析手法[1][2]に改良を加え Immersed Boundary 法[3]を用いた流体構造連成解析手法とした。構造物となる赤血球・血小板は超弾性体として扱っている[4][5]。さらに、血小板粘着を取り扱うために必要な GPIb alpha-VWF の結合に関する力学モデルを導入している。この力学モデルでは、遷移状態理論に基づいた反応速度定数[6][7]から確率的に各時間ステップにおける結合・解離を判定する。プログラムは OpenMP-MPI ハイブリッド並列プ
ログラムとなっており、MPI 並列には V-Sphere ライブラリ[8]を利用している。

図 2.30: 血流解析の時間発展の様子（Ht = 0.3, t = 90 ms）．左図:赤血球条件、右図:剛体球条件、赤:赤血
球、黄:血小板、青線:GPIb alpha - VWF 結合を表す。

56

Number of adherent platelets

500 Line
Solid: Red blood cell
Dashed: Undeformable sphere

400 Color
Ht = 0.3
Ht = 0.2
Ht = 0.1

300 Ht = 0
200
100
0
0

20

40
60
Time [ms]80100

図 2.31: 粘着血小板数の経時変化
表 2.1: ヘマトクリットと平均流速・相対実行粘性係数の関係
赤血球条件

剛体球条件

Ht

平均流速 [cm/s]

相対実行粘性係数

平均流速 [cm/s]

相対実行粘性係数

0

2.50

1.00

2.50

1.00

0.1

2.02

1.24

1.56

1.60

0.2

1.83

1.37

1.16

2.15

0.3

1.49

1.68

0.729

3.42

つづいて代表的な結果を述べる。本稿で述べる結果については、2048×512×512 の計算セル数を用
い、「京」 2048 ノードを利用して行った。解析は 0～30% のヘマトクリット値（Ht）に対して行い、そのときの粘着血小板数の変化に着目している。赤血球に対する解析と同時に等体積の剛体球に対する解析も行い、赤血球の変形の影響等についても併せて考察した。図 2.31 に、得られた結果の一例として、
赤血球（剛体球）
・血小板挙動の様子を示す。初期条件としてランダムに配置された赤血球と血小板が流れの発達とともにどのような振る舞いをしているか確認することができる。赤血球は初期形状から徐々に変形し、流れに沿った扁平形状となる。それとともに流路中央付近に集まる傾向が見られ、赤血球は壁面近傍よりも流路中央付近の方が、数が多い。剛体球条件に比べると変形の様子や流路中央付近への集積の様子がよくわかる。一方、血小板は相対的に大きな赤血球・剛体球の間を移動しながら、赤血球が流路中央付近に集まる傾向とは反対に壁面へと押し出されるように動く。その中でも壁と接触した血小板は GPIb alpha - VWF 結合により壁面に粘着する。表 2.1 に得られた数値計算結果から求めた平均流速および相対実行粘性係数を示す。赤血球・剛体球条件ともにヘマトクリットが高くなると、平均流速が遅くなり相対実行粘性係数が大きくなることがわかる。この現象は血流においては
Fahreus-Lindqvist 効果[9]としてよく知られた現象であり、本解析においても実血液で観察される現象が再現されていることが確認できた。数値計算結果として得られた相対実行粘性係数の値は多くの実験結果[10]に比べやや小さめの値となっており、赤血球モデルの扱いなどに改良の余地があることを示唆するものと考えられ、今後の課題として挙げられる。また、赤血球条件と剛体球条件とを比べると、同ヘマ
トクリット条件の場合には赤血球条件の方が、平均流速が高く、相対実行粘性係数が小さい値に留まっていることがわかる。これは、赤血球の変形によって血液の実行粘性が小さくなることを意味しており、

57

赤血球の変形能が血流に与える影響を示唆する結果である。図 2.31 に粘着血小板数の経時変化を示す。
赤血球条件・剛体球条件ともにヘマトクリット値が高い方が、粘着血小板数が多くなることがわかる。
これは、赤血球・剛体球が存在することで血小板が壁面と接触する機会が増えるためであり、その効果は高ヘマトクリット条件ほど大きくなるからである。赤血球が無い場合には接着する血小板はほとんど観察されない。またこの図では、赤血球条件の場合にグラフの傾き（単位時間当たりの粘着血小板数）が計算初期と後半とで大きく異なっていることもわかる。一方で剛体球条件の場合には赤血球条件ほど顕著な差は見られない。これは赤血球の初期配置の影響と考えられ、流れの発達とともに赤血球が変形することで壁面方向速度の変動分が小さくなっていくことを反映したものと考えられる。これらの数値
シミュレーション結果は、血小板粘着には赤血球（血小板以外の物体）が物理的に存在し血小板を壁面方向へと移動させる効果が誘起されることが重要であり、初期の血小板粘着に及ぼす赤血球の影響において生化学的な効果はそれほど重要ではないということを示している。これらの結果は、比較対象となる実験結果[11]と定性的にも定量的にも妥当な範囲で一致することを確認しており、一連の解析を通じて初期の血小板粘着における赤血球の役割の重要性を示すことができた。同時に、本解析モデルの妥当性についても確認された。
以上、本稿では割愛した部分も含めて本年度までに得られた成果をまとめると以下のようになる。本
血流シミュレーションは、赤血球の変形や中央部への集積などの様子、ヘマトクリットと相対実行粘性係数の関係（Fahreus-Lindqvist 効果）、赤血球の変形による圧力損失の低減効果、などの実血液でこれまで観測されている現象を再現していることを確認した。また、血小板粘着現象に赤血球が及ぼす影響について、粘着血小板数は高ヘマトクリットほど多くなり赤血球が無い場合には血小板粘着はほとんど起こらないこと、赤血球が物理的に存在し血小板を壁面方向に移動させる流れを誘起することが血小板粘着を促進させること、一連の結果が全血を用いた実験結果と定量的に概ね一致していることから初期の血小板粘着現象における赤血球の影響は生化学的な効果を通してではなく流体力学的な効果を通してである、という重要な知見を得ることができた。
参考文献[1] Sugiyama, K. et al., J. Comput. Phys., 230(2011), 596-627.
[2] Ii, S. et al., J. Biomech. Sci. Eng., 7(2012), 72-83.
[3] C. S. Peskin, J. Comput. Phys., 10(1972), 252-271.
[4] R. Skalak, et al., Biophys. J., 13(1973), 245-264.
[5] Modeling and Simulations of Capsules and Biological Cells, Edited by C. Pozrikidis, Chapman &Hall/CRC, 2001.
[6] S. Shiozaki, et al., J. Biomech. Sci. Eng., 7(2012), 275-283.
[7] J. F. Marko et al., Macromolecules, 28(1995), 8579-8770[8] Ono, K., et al., Lecture Notes in Comput. Sci. Engrg., 67(2009), 9-18.
[9] Fahraeus, R. et al., Am. J. Physiol., 96(1931), 562-568.
[10] A. R. Pries et al., Am. J. Physiol. -Heart Circ. Physiol., 263(1992), H1770-H1778.
[11] 田村典子他, 第 39 回日本バイオレオロジー学会年会, 2016.

58

(6) ヒト運動機能の神経筋骨格系シミュレータの開発（阪大 野村 G, 東大 高木 G）A) はじめに本研究テーマでは、ヒトの運動制御メカニズムおよび神経疾患に起因する運動障害の発現メカニズ
ムの解明を目的として、大規模計算ニューロメカニクス研究の基盤プラットフォームの開発を行っている。本研究では、昨年度までスーパーコンピュータ「京」のプロジェクトで開発された、全身筋骨
格-神経系統合シミュレータ[1]（脳神経系シミュレータ、筋線維の集合体として筋肉全体の振る舞いを再現する骨格筋シミュレータおよび全身骨格系運動シミュレータの統合シミュレータ）をベースとして、新たに収縮タンパクの確率的な挙動から筋骨格系の運動を再現するためのマルチスケール・マル
チフィジックス骨格筋シミュレータの開発を行った。現在までの成果として、収縮タンパクの確率的な挙動から 3 次元骨格筋モデルの収縮挙動を再現する連成シミュレーションが可能となった。
B) マルチスケール.マルチフィジックス骨格筋シミュレータ
骨格筋は複数の細長い筋線維の集合体である。骨格筋線維は疲労耐性や収縮速度など、収縮タンパ
クの一つであるミオシンの性質の違いによって、遅筋や速筋などに分類される。骨格筋の遅筋線維と速筋線維の割合は求められる機能によって異なり、一般的に抗重力筋の場合には遅筋が支配的であるが、例えば最新のヒト姿勢制御の理論によれば、腓腹筋の速筋線維の素早い収縮・弛緩が立位維持に重要な役割を果たすことが指摘されている。これは、運動制御の研究にとって、詳細な筋モデルの開発の必要性を示す一つの例である。骨格筋筋線維の収縮は脊髄にある運動ニューロンにより支配されており、脳からの運動指令が運動ニューロンに伝えられると筋線維は収縮して力を発揮し、最終的に、関節運動が生じることになる。
本研究では、このような骨格筋の構造的な特徴を考慮するため、計算モデルとして図 2.32 に示す 3次元連続体モデル内に複数の 1 次元筋線維モデルを埋め込んだ階層構造モデル[2]を採用した。

図 2.32: 3 次元有限要素（黒色枠）に組み込まれた 1 次元筋繊維モデル（赤色）(a) 3 次元連続体モデル
骨格筋や腱などの軟組織は超弾性、非圧縮性、大変形、非線形性などの特性を有する。そこで本研究では骨格筋の 3 次元連続体モデルとして、非圧縮性超弾性体の 3 次元混合型有限要素法による定式化[3]を採用した。
骨格筋は筋収縮により力を発揮するという特徴を持つ。従って、筋線維の応力状態 Stotal は次式で与えられる。
S total  S pas  S act(6.1)

ここで Spas は非圧縮性超弾性体の変形から定義される受動的な応力である。一方、Sact は筋収縮に59

よる力の寄与分で、1 次元筋線維モデルの計算から決定される応力である。
(b) 1 次元筋線維モデル（筋収縮モデル）
運動指令が脊髄の運動ニューロンを通して筋線維表面に伝えられると、神経筋接合部に活動電位が発生する。活動電位は筋線維に沿って長手方向に、また、T 管を通して筋線維内部に伝播する。筋線維内部に伝えられた活動電位はリアノジン受容体の構造変化を引き起こし、筋小胞体の Ca2+リリースチャネルを開放する。細胞体内に筋小胞体から Ca2+が放出されると、Ca2+は収縮タンパクの一つであるアクチンフィラメント上のトロポニン C と結合し、アクチンフィラメントとミオシンヘッドの結合サイトを露出させる。その結果、アクチンフィラメントとミオシンフィラメントの間にク
ロスブリッジが形成され、ミオシンヘッドの首振り運動により筋収縮が生じる。
本研究では、これらマルチスケール・マルチフィジックスな筋収縮のダイナミクスを再現するた
めに、i) 筋線維に発生する活動電位およびその伝播のシミュレーション、ii) 筋小胞体の Ca2+ダイナ
ミクスのシミュレーション、iii) 収縮タンパクのクロスブリッジダイナミクスのシミュレーションの連成による骨格筋筋収縮モデルを開発した。
i. 筋線維の活動電位モデル
筋線維の活動電位は運動ニューロンからの電気的シグナルによって引き起こされ、筋線維に沿って伝播する。骨格筋筋線維は互いに絶縁されているため、活動電位の伝播は次式の 1 次元
反応拡散方程式[4]で記述される。
Vm 

 Vm

 I ion (Vm )   eff  Am  Cm

x 

t



(6.2)

ここで、 𝜎eff および𝑉𝑚 は筋線維の電気伝導度および膜電位である。𝐴𝑚 は筋線維の体積に対す
る膜表面積の割合、𝐶𝑚 は細胞膜の単位面積あたりのキャパシタンスである。イオン電流

𝐼ion

は細胞膜を通過するイオン電流の総和であり、膜電位に依存する。ここでは活動電位のモデル
として、筋鞘および T 管を考慮した、Hodgkin-Huxley タイプの電気生理学モデル［5］を用いた。
ii. 筋小胞体の Ca2+ダイナミクスモデル
活動電位が T 管内に達すると、筋小胞体の Ca2+ 放出チャネルが開放し、Ca2+が細胞体に放出される。細胞体内の Ca2+濃度［Ca］は主に、筋小胞体からの Ca2+の放出と吸収により制御される。ここでは、Ca2+のダイナミクスモデルとして次式のモデル［6］を採用した。
d[Ca] [Ca]in  [Ca]up
dt

(6.3)

ここで、[Ca] および [Ca] はそれぞれ次式で定義される筋小胞体からの Ca2+の放出速度および吸inup

収速度である。
[Ca] 
[Ca]in  R0 1 
 (t )*  [Ca ] and

 [Ca ] 
[Ca ]up  Vmax  [Ca ]  K m 

60

(6.4)

ここで、P0 は 1 つの活動電位の到達χ(t）による[Xα]の増加量である。Vmax および Km は酵素
反応理論における Michaelis-Menten パラメータである。

iii. 収縮タンパクのクロスブリッジダイナミクスモデル
筋小胞体から放出された Ca2+はアクチンフィラメント上のトロポニン C に結合する。その結果、アクチンフィラメントが活性化し、アクチンフィラメントとミオシンヘッド間のクロスブ
リッジの形成が促され、ミオシンヘッドが首振り運動することにより収縮力が発生する。このようなタンパク質分子の運動は熱ゆらぎのもと、確率的に進行する。
本研究では、アクチンフィラメント上の各トロポニン／トロポミオシン（T/T）ユニットおよびミオシンフィラメント上の各ミオシンヘッドに対して、それぞれ状態遷移モデルを定義し（図2.33）、それぞれの状態遷移を直接的にモンテカルロ法でシミュレートすることでこれら収縮タ
ンパクの挙動を再現し、収縮力を算出する手法［7］を採用した。
筋収縮力 Sact はクロスブリッジダイナミクスモデルより、次式で定義される[7]。

S act 

Cf



ff

Cf 

with

k xb 1 ns nHM   LNH (im, is )Asr ns is 1 im1

(6.5)(a) T/T の 3 状態モデル

No binding

Weaklybinding

StrongbindingMHrotation(b) ミオシンヘッドの 4 状態モデル
12
図 2.33: トロポニン/トロポミオシンユニット（T/T）とミオシンヘッドの遷移モデル． kon[Ca ] , kon[Ca ] ,
1, koff2 は T/T ユニットの遷移速度, knpγn, kpnγ−n, gxb , fapp, gapp, hf , hb はミオシンヘッド
koffの遷移速度を表す
ここで、  は 3 次元連続体モデルで計算される筋線維のストレッチ、f は初期配置における筋線維方向単位ベクトルである。Asr は収縮力発生ユニットの断面積、kxb はミオシンアームの弾性係数である。ns はサルコメアのサンプル数、nMH はハーフサルコメアモデルにおけるミ61オシンヘッドの数である。 LNH は時間増分中のミオシンアームの平均的な伸びである。最終的に、筋収縮力 Sact を (1)式に代入することで、1 次元筋線維モデルと 3 次元連続体モデルを連成する。
C) 骨格筋筋線維のマルチスケールシミュレーション
筋線維の活動電位モデルのシミュレーション例を図 2.34 に示す。本解析では、時刻 1.0 から 1.5 msまで、筋線維中央部に入力電流を与えた。図 2.34(a) は筋鞘および T 管に生じた膜電位であり、典型的な活動電位を再現している。図 2.34(b)は電流注入後、0、0.5 および 1.0 ms 後の膜電位の分布であり、活動電位が筋線維長手方向に沿って進行する様子が再現されている。活動電位の伝播速度は 2.09
m/s であり、実験で観察される骨格筋の活動電位とほぼ同程度の伝播速度である。

(a) サルコメアと T 管の活動電位(b) 筋繊維方向の活動電位の伝播
図 2.34: 活動電位の変化
図 2.35 は筋小胞体の Ca2+ ダイナミクスモデルの計算例で、5、10 および 15 ms 間隔の異なる周波数で活動電位が到達した時の Ca2+濃度［Ca］変化である。活動電位の到達にしたがって［Ca］は徐々に増加し、高周波数の場合には、実験で観察されているような、より高い［Ca］が再現されている。

62

図 2.35: カルシウムイオンの濃度変化

図 2.36 は収縮タンパクのクロスブリッジダイナミクスモデルによる、収縮力の計算結果である。収縮タンパクの確率的な挙動から再現された収縮力は、骨格筋の筋線維レベルの実験から得られる収縮力をほぼ同程度で再現している。

Active stress [MPa]

0.15

0.10

0.05

0.00
0.0

0.5

1.0

1.5

2.0

2.5

Time [sec]図 2.36: 一定のカルシウムイオン濃度（1.0 μM）下において, 収縮タンパクのクロスブリッジダイナミクスモデ

ルから計算される収縮力

図 2.37 は 3 次元連続体モデルと 1 次元筋線維モデルの連成シミュレーション例である。収縮タンパ
クのクロスブリッジダイナミクスモデルから計算された収縮力により、3 次元連続体モデルに収縮変形が生じている。

63

図 2.37: 収縮タンパクのクロスブリッジダイナミクスモデルから計算された収縮力により生じる筋組織の変形
D) まとめ本研究では、ヒトの運動制御メカニズムおよび神経疾患に起因する運動障害の発現メカニズムの解明のため、力学的および生理学的に精緻な詳細筋骨格系マクロ−ミクロ階層統合モデルの開発を行った。計算モデルは骨格筋の構造的な特徴を考慮して、3 次元連続体モデル内に複数の 1 次元筋線維モデルを埋め込んだ階層構造モデルを採用した。また筋収縮モデルとして、筋線維における活動電位の発生と伝播を記述する反応拡散モデル、筋小胞体におけるカルシウムイオンの流入出に関するダイナ
ミクスモデルおよび収縮タンパクのクロスブリッジダイナミクスのモデルを連成したモデルを開発した。
今後、本計算モデルをヒト下腿の神経－筋骨格モデルに適用し、大規模並列化のもと、ヒト身体運動の神経制御メカニズムの解明を行う予定である。
参考文献[1] Yamamura et al., Supercomputational Life Science 2015.
[2] Röhrle et al, SIAM J. Sci. Comput., 30, pp. 2882–2904, 2008.
[3] Yamamura et al., J. Biomech. Sci. Eng., 9, 2014.
[4] Heidlauf et al., Comput. Math. Methods Med., Volume 2013.
[5] Shorten et al., J. Muscle Res. Cell Motil., 28, pp 293–313, 2007.
[6] Zahalak et al., J. Biomech. Eng., 112, pp. 52-62, 1990.
[7] Washio et al., Multiscale Model. Simul., 11, pp. 965–999, 2013.

64

(7) 日本語摩擦音の単純形状および実形状モデルに対する空力音響シミュレーション（阪大

玉川 G,

和田 G）
A) 研究の背景と目的
歯科治療において、口腔形状の異常により発生する発音(構音)障害は QOL の向上を考えるうえで大きな問題となっており、発音障害を克服する構音トレーニングを確立する必要があるとともに、歯科治療によって発生する発音障害を治療前にあらかじめ予測・予防することが求められている。発音の中でも摩擦音は、口腔内に小さな狭めを形成することにより発生するジェット流から音が発生するため、歯科治療による口腔形状の変化に大きな影響を受ける。また、日本語のサ行で用いられる摩擦音/s/と/ʃ/は先端の数ミリの形状の違いにより発音し分けており、上手く区別して発音できない症例が多数報告されている。しかし、摩擦音はジェット流による音源発生と口腔内での音の共鳴という複雑な物理現象が組み合わさって音が形成されるため、口腔形状の違いによる音の変化を歯科医師の判断のみで予測することは困難である。また、言語聴覚士はトライ＆エラーの訓練と経験から、発音障害の患者に対して構音トレーニングを行っており、口腔形状から正確に音を予測できるようになれば、より効率的・効果的に発音障害の治療が行えると考えられる。そこで本研究では、歯茎摩擦音の口腔モ
デルを用いて数値流体シミュレーションを行うことにより、発音時の口腔内に発生する流れ場および音場を予測し、実験による計測結果と比較検証した。また、被験者が発音している際の医療画像から構築した実形状モデルに対して、口腔形状の特徴点を抽出して単純化した単純形状モデルを作製し、比較検討を行った。
B) 方法
本研究で作製した口腔実形状および単純形状モデルを図 2.38 に示す。実形状は発音障害のない被験者が/s/を発音する際の CT 画像より口腔表面を抽出し、計算格子を構築した。また、実形状の 5 点（咽頭、狭窄流路、歯茎部、上下前歯、口唇）の断面積と鉛直方向の高さより矩形流路を構築した。単純モデル出口には顔の面を摸したバッフルを設置した。また、単純モデル内の舌を摸擬した部分を前後に移動させることにより、/s/と/ʃ/の口腔形状の違いを表現した。
上記のモデルを用いて、モデルから発生する音および流れを数値流体解析により予測し、実験による計測結果と比較した。数値流体解析では、3 次元圧縮性 Navier-Stokes 方程式を支配方程式とし、有限体積法により離散的に計算した。モデル内に発生する乱流を解析するため、圧縮性流体に有効なラージエディ―シミュレーションの一方程式モデル(4)を用いた。入口は一定流速、出口は無反射境界条件を与え、壁面は断熱剛壁とした。各形状モデルの計算格子数は約 4 千万であり、40 ノードを用いて 320CPU による並列計算を行った。得られたモデル内の圧力振動に関して、離散フーリエ変換により周波数成分を計算した。モデルへ流入する空気の流量は実験および数値解析ともに 300 cm3/s とした。
実験では、モデルをアクリル板により作製し、コンプレッサ(YC-4RS, Yaezaki)から空気を流入することにより音を発生させた。発生した音はマイクロフォン(Type 4939, Bruel & Kjaer)により遠方 30cm（地点 A）においてサンプル周波数 44.1 kHz で計測した。計測した音は離散フーリエ変換により周波数スペクトルに変換した。

65

図 2.38: 歯茎摩擦音/s/を発する口腔の(a)実形状モデルと単純形状モデル
C) 結果(a) 実験による数値解析の検証
図 2.39 に単純形状モデルにおける前歯の後流 0.5 mm における平均流速および流速変動の分布を示す。赤の点線は舌を前方に設置した/s/を摸擬した形状における流速分布であり、青の一点鎖線は舌を 7 mm 後方に位置させて/ʃ/を摸擬した形状における流速分布である。
両者の口腔形状において、
平均流速および流速変動の分布を数値解析により捉えられた。
図 2.40 に/s/と/ʃ/を摸擬した単純モデルから発生した遠方場 90 mm における音圧のスペクトルを示す。両者のスペクトルともに、実験により計測されたスペクトルの特徴となる 3 kHz から 5 kHzのピークと 10 kHz 以上の高周波数域の音圧レベルを捉えることができており、本研究の数値解析手法で摩擦音の発生を十分な精度で予測できることが分かった。

図 2.39: 上前歯から 0.5 mm 離れた位置における気流速度の平均値と二乗平均平方根

図 2.40: 口唇から 90 mm 離れた位置における音圧のスペクトル分布の計測値とシミュレーション結果

66

(b) 実形状モデルと単純形状モデルの比較
摩擦音/s/の実形状および単純形状モデルにおける正中矢状面の流速分布を図 2.41 に示す。両者の口腔形状ともに狭窄流路で流速の最大値が見られ、狭窄流路から発生したジェット流は上前歯に当たって大きく乱れた。また、上下前歯から出たジェット流は下口唇に当たり、口腔から出ていった。
この結果より、/s/を発音する際の口腔形状における流れ場の特徴を単純モデルで表現できていることが分かった。図 2.42 に上記の流れ場から発生した音圧の遠方場でのスペクトルを示す。実形状から発生した 5 kHz と 9 kHz の摩擦音の特徴的なピークを単純形状モデルで再現できていることが分かる。これらの結果より、摩擦音発音時の口腔内 5 点の断面形状より単純モデルを作製することによって、口腔内の流れ場および音場を表現できることが分かった。今後は、他の被験者の口腔形状にも同様の解析を適用することにより、被験者による発音の違いや発音障害を表現できるかどうかを調べる。

図 2.41: 実形状モデル(a)および単純形状モデル(b)を用いてシミュレーション解析で得られた歯茎摩擦音/s/発生時の声道中央矢状面における気流の速度分布

図 2.42: 実形状モデル(a)および単純形状モデル(b)を用いてシミュレーション解析で得られた音圧のスペ
クトル分布(c) 摩擦音/s/と/ʃ/の違い
単純形状モデルの舌部を/s/の位置から 7 mm 後方に位置させることにより、/ʃ/の発音を表現した。
この 2 つの形状で被験者の/s/と/ʃ/の音の違いを表現できることは、実験により確認した。図 2.43 にモデル内の流れ場を示す。/s/の形状(図 2.43(a))では狭窄流路から出たジェット流がスムーズに上下
前歯の隙間へ流れていったのに対し、/ʃ/の形状(図 2.43(b))では舌を後方に位置しているため舌と前歯の間に空間ができ、その空間で旋回流が起こって狭窄流路出口近傍で大きな乱れを発生させた。

67

図 2.44 に/s/と/ʃ/の特徴であるスペクトルのピークが発生する周波数 4.7 kHz と 3.5 kHz の圧力変動分布を示す。それぞれのピーク周波数において狭窄後流で大きな圧力変動が確認された。この狭窄
前方の口腔形状の違いから発生する圧力変動の変化が、/s/と/ʃ/の音の違いを生み出していることが明らかとなった。今後は、舌部を動かしながらシミュレーションを行うことにより、単語を発音する
際の口腔内の現象を再現し、文節内で上手く発音できる場合とできない場合の原因究明に取り組む。

図 2.43: 単純形状モデルの中央矢状面における速度分布。(a) LCA = 0 mm (/s/)、および (b) LCA = 7 mm (/ʃ/)

図 2.44: 単純形状モデルの中央矢状面における圧力分布。(a) LCA = 0 mm at 4.7 kHz (/s/)、および (b) LCA =7 mm at 3.5 kHz (/ʃ/)(8) 歯科用人工知能の研究開発（阪大 玉川 G）
ビックデータや AI といった現在経済界等を中心に第四次産業の育成が計画的に実行されている。一方
で、法則性の発見に関する注目はそれらに比べると低いと言わざるを得ない。AI に分類され主軸とされている深層学習に関する技術は、機械学習による対象の認識技術のブレークスルーによるものである。
肺がん患者の胸部レントゲン画像に異常所見があるかないかの判別精度が格段に向上したという成果が典型的である。一方で、病医院等の臨床現場では、そういった人間が見れば分かるような技術も重要ではあるが、さらに人間ができないような物事を見つけ出す能力を AI に期待している。技術的には推論技術の向上が期待されている。我々は、来るべき医療 AI 技術の進歩を予測し、電子カルテデータを利用した推論技術の研究開発を行ってきた。本項ではその詳細について解説する。
A) 電子カルテ診療録テキストの処理、分析技術の研究開発
診療録テキストから知識を抽出するためには前処理が必要である。特に多くの既存推論システムを利用するためには、文章を正規化し汎化しなければならない。そのために、論理型言語 Prolog で用いられている述語意論理に対応した述語項形式で文意を表現する方法を用いた。特に、文章の意味を正確に捉えるために、文節間の関係（事態間関係）の抽出を行った。
B) 事態間知識を用いた医療推論システムの研究開発
数十万文書からなる発話に関する診療録情報から事態間情報を抽出し、バスケット解析を行った。
バスケット解析は、“一緒によく買う商品”を抽出する手法として用いられているが、本研究では、多くの一般的な述語項の中から、”発音”に特異的な項を見つけ出すため、”ある述語と一緒にまれに使われる述語”を抽出する方法を用いた。その結果、“発音”に特徴的な述語項の事態間関係が抽出できることが分かった。
68

“発音”に関する診療録テキストから抽出した事態間関係を有向グラフとして可視化（G1）すると、多くは治療方法や疾患に特異的な内容となり、“発音”に関して専門的な背景知識を持った者であればグラフの理解が可能であるが、それ以外のものに理解は難しいものとなった。そこで、専門知識も同様の方法でグラフ化（G2）し、G1 と G2 の関係を共有する述語や項で接続する取り組みを行った。
その結果、特段専門知識が無くともグラフが示す内容が理解できることが分かった。
C) 今後の予定述語項の関係ネットワークを用いて、”発音”に関する推論を行い、これまで音声言語学、音響学的に明らかにされていない事象を仮説として獲得する。

図 2.45:“発音”に関する診療録テキストから抽出した事態間関係と“発音”の専門知識から抽出した事態間関係とを結合したネットワーク

69(9) 用語集
4.

4-2.

2.

離散化

(1)
符号付距離関数ステンシル

二光子顕微鏡
4.

4-2.

2.

Neumann 条件

(2)
4.
2.
(5)

4-2.

近似計算を容易にするため、音声や電波のような連続値データを離散化する。
領域内を内外境界線からの符号付き距離で定義する関数。複雑な形状は表現が困難。
情報の範囲・自由度。高次精度の離散化を構築する方法として,ステンシル（要素、セル）を広げることがある。
エネルギーの低い光を 2 回当てても、強いエネルギーの光で観察したのと同じように見える、生物にやさしい顕微鏡。
微分方程式における解の微分がとるべき境界条件。境界面で変数の勾配を与える境界条件。
血液のような粘性のある懸濁液が小口径の管を流れるときに、見Fahreus-Lindqvist 効果

かけの粘性が減少する現象をいう。直径が 0.3mm 以下の細管で顕著に現れる。

70

3.

心臓シミュレーションと分子シミュレーションの融合による基礎医学と臨床医学の架橋（ サ ブ 課 題 C）
以下では目標 1、目標 2 に大別して実施内容を説明する。
目標 1 は、心臓シミュレータ UT-Heart と分子シミュレータ CafeMol を融合させることにより、ミクロ・マクロ間の相互作用により病態が進行する心不全の解明と最適治療を可能とするマルチスケール心臓シミュレーションを実現することである。これまでの UT-Heart シミュレータにおいては、サルコメ
アモデルを構成するミオシン分子は言わば 1 自由度を持つヘッドおよびばねで表されるアーム部で近似され、ヌクレオチド状態の変化を考慮した確率的首振り運動のダイナミクスがモンテカルロシミュレー
ションにより表現されている。このサルコメアモデルにより生理学的に妥当な心拍動、血行動態が再現されただけでなく、協調性の役割などにつき重要な知見が得られた。しかしアクチン-ミオシン収縮タン
パクを構成するアミノ酸残基の変異に起因すると言われる肥大型心筋症のような心疾患のメカニズム解
明に切り込むには限界があることも又事実である。そこで図 3.1 に示されるように、ミオシンやアクチン分子を構成するアミノ酸一つ一つの自由度とタンパクの立体構造を考慮に入れた Cα-粗視化分子モデル
をモータータンパクとする CafeMol を用いた粗視化分子動力学シミュレーションと UT-Heart による心臓シミュレーションを連成させることにより、従来と次元を異にする心疾患の解明に挑戦する。

図 3.1:「京」で用いて来たサルコメアモデルとポスト「京」で開発を目指す新たな分子モデル
CafeMol による分子動力学(MD)は構造ベースのシミュレーションである。そのため、参照構造となる
ヒトのアクチン、ミオシンの立体構造情報が必須である。しかし、現時点における PDB(Protein Data
Bank)等のデータベースにおいては、特にミオシンタンパクの各ヌクレオチド状態(ATP, ADP, ADP*pi,
Rigor)における立体構造データが不足している。従って、ここではヒトと配列相同性の高い生物種のミオシン構造を鋳型にしたホモロジーモデリングによりヒトのアクトミオシンモデルを構築した（図 3.2、図3.3、表 3.1）。配列相同性が十分高くないものもあるが、今後新たなデータが入手可能になった時点で更新する予定である。

71

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

図 3.3: モデラーにより構築された 300 サンプル中ベストのミオシン構造（Rigor 状態）
表 3.1: ホモロジーモデリングの結果（上:アクチン,下:ミオシン）72このように構築したモデル構造を用い CafeMol によって図 3.4 に示すような ATP 加水分解反応に伴うミオシンの首振り運動が実現可能か否か、フィージビリティスタディを行った。

図 3.4: ヌクレオチド状態の変化に伴うミオシンパワーストローク
その結果、モデル構造に基づく粗視化 MD でも首振り運動の再現は可能であり、首振りの際にミオシンが発生する力のオーダーも in-vitro の 1 分子計測で測定される大きさと整合性があることが分かった（図3.5、図 3.6）。

図 3.5: 分子シミュレーション結果（図 3.4 の各ヌクレオチド状態に対応）

図 3.6: 分子シミュレーションから計測される各方向の発生力と実験結果
分子シミュレーションと心臓シミュレーションの連成法については、今年度も UT-Heart チームと京大・高田教授、理研・金田研究員と議論を継続した。その結果、自由エネルギー地形から各ヌクレオチ
ド状態を離散化した確率モデルに置き換えることなく、直接的な連成を試みることとした。しかし、そのためには分子シミュレーションの時間ステップを大きく取る必要があり、CafeMol を構成する各ポテンシャル項の何れが時間ステップを律速しているかの調査、また現象の本質を損なうことなく、いかにそれを緩和し得るかの手法の検討を行った。この検討は次年度も引き続き行い計算時間短縮を目指す。
一方、心不全治療と外科手術の急性効果の検討も実施した。前年度の検討に基づいて、自治医科大学

73さいたま医療センターから成人心不全症例 4 例について心臓再同期療法施行前のデータの提供を受け、解析を進めた。この間、心臓再同期療法のシミュレーションプロトコールを整備し、調整パラメータの数を臨床データから計算可能なものに限定することで急性期の治療効果予測の精度向上を図った。小児
先天性心疾患については岡山大学からの 4 症例（両大血管右室起始症、三尖弁閉鎖症、左室低形成 2 例）について術前モデルの作製を進めたが、これに関係して、僧帽弁モデルおよび三尖弁モデルの改良を進めた。これらは修正大血管転移症における体血管系右室のモデリングにおいて重要となる。年度の後半
からは自治医科大学附属病院との協力を開始し、急性期の症例検討を進めるとともに、リモデリングに関する基礎実験を開始した。
臓器モデルの高度化に関しては、協力関係にあるアリゾナ大学 Granzier 教授との討論において受けた指摘に基づいてサルコメアの構造を修正し、より現実を反映したものとした。この結果、図 3.7 に示すように心室内に不均一に病変が存在する場合には、病変部のサルコメアが過伸展され得るとの知見が得られ、リモデリングのメカニズム解明を通じて新たな治療戦略の開発に繋がる可能性が出て来た。

①① SL Force② SL Force
60
50
40
2.5
30
20
2
10
0
0 0.2 0.4 0.6 0.8 1

Time [s]Force [kPa]

3

SL [m]

60
50
40
2.5
30
20
2
10
0
0 0.2 0.4 0.6 0.8 1

Force [kPa]

3

SL [m]②

Time [s]図 3.7:健常部①と病変部②におけるサルコメア動態の比較
目標 2 は、重点課題 1（奥野恭史・重点課題責任者）の寺田 透・東大准教授との協同により、心臓シ
ミュレータ UT-Heart と心筋細胞イオンチャネルの分子シミュレーションを融合し、創薬における候補
物質の分子構造と遺伝子情報のみから不整脈のリスク評価を可能とするシステムを開発することを目的とする。具体的には、既に開発済の図 3.8 に示される UT-Heart と細胞パッチクランプ実験と組み合わせた「ハイブリッド心毒性スクリーニングシステム」の細胞薬理実験（赤破線内）の部分を分子シミュレーションに置き換えることになる。

74

図 3.8: UT-Heart と細胞パッチクランプ実験と組み合わせた「ハイブリッド心毒性スクリーニングシステ
ム」
図 3.9 はカリウムイオンチャネルの立体構造を横（細胞膜切断面）および上（細胞外）から見たものである。青色部は電位センサードメインの立体構造、それによって囲まれる薄茶色部はポアドメインの立体構造である。電位センサードメインを構成する 4 本のヘリックスのうちの一つ（S4）の正に荷電した残基が電位センサーの役割を果たし、脱分極状態では S4 が細胞外側に移動することによりチャンネルが開く。また分極状態では S4 が細胞質側に移動しチャネル閉まる。ただし現時点ではヒト心筋の電位依存性イオンチャネルについては立体構造が決定されていないため、ここではラットのカリウムチャネル
Kv1.2 の S3 の C 末端側半分と S4 をラットの Kv2.1 に置換した（図中にピンク色で表示）キメラの立体構造（PDB ID: 2R9R）を作成し、これを計算対象とする。もちろん今後新たな構造データが得られればモデルは更新する。

図 3.9: カリウムイオンチャネルの立体構造（左図は横（細胞膜切断面）からの概観,右図は上（細胞外）
からの概観）
75

分子シミュレーション開発は二つのパートに大別される。第 1 は、電位センサードメインの立体構造変化に関するシミュレーション手法の開発、第 2 は、不整脈を誘発する薬剤と心筋イオンチャネルとの間の相互作用を予測するためのポアドメインでのドッキングシミュレーション手法の開発である。
寺田らは第 1 のパートについては、Jensen らの研究（Jensen et al. Science, 336, 229–233, 2012）に基づき、先ず、立体構造変化を説明する反応座標（集団変数）を選択し、立体構造変化過程のスナップ
ショットを用いて立体構造変化の初期経路を生成する。これを、String 法を用いて最小自由エネルギー
経路に最適化したのち、Milestoning 法により、この経路に沿った立体構造変化の速度定数を求める、という方針で開発を進めている。また第 2 のパートについては、薬剤分子は狭いチャネルへ入り込んで結合しており、粗視化分子動力学シミュレーションを用いて現象を再現するのは困難であるため、第 1 のパートと同様、全原子モデルを用いたシミュレーションを行う。具体的には、検討対象の薬剤について
ドッキングシミュレーションを実施し薬剤の結合自由エネルギーを評価することにより当該チャネルへの阻害率(IC50 等)を推定する方針で開発を進めている。
以上の分子シミュレーション開発方針に基づき、心臓シミュレーションとの連成は図 3.10 に示すような方法をとることとした。この方法によれば、新たな薬剤に対してチャネル阻害率を評価するドッキン
グシミュレーションは毎回行う必要があるものの、膜電位に依存するチャネルの構造変化は確率モデルに置き換えているため、分子シミュレーションを再度行う必要はない。多数の候補薬に対する心毒性評
価を効率的に行うことが可能な現実的連成法であると言える。

図 3.10: イオンチャネル分子シミュレーションと UT-Heart の連成方法
なお、昨年度来継続中の課題である、複数のイオン電流（遅延整流性カリウム電流 IKr、IKs、ナトリ
ウム電流 INa、INaL、カルシウム電流 ICaL）の重畳抑制がどのように不整脈発生に影響を及ぼすかについての網羅的調査は今年度終了し、現在 PNAS への投稿論文のリバイズ中である。

76

用語集
連成
分子動力学 MD (Molecular
Dynamics)

同時に進行する複数の現象を解析するとき、必要な情報を相互に渡しながら進行させること。
原子間ポテンシャルを定義して、多数の原子や分子の運動を解く方法。
“ガラス（の中）で”という意味。試験管や培養器などの中でヒトや動in-vitro

物の組織を用いて、体内に似せた人工環境において薬物の反応を検出す
る試験のこと。対となるのは in vivo（イン・ビボ）、“生体内で（の）”。

77

4.

プロジェクトの総合的推進
平成 28 年度からの本格実施フェーズが円滑に進行するよう参画各機関の連携・調整に努め、業務計画を策定・審議し、また成果を積極的に公開した。プロジェクト運営のためにテレビ会議、電話会議による運営委員会を開催し、全体の進捗状況を確認しつつサブ課題参画各機関による協力を推進した。前段
階からの成果を進化させるべく目標と体制を最適化して計画を立案・策定した。外部有識者による諮問委員会を開催し重要事項を審議に掛けた。
得られた成果は学会で発表するほか、シンポジウムやワークショップで紹介した。また重点課題広報
責任者ワーキンググループにおいても成果取扱い改善を提言し、全重点課題を含むシンポジウム等においてもポスター掲示などで成果をアピールした。さらにホームページ運営では、随時成果の掲載や英語ページ立ち上げ等、コンテンツの更なる拡充を行った。またニュースレター（電子版）を発刊し、ホームページに掲載した。
ポスト「京」重点課題推進 WG、重点課題連携推進会議、HPCI 計画推進委員会、人材育成タスフォー
ス、ポスト「京」開発の概要等に関する説明会、計算科学発展のための検討会等、全重点課題に関する会議において連携・調整を行った。ソフトウェア普及においては、重点課題 2 のターゲットアプリケー
ションである Genomon2 の講習会を開催した。人材育成の一環として「生命情報と人工知能セミナー」を開催した。
(1) 運営委員会
プロジェクト運営を円滑かつ効率的に行うために、重点課題 2 責任者、サブ課題責任者などで構成した運営委員会を実施している。研究開発やプロジェクトの総合的推進の活動進捗の確認や、課題についての調整等を、サブ課題参画各機関と協調して行った。
平成 28 年度は以下の日程で開催した。
第 1 回 5 月 20 日(金)東京大学医科学研究所にて開催

第 2 回 6 月 17 日(金)東京大学医科学研究所にて開催

第 3 回 7 月 15 日(金)東京大学医科学研究所にて開催

第 4 回 8 月 19 日(金)東京大学医科学研究所にて開催

第 5 回 9 月 16 日(金)東京大学医科学研究所にて開催

第 6 回 10 月 28 日(金)東京大学医科学研究所にて開催

第 7 回 11 月 18 日(金)東京大学医科学研究所にて開催

第 8 回 12 月 16 日(金)東京大学医科学研究所にて開催

第 9 回 1 月 23 日(月)国立科学博物館講堂にて開催

第 10 回 2 月 17 日(金)東京大学医科学研究所にて開催

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

78

に反映させた。
A) 諮問委員（◎委員長）
礒田 治夫 （名古屋大学脳とこころの研究センター 教授）
梶谷 文彦 （川崎医療福祉大学 名誉教授)
高井 義美 （神戸大学大学院医学研究科 教授）◎樋口 知之 （情報・システム研究機構

統計数理研究所 所長）

B) 開催状況2016 年 12 月 19 日(月)研究開発の進捗報告,研究開発内容の審議および助言

2017 年 3 月 29 日（水）審議および助言に対する対応状況の説明,平成 28 年度成果の報告(3) シンポジウムおよびワークショップの開催
他機関と連携したシンポジウムや学会においてワークショップ等を開催した。
平成 28 年度に行った活動概要は以下の通り（共催,協力,協賛,後援を含む）。
A) シンポジウム(a) 名称：その予防・医療,時代遅れです－ポスト「京」が必要な訳－
開催日：2017 年 1 月 23 日(月)場所：国立科学博物館
参加者数：95 名
主催：ポスト「京」重点課題 2

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

協賛：「システム癌新次元」文部科学省科学研究費補助金 新学
術領域研究 (複合領域 4701) がんシステムの新次元俯瞰と攻略
東京大学医科学研究所の古川洋一教授を基調講演者として招聘し、ゲノム解析から明らかになった若年性腫瘍の背景となるゲノム変異や、治療抵抗性となった腫瘍のゲノム変異などについての講演を行った.サブ課題Ａ、Ｂ、Ｃからはサブ課題責任者が登壇し、研究開発の最新状況について発表を行った。
(b) 名称：KOBE HPC Summer School 2016
開催日：2016 年 8 月 1 日(月)～8 月 5 日（金）場所：理化学研究所計算科学研究機構 6 階講堂
主催：神戸大学計算科学教育センター、理化学研究所計算科学研究機構、他 3 機関
後援：ポスト「京」重点課題 9 課題、高度情報科学技術研究機構(c) 名称：第 3 回 大型実験施設とスーパーコンピュータとの連携利用シンポジウム
開催日：2016 年 9 月 1 日 (木)場所：東京・秋葉原 UDX 6 階カンファレンス
主催：高輝度光科学研究センター、総合科学研究機構、高度情報科学技術研究機構、
ポスト「京」重点課題 579

協賛：理化学研究所計算科学研究機構、ポスト「京」重点課題 2－4・6－7・9、他 7 機関
協力：ポスト「京」重点課題 8
後援：文部科学省、兵庫県、茨城県(d) 名称：第 3 回「京」を中核とする HPCI システム利用研究課題 成果報告会
開催日：2016 年 10 月 21 日(金)場所：東京・コクヨホール
主催：高度情報科学技術研究機構
協賛：HPCI コンソーシアム、スーパーコンピューティング技術産業応用協議会
協力：理化学研究所計算科学研究機構、国立情報学研究所、他 11 機関
後援：計算科学振興財団、ポスト「京」重点課題 9 課題(e) 名称：RIKEN AICS HPC Youth Workshop 2016
開催日：2016 年 11 月 28 日(月)～11 月 30 日(水)場所：理化学研究所計算科学研究機構 6 階講堂
主催：理化学研究所計算科学研究機構
後援：ポスト「京」重点課題 9 課題、計算科学振興財団、HPCI コンソーシアム、
高度情報科学技術研究機構(f) 名称：1st International Symposium on Research and Education of Computational Science
開催日：2016 年 11 月 29 日(火)～11 月 30 日(水)場所：Koshiba Hall, School of Science, The University of Tokyo
主催：The Computational Science Alliance of the University of Tokyo
協賛：HPCI Consortium, Joint Usage/Research Center for Interdisciplinary Large-scale
Information Infrastructures, Post-K Computer Priority Issue 1-3,5-9(g) 名称：人工知能とスーパーコンピュータでがんにチャレンジ
開催日：2016 年 12 月 25 日(日)場所：東京・お台場 日本科学未来館 7 階 未来館ホール
主催：「システム癌新次元」文部科学省科学研究費補助金 新学術領域研究 (複合領域 4701)「がんシステムの新次元俯瞰と攻略」
協賛：ポスト「京」重点課題 2(h) 名称：The 7th AICS International Symposium Emerging numerical techniques for exascale
and post-Moore era
開催日：2017 年 2 月 23 日（木）～24 日(金)場所：RIKEN Advanced Institute for Computational Science
主催：RIKEN Advanced Institute for Computational Science80

後援：Priority issues (9 issues) on Post-K computer, RIKEN Interdisciplinary Theoretical andMathematical Sciences Program, Research Organization for Information Science andTechnology
協賛：FOCUS Establishing Supercomputing Center of Excellence(i) 名称：見える化シンポジウム 2017 ～シミュレーションの価値～
開催日：2017 年 3 月 11 日(土)場所：東京・日本橋ライフサイエンスハブ
主催：理化学研究所計算科学研究機構
協力：ポスト「京」重点課題 9 課題、高度情報科学技術研究機構、豊橋技術科学大学(j) 名称：KOBE HPC Spring School 2017
開催日：2017 年 3 月 13 日（月）～ 3 月 15 日（水）場所：神戸大学計算科学教育センターセミナー室 208
主催：神戸大学計算科学教育センター、理化学研究所計算科学研究機構、東京大学情報基盤セン
ター、兵庫県立大学大学院シミュレーション学研究科、計算科学振興財団
後援：ポスト「京」重点課題 9 課題、高度情報科学技術研究機構(k) 名称：International Symposium on IMSUT Joint Research Project "Deciphering Genome Big
Data"
開催日：2017 年 3 月 21 日(火)場所：Tomy Conference Hall in IMSUT Hospital Building
主催：Institute of Medical Science, the University of Tokyo (IMSUT)
共催：Medical Genomics Research Initiative, The University of Tokyo
協賛：Priority Issue on Post-K computer (Integrated Computational Life Science to Support
Personalized and Preventive Medicine)B) ワークショップ(a) 名称：ポスト「京」重点課題 2 ワークショップ
開催日：2016 年 12 月 19 日(月)場所：東京大学医科学研究所
参加者数：38 名
研究開発の進捗確認および報告の場、またプロジェクト参加機関内での実質的な相互連携の推進を目的として、プロジェクト参加研究者を対象にワークショップを開催した。参加研究者をはじめ諮問

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

委員からの質疑応答では活発な議論が交わされた。
(b) 名 称 ： International Workshop on Integrated Computational Life Sciences towards

81Personalized and Preventive Medicine
開催日：2016 年 11 月 1 日(火)参加者数：約 20 名
発表タイトルと発表者（発表者は Bold 文字で表記）：
・Comparative analysis of transformation methods for gene expression profi les in breast
cancer datasets (Yoshiaki Sota, Shigeto Seno, Yoichi Takenaka, Shinzaburo Noguchi and
Hideo Matsuda)
・ Patient-specific blood flows simulation on cerebral aneurysm based on physicallyconsistency feedback control (Mohd Azrul Hisham Mohd Adib, Satoshi Ii, Yoshiyuki
Watanabe and Shigeo Wada)・Development of Multi-scale Musculo-skeletal Simulator (Naoto Yamamura, Shu Takagi andTaishin Nomura)・ A Computational Approach for Blood Flow Analysis in the Densely Coiled Cerebral
Aneurysm (Tomohiro Otani, Satoshi Ii, Tomoyoshi Shigematsu, Toshiyuki Fujinaka,
Masayuki, Hirata, Tomohiko Ozaki and Shigeo Wada)
オーガナイザー：松田 秀雄（ポスト「京」重点課題 2 サブ課題Ｂ責任者）
学 会 名 称 ： The 16th annual IEEE International Conference on Bioinformatics andBioengineering
開催期間：2016 年 10 月 31 日(月)～11 月 2 日(水)場所：スプレンダーホテル台中、台中市、台湾(4) ポスト「京」重点課題広報責任者ワーキンググループ
ポスト「京」の研究成果やシミュレーションの重要性・必要性を、重点課題全体として効果的に行うこと、また広報に関する重点課題間の連携推進等を目的として、広報責任者ワーキンググループが開催された。理化学研究所計算科学研究機構・重点課題 9 課題・高度情報科学技術研究機構の広報責任者および文部科学省で構成され、重点課題全体で実践する成果の普及活動等、計算科学の発展に資する議論を行った。
平成 28 年度に開催したワーキンググループは以下の通り.
第 1 回 4 月 28 日(木)理研東京連絡事務所会議室,理研 AICS（TV 会議）にて開催

第 2 回 7 月 7 日(木)理研東京連絡事務所会議室,理研 AICS（TV 会議）にて開催

第 3 回 8 月 4 日(木)理研東京連絡事務所会議室,理研 AICS（TV 会議）にて開催

第 4 回 9 月 29 日(木)理研東京連絡事務所会議室,理研 AICS（TV 会議）にて開催

第 5 回 12 月 14 日(水) 理研東京連絡事務所会議室,理研 AICS（TV 会議）にて開催
第 6 回 1 月 26 日(木)理研東京連絡事務所会議室,理研 AICS（TV 会議）にて開催(5) ホームページ（日本語・英語）
本年度からのプロジェクト本格稼働に伴い、サイト構成やコンテンツの更なる拡充を行った。またホ82

ームページの英語版を開設した。

図 4.3:日本語版

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

図 4.5: Vol.1（9 月）

図 4.6: Vol.2（2 月）

図 4.7: パンフレット（12 月）(7) ロゴマーク制作
全サブ課題の研究内容に共通する「ヒト」をモチーフにし、ポスト「京」重点課題 2 の活動を象徴したロゴマークを作成した。

図 4.8: ロゴマーク（日本語版と英語版）

83

(8) セミナー,講習会等の開催
A) 人材育成セミナー
情報生命科学の最先端で活躍する研究者を招聘し、将来の研究に役立つ実践的な知識を得ることを目指した一般公開セミナーを開催し、産業技術総合研究所人工知能研究センター（ＡＩＲＣ）と東京大学大学院新領域創成科学研究科メディカル情報生命専攻と連携し、最先端の研究紹介を行った。このセミナーは東京大学大学院 新領域創成科学研究科「情報生命科学特別講義 Ⅲ」を兼ねており、各キャンパスに TV 会議システムで配信した。※第 2 回（2016 年 10 月 14 日）以降が一般公開
日程：2016 年 10 月 7 日～2017 年 1 月 27 日（14:55～16:40）毎週金曜日 全 12 回
開催場所：産総研臨海副都心センター別館 10F 人工知能研究センター会議室 2 ほか 3 キャンパス）表 4.1「生命情報と人工知能セミナー」（情報生命科学特別講座Ⅲ）セミナーの詳細
10 月 7 日

第1回

新しい生命科学の方法論を目指して
辻井潤一（産業技術総合研究所人工知能研究センター）

10 月 14 日

第2回

創薬インフォマティクスによる分子設計
広川貴次（産業技術総合研究所創薬分子プロファイリング研究センター）

10 月 28 日

第3回

機械学習とデータ駆動型ドラッグリポジショニングによる創薬
山西芳裕（九州大学生体防御医学研究所）

11 月 4 日

第4回

確率潜在意味解析とベイジアンネットワークによる人と相互理解する次世代人工知能技術
本村陽一（産業技術総合研究所人工知能研究センター）

11 月 11 日

第5回

次世代シークエンス技術がかえる生命情報解析鈴木穣（東京大学新領域創成科学研究科）

11 月 18 日

第6回

データベース統合
高木利久（東京大学理学系研究科）

11 月 25 日

第7回

人工知能・シミュレーション技術の創薬応用
奥野恭史（京都大学医学系研究科）

12 月 9 日

第8回

Finding homologies and rearrangements in huge bio-sequence datasets
Martin Frith（産業技術総合研究所人工知能研究センター・東京大学新領域創成科学研究科）

12 月 16 日

第9回

創薬等支援に向けたタンパク質基質結合部位の大規模比較解析
富井健太郎（産業技術総合研究所人工知能研究センター）

1 月 13 日

第 10 回

がんを理解するための大規模ゲノムデータ解析と人工知能応用
宮野悟（東京大学医科学研究所）

1 月 20 日

第 11 回長鎖ノンコーディング RNA の機能の解明に向けたバイオインフォマティクス技術
浜田道昭（早稲田大学理工学術院）

1 月 27 日

第 12 回大量・複雑なデータの理解に向けた新たな統計手法
津田宏冶（東京大学新領域創成科学研究科）

84

B) 講習会
主にがんのゲノム・トランスクリプトーム解析に興味のある実験研究者、医学研究者を対象に、がんの後天的変異・構造変異・融合遺伝子を検出する一般的原理を解説し、このような解析におけるGenomon2 の特長について講習会を行った。また、第一線でがんゲノム解析研究に従事する医学研究者による最先端のがんゲノム研究を紹介した。
開催場所：東京大学医科学研究所
表 4.2 講習会の詳細
5 月 24 日

がんゲノムシークエンス解析の原理と Genomon2 の紹介

Genomon2 Tutorial

Genomon2 を利用したがんゲノム解析の実際

1:00-4:00pm

ハンズオンセミナー
6 月 24 日

Genomon2 Tutorial

がんゲノムシークエンス解析の原理と Genomon2 の紹介

1:00-5:00pm

実践編

Genomon2 ハンズオンセッション

4－3.活動（研究会等）
平成 28 年度に実施した実績についての概要は、別添の様式第 21 に記載.
様式第 21 における学会誌・雑誌等における論文掲載と学会等における口頭・ポスター発表以外の活動について以下に示す。
(1) 受賞等№1

名称
高松宮妃癌研究基金
学術賞

2

上原賞

3

中谷賞

4

受章者氏名

授賞機関(学会名等)
公益財団法人小川 誠司

2016 Poster Award

高松宮妃癌研2017 年 2 月 23 日

究基金
公益財団法人

小川 誠司、 宮野 悟

大賞

受賞時期

上原記念生命科学財団
公益財団法人

杉浦 清了

2016 年 12 月 20 日

中谷医工計測2017 年 2 月 24 日

技術振興財団

T. Otani, S. Ii, T.

13th interdisciplinary

Shigematsu, T.

cerebrovascular symposium

Fujinaka, T. Ozaki, M.

international intracranial

Hirata, S. Wada

stent meeting, 神戸

2016 年 11 月 26 日

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

タイトル

媒体・方法

媒体種類

時期

1

大規模データ解析で生命システムを明らかに

ポスト「京」重点課題２ニュースレターVol.2

Web 上 広報誌

2017 年 2 月 1 日

85

2

3
4
5

6

7

8

骨髄異形成症候群におけるクロー
ン進化の解明について（タイトルは各社により異なる）
白血病発症リスク高める遺伝子
京大・東大など発見
京大・東大など、白血病発症リス
ク高める遺伝子特定
骨髄異形成症候群におけるクロー
ン進化は 2 次性白血病への進展と関連する
Where and When Did it ALL Go
Wrong?
DicipheringtheMolecular Pathophysiology of
Disease Progression in MDS.
Investigating Clonal Dynamics
of Different Mutations in
Myelodysplastic Syndromes.
がん細胞が免疫から逃れるメカニ
ズムの解明− 免疫チェックポイント阻害剤のプレスリリース

TV、新聞2016 年 12 月 20日

日本経済新聞（朝刊 42 面）

新聞

2016/12/21

日本経済新聞Web （ 新聞）Web

2016/12/21

The Hematologist （ 2017 年
3/4 月号 11 ページ）

雑誌

2017/3/1

ASH Clinical News （2017 年 3月号 24/25 ページ）

雑誌

2017/3/1

プレスリリース

TV、新聞2016 年 5 月 24 日

ライフサイエンス 新着論文レビュー

2017/1/13

効果予測への応⽤に期待 −9遺伝⼦異常でがんへの免疫抑制

NHK 関⻄ニュースTV

2016 年 5 月 24 日

10

がん、免疫の攻撃回避仕組み解明

NHK 総合⼤阪

TV

2016 年 5 月 24 日

京都⼤など、英科学誌に発表

NHK ニュース おはよう⽇本

⾼額抗がん剤の効き⽬⾒分ける可

毎⽇放送

TV

2016 年 5 月 24 日

新聞

2016 年 5 月 24 日

⽇本経済新聞（朝刊 34⾯）

新聞

2016 年 5 月 24 日

毎⽇新聞（朝刊 28⾯）

新聞

2016 年 5 月 24 日

11

能性
12

京⼤グル―プ発表

⾼額抗がん剤、効き⽬予測 京⼤、 朝⽇新聞（朝刊 3⾯）

14

たんぱく質の⽬印発⾒免疫細胞の攻撃回避の仕組み解明
英科学誌に発表
がん細胞：「オプジーボ」の適性判

15

別京⼤⼤学院チーム解明
新型抗がん剤、有効な患者絞り込み
がん治療薬、効くかどうか判別 京

京都新聞（朝刊 28⾯）

新聞

2016 年 5 月 24 日

16

⼤、⼿法確⽴つながる成果
がんの免疫回避メカニズムを解明

⽇刊⼯業新聞（25⾯）

新聞2016 年 5 月 24 日

17

がんで免役⼒低下

産経経済新聞（朝刊 24⾯）

新聞

2016 年 5 月 24 日

中⽇新聞（朝刊 3⾯

新聞

2016 年 5 月 24 日

四国新聞（朝刊 3⾯）

新聞

2016 年 5 月 24 日

東奥⽇報（朝刊 20⾯）

新聞

2016 年 5 月 24 日

13

部解明
18
19

仕組みの⼀

京⼤チーム

京⼤など
遺伝⼦解析で⾼額のがん新薬効く患者特定へ
がんの免疫回避メカニズム解明
京都⼤、「遺伝⼦に異常」

20

がんの免疫回避判明

京⼤など

遺伝⼦特定領域に異常86

21

免疫の攻撃回避、がんの作⽤解明

秋⽥魁新聞（朝刊 24⾯）

新聞

2016 年 5 月 24 日

京⼤などのチーム
22

23
24
25

⾼額抗がん剤、効き⽬予測 京⼤、 朝⽇新聞、Google ニュース、 Web （ 新NEWSPICKS、CURASAW、2 チャン 聞）
たんぱく質の⽬印発⾒ネルニュース速報、Yahoo ニュース
Web （ 新がん免疫薬、効き⽬の⽬印 京⼤な ⽇本経済新聞
聞）ど発⾒2016 年 5 月 24 日

がん細胞：免疫細胞の攻撃回避の仕組み解明
免疫回避の仕組み解明＝がん治

Web （ 新聞）
Web （ 新聞）2016 年 5 月 24 日

⽇刊⼯業新聞Web （ 新聞）2016 年 5 月 24 日

がん、免疫の攻撃回避仕組み解明

東京新聞、どうしんウェブ、

2016 年 5 月 24 日

京都⼤など、英科学誌に発表

⻄⽇本新 聞、神奈川新聞カナ
コロ、共同通信、

Web
（新聞）

Web

2016 年 5 月 24 日

産経ニュース、産経 WEST、終活
WEB ソナエWeb2016 年 5 月 24 日

マイナビニュース、excite ニュース、 goo ニュース、ニュースWeb

2016 年 5 月 24 日

Web2016 年 5 月 24 日

毎⽇新聞
時事通信

療、効果的に−京⼤など
26

がんの免疫回避メカニズムを解明−京⼤27

2016 年 5 月 24 日

2016 年 5 月 24 日

⾼知新聞、秋⽥魁新報社、
47NEWS、 静岡新聞、河北新報、
⻑崎新聞、Web 東奥、下野新聞、
Newspect、⼭陽新 聞、⼭陰中
央新報、デーリー東北新聞 社、
四国新聞社、徳島新聞、佐賀新
聞、上⽑新聞、福島⺠報、⼤分
合同新 聞、⼭形新聞、信濃新聞、千葉⽇報、 中国新聞、中⽇新聞、愛媛新聞、⽇本 海新聞、茨城新聞、神⼾新聞、北⽇
本 新聞、スマートマイズニュ
28

遺伝⼦異常でがんへの免疫抑制 -ース、信濃 毎⽇新聞
NHK オンライン

NHK 関⻄ NEWS WEB 2930

がん細胞遺伝⼦に特定の異常→免疫の攻撃⼒低下 京⼤チーム、メカ
ニズ
ム解明 スパコン「京」も駆使
がんの免疫チェックポイント阻害薬の効果を予測するマーカー 京⼤ら 発⾒31＜がん細胞＞免疫細胞の攻撃回避の仕組み解明

記録アンテ ナ、⽇本統合医療推奨協会
Biglobe ニュース、You tube

87

32

免疫回避の仕組み解明＝がん治33

療、効果的に−京⼤など
がん細胞の免疫監視回避のメカニ

時事ドットコムニュース、アメーバニュースWeb

2016 年 5 月 24 日

QlifePro 医療 NEWSWeb

2016 年 5 月 27 日

JCAST ヘルスケア

Web2016 年 5 月 29 日

ズム解明−北⼤ほか
34

がん遺伝⼦の異常で免疫の攻撃を回避

京⼤などチームが発⾒

35

SiGN-BN Para-OS

プログラム公開Web

2017 年 2 月 20 日

1

生体力学シミュレーションの可視
化事例

Web

2017 年 1 月 1 日

2

歯茎摩擦音/s/の空力音響シミュ
レーションの可視化事例

Web2016 年 11 月 1 日

3歯茎摩擦音/s/の空力音響シミュ
レーションと実験による検証

大阪大学サイバーメディアセンターHP(http://vis.cmc.osaka-u.ac.jp/cat_gallery/usercase03/)
大阪大学サイバーメディアセンターHP(http://vis.cmc.osaka-u.ac.jp/cat_gallery/usercase02/)
大阪大学サイバーメディアセンターCYBERMEDIA HPC JOURNAL
利用者研究報告

広報誌

2016 年 8 月 1 日

1

WEB 雑誌 日経バイオテク

Web 上 広報誌

2017 年 3 月 23 日

2Bridging
ContinuumandMolecular Models in a Beating
Heart Simulation

富士フイルムと東大発ベンチ
ャーの心臓シミュレーター、
2017 年夏にも治験申請 手術
シミュレーションシステムのアプリケーションとして開発へ(https://bio.nikkeibp.co.jp/atcl/news/p1/17/03/20/0247
3/)
SIAM NEWS DECEMBER 2016

Web 配 信
（動画な
し） (鷲
尾 巧)

2016 年 12 月 1 日

４－４.実施体制
１.業務主任者（受託者（委託先））
役職・氏名 国立大学法人東京大学 医科学研究所
E-メールアドレス： miyano@ims.u-tokyo.ac.jp
TEL：03-5449-5615 FAX: 03-5449-5442
88

教授 宮野 悟（再委託先）役職・氏名 国立大学法人大阪大学 大学院基礎工学研究科 教授

和田 成生

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： 0285-58-7006 FAX：0285-44-5019（再委託先）
役職・氏名 国立大学法人岡山大学 大学院医歯薬学総合研究科 教授
佐野 俊二
E-メールアドレス：s_sano@cc.okayama-u.ac.jp
TEL：086-235-7357

FAX：086-235-7431

89２.業務項目別実施区分
業

務

項

目

実

施 場 所

担

当 責 任 者

（１）大量シーケンスによ 東京都港区白金台4-6-1国 ヒトゲノム解析センター長
るがんの個性と時立大学法人東京大学 医間的・空間的多様科学研究所ヒトゲノム解

宮野 悟

性・起源の解明（サ 析センター
ブ課題 A）
京都府京都市左京区吉田大学院医学研究科 教授

近衛町

小川 誠司

国立大学法人京都大学大学院医学研究科（２）データ同化生体シミ 大阪府豊中市待兼山町1-3 大学院基礎工学研究科 教授ュレーションによ国立大学法人大阪大学

る個別化医療支援

大学院基礎工学研究科

和田 成生（サブ課題B）
東京都文京区本郷 7-3-1

大学院工学系研究科 教授

東京大学大学院工学系研

高木 周

究科
大阪府吹田市山田丘1-1

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

国立大学法人大阪大学大松田 秀雄

学院情報科学研究科
大阪府豊中市待兼山町1-3 大学院基礎工学研究科 教授
国立大学法人大阪大学

野村 泰伸

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

歯学部附属病院 准教授

国立大学法人大阪大学歯玉川 裕夫

学部附属病院（３）心臓シミュレーショ 東京都世田谷区野沢3－25 代表取締役会長
ンと分子シミュレ－8

久田俊明

ーションの融合に

株式会社UT-Heart研究所

よる基礎医学と臨
床医学の架橋栃木県下野市薬師寺（サブ課題C）

3311-190

学長 永井 良三

学校法人自治医科大学
岡山県岡山市北区鹿田町

大学院医歯薬学総合研究科

2-5-1

教授 佐野 俊二

国立大学法人岡山大学大学院医歯薬学総合研究科（４）プロジェクトの総合 東京都港区白金台4-6-1国 ヒトゲノム解析センター長
的推進

立大学法人東京大学 医科学研究所ヒトゲノム解
析センター

91

宮野 悟

３.経理担当者（受託者（委託先）)
役職・氏名 国立大学法人東京大学 医科学研究所

事務部 研究支援課

外部資金戦略チーム 主任 齋藤 利弥
E-メールアドレス： t-gshikin@ims.u-tokyo.ac.jp
TEL：03-5449-5138

FAX： 03-6409-2017（再委託先）
役職・氏名 国立大学法人大阪大学 大学院基礎工学研究科
研究協力係 係長 武中 典子
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（再委託先）
役職・氏名 国立大学法人岡山大学 大学院医歯薬学総合研究科等事務部
総務課 主査 岩﨑 哲也
E-メールアドレス：iwasaki-t@adm.okayama-u.ac.jp
TEL：086-235-7014

FAX：086-235-7046

４.知的財産権の帰属「知的財産権は乙に帰属することを希望する.」５.委託契約書の定めにより甲に提出することとされている著作物以外で委託業務により作成し,甲に納入する著作物の有無（ 有 ・ 無 ）92