特許 5946127 シミュレーション方法、シミュレーション装置、およびシミュレーションプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】5946127

(24)【登録日】2016年6月10日

(45)【発行日】2016年7月5日

(54)【発明の名称】シミュレーション方法、シミュレーション装置、およびシミュレーションプログラム

(51)【国際特許分類】

   A61B 5/00 20060101AFI20160621BHJP

   A61B 6/00 20060101ALI20160621BHJP

   A61B 6/03 20060101ALI20160621BHJP

【ＦＩ】

   A61B5/00 G

   A61B6/00 360Z

   A61B6/03 360Z

【請求項の数】9

【全頁数】35

(21)【出願番号】特願2012-109067(P2012-109067)

(22)【出願日】2012年5月11日

(65)【公開番号】特開2013-233369(P2013-233369A)

(43)【公開日】2013年11月21日

【審査請求日】2015年2月4日

(73)【特許権者】

【識別番号】000005223

【氏名又は名称】富士通株式会社

(73)【特許権者】

【識別番号】504137912

【氏名又は名称】国立大学法人 東京大学

(74)【代理人】

【識別番号】100092152

【弁理士】

【氏名又は名称】服部 毅巖

(72)【発明者】

【氏名】岩村 尚

(72)【発明者】

【氏名】久田 俊明

(72)【発明者】

【氏名】杉浦 清了

(72)【発明者】

【氏名】鷲尾 巧

(72)【発明者】

【氏名】岡田 純一

【審査官】 伊知地 和之

(56)【参考文献】

【文献】 特表２０１３−５３４１５４（ＪＰ，Ａ）

【文献】 国際公開第２００５／０８３６１５（ＷＯ，Ａ１）

【文献】 特表２００３−５０８１５２（ＪＰ，Ａ）

【文献】 米国特許出願公開第２０１２／０１２１１５１（ＵＳ，Ａ１）

【文献】 特表２０１２−５１９９０２（ＪＰ，Ａ）

【文献】 特表２０１３−５０１２９０（ＪＰ，Ａ）

【文献】 H.J.Kim et al.，"Patient-Specific Modeling of Blood Flow and Pressure in Human Coronary Arteries"，Annuals of Biomedical Engineering，２０１０年１０月，Vol.38, No.10，pp.3195-3209

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ａ６１Ｂ ５／００ − ５／０３

Ａ６１Ｂ ６／００ − ６／１４

Ｇ０６Ｔ １／００

Ｇ０６Ｔ ７／００ − ７／６０

ＣＳＤＢ（日本国特許庁）

(57)【特許請求の範囲】

【請求項1】

コンピュータが、
複数の節点間を接続して得られる複数の要素により器官の形状を定義した形状モデルと、前記器官内の閾値以上の径の血管による血管ネットワークを定義した第１の血管ネットワークモデルとに基づいて、前記器官の動きに応じた前記第１の血管ネットワークモデルの血管内の血液の状態を解析し、
前記器官内の前記閾値未満の径の血管による血管ネットワークを定義した第２の血管ネットワークモデルが、前記複数の節点それぞれに、１つまたは複数接続されており、前記第１の血管ネットワークモデルの血管内の血液の状態の解析で得られた、前記複数の節点それぞれにおける血液の状態を示す情報を用いて、前記複数の節点にそれぞれ接続された第２の血管ネットワークモデルの血液の状態を、解析することを特徴とするシミュレーション方法。

【請求項2】

前記第２の血管ネットワークモデルは、毛細血管、および毛細血管と前記第１の血管ネットワークモデルに定義された血管ネットワークとを接続する中間層の血管で構成され、該中間層の血管ネットワークは、毛細血管を挟んで、動脈側と静脈側とが対称となるように定義されていることを特徴とする請求項１記載のシミュレーション方法。

【請求項3】

前記器官の解剖学的特徴として、節点の位置の血管密度が高いほど、該節点に接続される第２の血管ネットワークモデルの数を多くすることを特徴とする請求項１または２に記載のシミュレーション方法。

【請求項4】

節点に接続された第２の血管ネットワークモデルの血液の状態の解析では、前記第１の血管ネットワークモデルの血管内の血液の状態の解析で得られた、該節点における圧力に関する情報を用いて解析を行うことを特徴とする請求項１乃至３のいずれかに記載のシミュレーション方法。

【請求項5】

前記器官は心臓であり、前記器官の動きは、心筋の収縮・弛緩による動きであることを特徴とする請求項１乃至４のいずれかに記載のシミュレーション方法。

【請求項6】

心外膜側の節点に接続される第２の血管ネットワークモデルの数よりも、心内膜側の節点に接続される第２の血管ネットワークモデルの数が多いことを特徴とする請求項５記載のシミュレーション方法。

【請求項7】

前記コンピュータは、
複数の演算処理装置を有し、
前記シミュレーション方法は、
前記複数の節点にそれぞれ接続された第２の血管ネットワークモデルの血液の状態を、前記複数の演算処理装置を用いてそれぞれ解析することを特徴とする請求項１乃至６のいずれかに記載のシミュレーション方法。

【請求項8】

複数の節点間を接続して得られる複数の要素により器官の形状を定義した形状モデルと、前記器官内の閾値以上の径の血管による血管ネットワークを定義した第１の血管ネットワークモデルとに基づいて、前記器官の動きに応じた前記第１の血管ネットワークモデルの血管内の血液の状態を解析する第１の解析部と、
前記器官内の前記閾値未満の径の血管による血管ネットワークを定義した第２の血管ネットワークモデルが、前記複数の節点それぞれに、１つまたは複数接続されており、前記第１の血管ネットワークモデルの血管内の血液の状態の解析で得られた、前記複数の節点それぞれにおける血液の状態を示す情報を用いて、前記複数の節点それぞれに接続された第２の血管ネットワークモデルそれぞれの血液の状態を、解析する第２の解析部を有することを特徴とするシミュレーション装置。

【請求項9】

コンピュータに、
複数の節点間を接続して得られる複数の要素により器官の形状を定義した形状モデルと、前記器官内の閾値以上の径の血管による血管ネットワークを定義した第１の血管ネットワークモデルとに基づいて、前記器官の動きに応じた前記第１の血管ネットワークモデルの血管内の血液の状態を解析させ、
前記器官内の前記閾値未満の径の血管による血管ネットワークを定義した第２の血管ネットワークモデルが、前記複数の節点それぞれに、１つまたは複数接続されており、前記第１の血管ネットワークモデルの血管内の血液の状態の解析で得られた、前記複数の節点それぞれにおける血液の状態を示す情報を用いて、前記複数の節点それぞれに接続された第２の血管ネットワークモデルそれぞれの血液の状態を、複数の演算装置による並列処理で解析させることを特徴とするシミュレーションプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は血行動態のシミュレーション方法、シミュレーション装置、およびシミュレーションプログラムに関する。

【背景技術】

【0002】

近年における計測技術の発展に伴い、生物のタンパク質、細胞、器官ごとに、膨大な知見が得られている。なお、タンパク質、細胞、器官は、それぞれ計測対象のスケールが異なる。そして、計測技術では、知見はスケールごとに得られるものである。そのため計測技術だけでは、スケールが異なる計測対処物間の関連、例えば細胞と器官の挙動が互いにどのように影響するか、については知見が十分に得られていない。そこで、性能向上が著しい計算機を用いたコンピュータシミュレーションにより、スケールが異なる計測対処物間の関連に関する知見不足を補う研究が活発になっている。このような研究により、従来の測定・実験だけでは観察困難であった現象が解明されつつある。

【0003】

例えばコンピュータシミュレーションによる研究が盛んな器官の１つとして、心臓がある。心臓は全身に血液を送り出すという生命活動の根幹を担う重要な器官である。そして心臓自身もまた安定した血液の供給を必要とする。血液供給の役割を担うのは、心臓の内部を走行する血管である、冠循環である。冠循環の異常に基づく病態は、虚血性心疾患と総称される。

【0004】

虚血性心疾患は、主に冠動脈の動脈硬化によって引き起こされる。動脈硬化については、生化学、細胞生物学的側面を中心に研究が進められている。その一方で、虚血性心疾患の病態形成には、冠循環系の構造およびその制御といった、機能的な要因が大きく関わっている。例えば、冠循環系の血流は心収縮の影響を強く受けると考えられている。実際に狭心症を含めた病態では、最も組織圧が上昇すると考えられる心内膜側（心腔側）で虚血が生じやすいことが知られている。組織圧は、力学解析で求めることができる。また、虚血性心疾患の病態を把握する上では、粥腫の形成に対する血流（ずり応力）も重要な観測対象であり、このような血流も力学解析で求めることができる。そのため、虚血性心疾患の病態の解明には力学解析が重要である。

【0005】

このような冠循環特有の血行動態は基礎医学の面からも強い関心を集め、これまでに多彩な実験的研究がなされ多くの知見が集積されてきたが、心臓の構造や動きは複雑であり、不明な点も多く残されている。例えば、冠血管は主幹部から毛細血管にいたるまで、その径はμｍからｍｍへと約１０００倍ものスケール差を有する。また心臓は収縮、弛緩を繰り返し大きく拍動するため冠血管の観察には大きな制約がある。そのため、特に壁内を走行する血管内の血行動態については、未だ不明な点が多い。コンピュータ上に作成された冠循環モデルを利用したコンピュータシミュレーションは、冠血管観察における技術的な制約を越えた現象の把握を可能とし、心臓のメカニズム解明に迫る有効な手段である。

【0006】

現在では、冠循環シミュレーションに用いるため、冠動脈、冠静脈、毛細血管の血管径や長さ、分岐情報が詳細に計測され、データベース化されている。このデータベースを用いて、さまざまなシミュレーション手法が考えられている。例えば、心臓の冠循環シミュレーションにおいて生ずる連立一次方程式の並列反復解法として、冠循環系における血管径の大きな違いから生じる悪条件性を解消し、安定した収束性を維持する並列前処理が提案されている。

【先行技術文献】

【非特許文献】

【0007】

【非特許文献1】G.S.Kassab, C.A.Rider, N.J.Tang and Y.C.Fung, “Morphometry of pig coronary arterial trees”, American Physiological Society Heart and Circulatory Physiology, 1993, 265:H350-H365.

【非特許文献2】G.S.Kassab, D.H.Lin and Y.C.Fung, “Morphometry of pig coronary venous system”, American Physiological Society Heart and Circulatory Physiology, 1994, 267: H2100-H2113.

【非特許文献3】G.S.Kassab and Y.C.Fung, “Topology and dimensions of pig coronary capillary network”, American Physiological Society Heart and Circulatory Physiology, 1994, 267:H319-H325.

【非特許文献4】鷲尾，岡田，久田，“心臓冠循環シミュレーションにおける並列前処理法について”，計算工学講演会論文集，計算工学会，2009年5月，Vol.14

【非特許文献5】Cortassa S, Aon M A, O'Rourke B, Jacques R, Tseng H-J, Marban E, Winslow R L, “A computational model Integrating electrophysiology, contraction, and mitochondrial bioenergetics in the ventricular myocyte”, Biophys J 91, pp. 1564-1589, 2006

【発明の概要】

【発明が解決しようとする課題】

【0008】

しかし、従来技術では、冠動脈・冠静脈、さらには毛細血管にいたるまでの血管形態を用い、かつ心臓収縮と連成させた血行動態解析は困難である。その主な原因は、毛細血管までも忠実にデータ化することによりデータ量が膨大になり、現実的な時間でのシミュレーション実行が困難であることによる。このため、従来は毛細血管レベルの血管を多孔質体で代替する手法が用いられていた。しかしながら、従来の手法では、・毛細血管レベルの圧力下降の解析、・微小循環の構造が冠血流に与える影響、・代謝機能と微小循環との連成解析などの解析は困難である。なお、代謝機能とは、心筋細胞が毛細血管内の血流から供給される酸素・栄養を消費し、アデノシン三リン酸（ＡＴＰ：Adenosine TriPhosphate）と呼ばれる細胞の諸活動に必要な物質を産生する機能である。

【0009】

以上、心臓の冠循環シミュレーションについて説明したが、心臓以外の器官においても、血流は、器官の動きに連動して変化する。そのため、心臓以外の器官であっても、毛細血管にいたるまでの血管形態と、その器官の動きとの連成解析による血行動態シミュレーションの困難性が存在する。

【0010】

１つの側面では、本発明は、効率的な血行動態シミュレーションを可能とするシミュレーション方法、シミュレーション装置、およびシミュレーションプログラムを提供することを目的とする。

【課題を解決するための手段】

【0011】

１つの案では、コンピュータが、複数の節点間を接続して得られる複数の要素により器官の形状を定義した形状モデルと、前記器官内の閾値以上の径の血管による血管ネットワークを定義した第１の血管ネットワークモデルとに基づいて、前記器官の動きに応じた前記第１の血管ネットワークモデルの血管内の血液の状態を解析し、前記器官内の前記閾値未満の径の血管による血管ネットワークを定義した第２の血管ネットワークモデルが、前記複数の節点それぞれに、１つまたは複数接続されており、前記第１の血管ネットワークモデルの血管内の血液の状態の解析で得られた、前記複数の節点それぞれにおける血液の状態を示す情報を用いて、前記複数の節点それぞれに接続された第２の血管ネットワークモデルの血液の状態を、解析することを特徴とするシミュレーション方法が提供される。

【発明の効果】

【0012】

１態様によれば、毛細血管レベルの圧力下降の解析、微小循環の構造が心臓に与える影響、代謝と微小循環との連成解析などが実現でき、効率的な血行動態シミュレーションが可能となる。

【図面の簡単な説明】

【0013】

【図1】第１の実施の形態に係る並列計算システムの機能構成の一例を示す図である。

【図2】太い冠動脈・静脈の血管ネットワークモデルの一例を示す図である。

【図3】微小循環モデルの一例を示す図である。

【図4】対称性モデルの一例を示す図である。

【図5】ミクロモデルの一例を示す図である。

【図6】圧力節点におけるマクロモデルとミクロモデルとの接続法を示す図である。

【図7】２分岐血管モデルの一例を示す図である。

【図8】第２の実施の形態のシステム構成例を示す図である。

【図9】ＣＡＤシステムのハードウェア構造の一例を示す図である。

【図10】計算ノードのハードウェア構成の一例を示す図である。

【図11】第２の実施の形態の各装置の機能を示すブロック図である。

【図12】有限要素モデルのデータ構造の一例を示す図である。

【図13】マクロモデルのデータ構造の一例を示す図である。

【図14】ミクロモデルのデータ構造の一例を示す図である。

【図15】冠循環シミュレーションの手順を示すフローチャートの一例である。

【図16】冠循環へのレベル付けの一例を示す図である。

【図17】マクロモデル作成処理の手順の一例を示すフローチャートである。

【図18】心表面冠血管モデルの表示例を示す図である。

【図19】中隔枝を付与した心表面冠循環モデルの一例を示す図である。

【図20】ミクロモデル作成処理の手順の一例を示すフローチャートである。

【図21】微小循環モデルの作成過程の一例を示す図である。

【図22】対称性モデルの末端に微小循環モデルを１２８個並列接続したミクロモデルの一例を示す図である。

【図23】血管密度比の設定結果を表した図である。

【図24】血管径の違いを示す図である。

【図25】心臓の冠循環シミュレーションの手順の一例を示すフローチャートである。

【図26】第４の実施の形態におけるシミュレーション手順の一例を示すフローチャートである。

【図27】第５の実施の形態におけるシステム構成の一例を示す図である。

【発明を実施するための形態】

【0014】

以下、本実施の形態について図面を参照して説明する。なお各実施の形態は、矛盾のない範囲で複数の実施の形態を組み合わせて実施することができる。
〔第１の実施の形態〕
まず、第１の実施の形態について説明する。第１の実施の形態は、心臓などの器官の解剖学的特徴を活かすことにより、マルチスケール解析の分散処理手法を適用し、従来の手法よりも高速に計算することを可能としたものである。例えば心臓には、微小循環系どうしは直接接続しておらず、その上流側或いは下流側の分岐を通じてのみ互いに繋がっているという解剖学的特徴がある。この解剖学的特徴を活かすことで、シミュレーションの高速化が可能となる。さらに第１の実施の形態では、微小循環と太い血管の中間スケールの血管に対しては対称性を仮定する。これにより、計算量を大幅に削減することができる。

【0015】

図１は、第１の実施の形態に係る並列計算システムの機能構成の一例を示す図である。シミュレーション装置Ｘは、記憶手段１と複数の演算装置２，３，４，５，・・・とを有している。図１に示したシミュレーション装置Ｘは、複数の演算装置２，３，４，５，・・・による並列計算を行えば、並列計算システムとして機能する。

【0016】

記憶手段１は、形状モデル１ａ、第１の血管ネットワークモデル１ｂ、および第２の血管ネットワークモデル１ｃを記憶する。形状モデル１ａは、複数の節点間を接続して得られる複数の要素により器官の形状を定義した情報である。第１の血管ネットワークモデル１ｂは、器官内の閾値以上の径を有する血管による血管ネットワークを定義した情報である。第２の血管ネットワークモデル１ｃは、器官内の閾値未満の径を有する血管による血管ネットワークを定義した情報である。なおシミュレーション対象の器官が心臓であれば、第１の血管ネットワークモデル１ｂと第２の血管ネットワークモデル１ｃとで示される血管ネットワークは、心臓の冠循環ネットワークを示していることとなる。

【0017】

複数の演算装置２，３，４，５，・・・のうちの１つの演算装置２は、第１の解析部２ａを有する。第１の解析部２ａは、形状モデル１ａと第１の血管ネットワークモデル１ｂとに基づいて、器官の動きに応じた第１の血管ネットワークモデル１ｂの血管内の血液の状態を解析する。これにより、形状モデル１ａの複数の節点における、圧力などの血液の状態を示す情報が得られる。なお器官が心臓であれば、器官の動きは、心筋の収縮・弛緩による動きである。なお、器官の動き（例えば心臓の収縮・弛緩）の影響を直接的に受けるのは、第２の血管ネットワークモデル１ｃである。第１の解析部２ａは、例えば、第１の血管ネットワークモデル１ｂの血流の解析と、器官の動きの解析を合わせて、複数の節点における圧力などの物理量の状態を解析する。

【0018】

第１の解析部２ａを有する演算装置２以外の演算装置３，４，５，・・・は、それぞれ第２の解析部３ａ，４ａ，５ａ，・・・を有する。第２の解析部３ａ，４ａ，５ａ，・・・は、節点における血液の状態を示す情報を用いて、その節点に接続された第２の血管ネットワークモデル１ｃの血液の状態を解析する。例えば、形状モデル１ａの複数の節点のそれぞれに、第２の血管ネットワークモデル１ｃが１つまたは複数接続されている。節点に接続された第２の血管ネットワークモデル１ｃは、第２の解析部３ａ，４ａ，５ａ，・・・それぞれに、解析対象として割り当てられる。そして第２の解析部３ａ，４ａ，５ａ，・・・は、割り当てられた第２の血管ネットワークモデル１ｃが接続された節点の圧力の情報を用いて、その第２の血管ネットワークモデル１ｃの血液の状態を解析する。すなわち複数の演算装置３，４，５，・・・が、複数のミクロモデルの血液の状態を、並列処理で解析する。

【0019】

なお、演算装置２，３，４，５，・・・は、シミュレーション装置Ｘが有する演算処理装置としてのＣＰＵ（Central Processing Unit）、またはマルチコアＣＰＵ内の個々のコアである。また、記憶手段１は、シミュレーション装置Ｘが有するＲＡＭ（Random Access Memory）またはハードディスクドライブ（ＨＤＤ:Hard Disk Drive）などの記録媒体である。

【0020】

また、図１に示した各要素間を接続する線は通信経路の一部を示すものであり、図示した通信経路以外の通信経路も設定可能である。
以下、第１の実施の形態に係るシミュレーション装置Ｘで心臓の冠循環シミュレーションを行う場合を例にとり、血行動態シミュレーションについて詳細に説明する。なお、以下の説明では、第１の血管ネットワークモデル１ｂをマクロモデル、第１の血管ネットワークモデル１ｃをミクロモデルと呼ぶこととする。また第２の解析部３ａ，４ａ，５ａ，・・・は、節点における圧力に関する情報を用いて、ミクロモデルの血液の状態を解析するものとする。そこで以下の説明では、ミクロモデルが接続される節点を、圧力節点と呼ぶ。さらに、説明で使用する用語を、以下のように定義する。

【0021】

＜用語の定義＞
−太い冠動脈・静脈：直径約１００μｍ以上の血管（図２参照）
−中間スケールの血管：直径約２０μｍ−１００μｍの血管（図４参照）
−細動脈・細静脈：直径約１０μｍ−２０μｍの血管（図３参照）
−毛細血管：直径約５−７μｍの血管で、実際に心筋細胞との酸素・栄養素の交換が行われる場（図３参照）
−微小循環系：細動脈、毛細血管、細静脈より構成される心筋細胞との酸素・栄養素のやりとりを行う循環系（図３参照）
＜マルチスケールでの血管ネットワークモデル作成＞
第１の実施の形態では、太い冠動脈・静脈から毛細血管までを含むすべてのスケールの血管ネットワークを考える。すべてのスケールの血管を含む血管ネットワークは、大動脈から発する太い冠動脈から分岐を繰り返して小動脈、細静脈、最終的に毛細血管にいたる。さらに血管ネットワークは、毛細血管は細静脈で集合し、集合を繰り返しながら最終的に太い冠静脈となり右心房に達する。なお、血管が分岐しているとき、その分岐点から太い方の血管を、上位の血管、細い方の血管を下位の血管と呼ぶこととする。

【0022】

このような血管ネットワークをモデリングするとき、まず、太い冠動脈の血管ネットワークと太い冠静脈の血管ネットワークモデルとを別々に作成する。
図２は、太い冠動脈・静脈の血管ネットワークモデルの一例を示す図である。図２中、左側が太い冠動脈のネットワークモデル１１であり、右側が太い冠静脈のネットワークモデル１２である。太い冠動脈のネットワークモデル１１と太い冠静脈のネットワークモデル１２とを合わせた血管ネットワークモデルが、マクロモデル１０である。

【0023】

次に、細動脈・細静脈と毛細血管とで構成される微小循環系を、解剖学的データに基づいた微小循環モデルで表現する。
図３は、微小循環モデルの一例を示す図である。図３の例では、微小循環モデル２０の左側の血管ネットワークが細動脈２１であり、右側の血管ネットワークが細静脈２２である。そして、細動脈２１と細静脈２２とを接続する血管ネットワークが、毛細血管２３である。

【0024】

さらに微小循環モデル２０とマクロモデル１０との間のスケールに属する中間の血管ネットワークモデルを作成する。第１の実施の形態では、中間の血管ネットワークは、動脈系と静脈系との間で対称性を有するものと仮定する。以下、動脈系と静脈系との間で対称性を有するものと仮定することで得られる、マクロモデルと微小循環モデルとの間の血管ネットワークモデルを、対称性モデルと呼ぶ。

【0025】

図４は、対称性モデルの一例を示す図である。対称性モデル３０には、動脈系中間モデル３１と静脈系中間モデル３２とが含まれる。動脈系中間モデル３１と静脈系中間モデル３２とは、毛細血管を挟んで対称となるように作られている。そのため動脈系中間モデル３１と静脈系中間モデル３２とでは、上位の血管から分岐する２つの下位の血管の血管径、長さ、その後の分岐形態が同様である。このように対称性を仮定することにより、血管ネットワークの本質的な特性を失うことなく、計算量を大幅に削減することが可能となる。

【0026】

作成された対称性モデル３０の上端はマクロモデル１０に接続される。また対称性モデル３０の末端は、微小循環モデル２０に接続される。微小循環モデルと対称モデルを包含する血管ネットワークモデルが、ミクロモデルである。

【0027】

図５は、ミクロモデルの一例を示す図である。ミクロモデル４０は、動脈系中間モデル３１、静脈系中間モデル３２、および複数の微小循環モデル２０ａ，・・・，２０ｎを含んでいる。

【0028】

次に、マクロモデル１０とミクロモデル４０とを接続する。マクロモデル１０とミクロモデル４０との接続点は、仮想的に心筋の圧力節点ｋと同じ位置にあると考える。圧力節点ｋは、心臓の有限要素モデルを構成する各多面体要素の頂点である。そして、ミクロモデル４０が圧力節点ｋにおいて心筋から受ける外圧を、心筋節点圧Ｐ_kとする。

【0029】

圧力節点ｋへのミクロモデル４０の接続数は、圧力節点ごとの心筋節点集中体積Ｖ_k^muscleより決定する。各圧力節点の心筋節点集中体積Ｖ_k^muscleは、要素の体積を示す値をその要素を構成する圧力節点に分配し、圧力節点ごとに、その圧力節点に分配された値を合計した値である。例えば４面体の要素の場合、その要素の体積の１／４の値を、その要素を構成する４つの圧力節点それぞれの心筋節点集中体積Ｖ_k^muscleに分配する。圧力節点ｋが、ｍ個（ｍは１以上の整数）の要素で共有されている場合、その圧力節点ｋの心筋節点集中体積Ｖ_k^muscleは、ｍ個の要素から分配された値の合計となる。

【0030】

なお心臓の解剖学的特徴として、外側より内側の方が、血管密度が高いことが分かっている。そこで、圧力節点ｋへのミクロモデル４０の接続数の決定に際し、貫壁性（心臓の外壁からの深さ）に応じた、心筋容積に対する細血管の密度の違いを反映させる。例えば圧力節点の心筋節点集中体積Ｖ_k^muscleに、貫壁性に応じた血管密度比を乗算する。血管密度比は、心臓の内部に近い圧力節点ほどを高い値とする。

【0031】

次に、心臓容積に対するミクロモデル血管の総容積の割合より、圧力節点に単位体積当たりに含まれるミクロモデル血管容積ρ_kを得る。ミクロモデル１個の血管総容積をＶ^microとすると、マクロモデルの末端の圧力節点ｋに接続されるミクロモデルの数ｎ_kは次式で決定できる。

【0032】

【数1】

【0033】

圧力節点に接続されるミクロモデルの数ｎ_kが確定すると、マクロモデルとミクロモデルの動脈側、静脈側接続部において次式の流量保存則が成立する。Ｑγ_macroは、マクロモデルから、ミクロモデルに送出される血液の流量を示す。Ｑγ_microは、ミクロモデルから、マクロモデルへ送出される血液の流量を示す。γは動脈と静脈とを区別する符号である。

【0034】

【数2】

【0035】

図６は、圧力節点におけるマクロモデルとミクロモデルとの接続法を示す図である。有限要素モデルの要素５１は、複数の圧力節点間を線分で接続して得られる多面体である。要素５１の圧力節点ｋ付近を通過するマクロモデル１０の血管に対して、要素５１の圧力節点ｋを介して、ｎ_k個のミクロモデル４０が接続されている。なお、図６に示されているマクロモデル１０の血管ネットワークは、動脈系のネットワークであるものとする。

【0036】

すると、マクロモデルから、圧力節点を介してミクロモデルに送出される血液の流量Ｑγ_macroは、ｎ_k個のミクロモデル４０に分配される。また、ミクロモデル４０から、圧力節点を介してマクロモデルへ送出される血液の流量Ｑγ_microは、ｎ_k個分が合成され、マクロモデル１０に流入する。なお、例えばミクロモデル４０において、マクロモデル１０から血流が流入している場合、流量Ｑγ_microは負の値となる。そのため、流量Ｑγ_macroと流量Ｑγ_microのｎ_k倍との加算値は、式（２）に示したように０になる。

【0037】

＜計算手法＞
次に、心臓の冠循環シミュレーションの計算手法について説明する。なお、以下の説明は、Newton-Raphson法による非線形解析を行う場合の例である。

【0038】

第１の実施の形態では、マルチスケール解析の考え方を用いて並列計算により解く。全体のマトリクス（行列）は式（３）で与えられる。

【0039】

【数3】

【0040】

Ａ_Micは血管ミクロモデルより作成される係数マトリクスである。Ｂ_Macは血管マクロモデルにより作成される係数マトリクスである。Ｃはミクロ部とマクロ部の接続を表すマトリクスである。Ｃ^Tは、Ｃの転置マトリクスである。Δμ_Micは、非線形解析におけるミクロモデルの圧力の増分である。Δμ_Macは、非線形解析におけるマクロモデルの圧力の増分である。ｒ_Micは、ミクロモデルの非線形解析で生じる残差である。ｒ_Macは、マクロモデルの非線形解析で生じる残差である。

【0041】

式（３）の係数マトリクスは、以下のように分解できる。

【0042】

【数4】

【0043】

式（４）上式に基づき以下の３ステップでミクロ、マクロの圧力の増分を計算する。

【0044】

【数5】

【0045】

【数6】

【0046】

【数7】

【0047】

式（５）に示す第１のステップでは、ミクロモデルの計算式を、マクロの圧力の増分を除いた形式に置き換えている。これにより、ミクロモデルにおける圧力の仮の増分Δμ_Mic（μはチルダ付き）が得られる。式（６）に示す第２のステップでは、ミクロモデルで得られた圧力の仮の増分を用いて、マクロモデルの圧力の増分を計算している。式（７）に示す第３のステップでは、マクロモデルの圧力の増分を用いて、ミクロモデルの本来の圧力の増分値を計算している。

【0048】

以上のように、マクロモデルのマトリクスに入り込むミクロモデルの自由度を消去することで、マクロモデルとミクロモデルを分けて計算することが可能となる。さらにミクロモデル同士は完全に独立しているため、並列計算においてはミクロモデルを各演算装置が独立して求解することが可能である。その結果、効率的な計算が実現できる。

【0049】

＜対称性を考慮することによる計算量の削減＞
既に述べたように、第１の実施の形態では、毛細血管とマクロモデルとの中間層にあたる中間スケールの血管は、対称性を仮定した対称性モデル３０を用いて表現している。このような対称性モデル３０を用いることで、計算量の削減が可能である。以下、最も簡単な２分岐の血管（３ 要素、４節点）を例に、対称性の考え方を説明する。

【0050】

図７は、２分岐血管モデルの一例を示す図である。図７に示す２分岐血管モデル４１は、動脈系の血管に関するモデルである。２分岐血管モデル４１には、上位の血管が２つの下位の血管に分岐している。各血管の端部の数字は、節点番号である。ここで、上位の血管には、流量Ｑ₁の血流が流れ込むものとする。また２つの下位の血管からは、それぞれ流量Ｑ₂、Ｑ₃の血流が流れ出すものとする。

【0051】

ここで、各血管のコンダクタンスを定義する。コンダクタンスは、血液の流れやすさを示す値であり、抵抗の反対の概念である。コンダクタンスが大きい血管ほど、血液が流れやすい。以下、ｉ番目（１は１以上の整数）の節点と、ｊ番目（ｊは１以上の整数）との間の血管のコンダクタンスを、「ｋ_ij」と表記する。例えば、節点番号「１」の節点と節点番号「２」の節点との間の血管のコンダクタンスは、ｋ₁₂と表記される。

【0052】

簡単のために血管のコンダクタンスのみを考慮して、２分岐の血管の血流は、以下の式（８）で表される。

【0053】

【数8】

【0054】

この式（８）に基づき剛性マトリクスを作成すると式（９）のようになる。

【0055】

【数9】

【0056】

Ｑ_iは節点ｉに流入する血液の流量である。Ｆ_iは節点ｉに接続する節点の集合である。μ_iは節点ｉの血圧である。Ｖは血管の容積である。
今、上位の血管から分岐している２本の下位の血管に対称性を仮定する。つまり下位の血管の血管径、長さ、その後の分岐形態がまったく同じであると考える。すると、
ｋ₂₃＝ｋ₂₄
μ₃＝μ₄
Ｑ₃＝Ｑ₄
が成立する。そのため、 以下のようにマトリクスを縮退させることができる。

【0057】

【数10】

【0058】

これにより、計算量の削減が可能となる。
以上説明したように、第１の実施の形態では、血管ネットワークの解剖学的特徴を利用してマルチスケール解析を導入している。すなわち、微小循環系同士は、より太い血管を介してのみ接続されているという特徴に基づいて、微小循環系を含むミクロモデル４０を定義している。そして形状モデルの節点それぞれに、マクロモデル１０にミクロモデル４０を接続している。これにより、毛細血管を含む微小循環モデル２０が、統一したミクロモデル４０で解析できるため、データ量が少なくなる。また、マクロモデル１０に対するミクロモデル４０の接続位置が、圧力節点に制限されていることで、第１の解析部２ａが、複数の第２の解析部３ａ，４ａ，５ａ，・・・から渡される情報の位置に関する自由度が少なくなり、処理量も少なくなる。しかも多数の節点に接続されたミクロモデルは互いに独立して解析できるため、並列処理が可能である。

【0059】

このようにして、第１の実施の形態では、並列処理に適した計算アルゴリズムを導出し、シミュレーション計算時間の短縮が可能となっている。その結果、これまで実現されていなかった、血管の形態を毛細管レベルまで再現し、かつ心筋の拍動と連成したシミュレーションを可能となる。

【0060】

さらに、第１の実施の形態では、微小循環と太い血管の中間層の血管に対して対称性があるものと仮定することで計算量を削減している。その結果、シミュレーションのさらなる高速化が図られている。

【0061】

〔第２の実施の形態〕
次に第２の実施の形態について説明する。第２の実施の形態は、マルチコアＣＰＵを有する計算ノードを複数使用した並列計算システムで、心臓の冠循環シミュレーションを行うものである。

【0062】

図８は、第２の実施の形態のシステム構成例を示す図である。第２の実施の形態では、ＣＡＤ（Computer Aided Design）システム１００が設けられている。ＣＡＤシステム１００は、ユーザからの操作に基づき、心臓の立体構造や冠循環血管ネットワークをモデリングする。またＣＡＤシステム１００は、ネットワーク６１を介して管理ノード２００に接続されている。ＣＡＤシステム１００は、心臓の立体構造や冠循環血管ネットワークのモデルを、管理ノード２００に送信し、そのモデルに基づく冠循環シミュレーションの実行を指示する。

【0063】

管理ノード２００はネットワーク６２を介して、ストレージ装置３００や、複数の計算ノード４００，５００，６００，・・・に接続されている。管理ノード２００は、ＣＡＤシステム１００から受信した心臓の立体構造や冠循環血管ネットワークのモデルを、ストレージ装置３００に格納する。また管理ノード２００は、ＣＡＤシステム１００から冠循環シミュレーションの実行指示を受信すると、冠循環シミュレーションを実行するためのジョブを、計算ノード４００，５００，６００，・・・に投入する。ここでジョブの投入処理では、例えば管理ノード２００から計算ノードへ、モデルデータの位置情報、ジョブを実行するためのプログラム、ジョブの実行に使用する各種パラメータなどが送信される。

【0064】

ストレージ装置３００は、ＣＡＤシステム１００で作成されたモデルを記憶する。ストレージ装置３００には、複数のハードディスク装置が実装されている。またストレージ装置３００は、ハードディスク装置に代えて、ＳＳＤ（Solid State Drive）を有していてもよい。さらにストレージ装置３００には、ＲＡＩＤ（Redundant Array of Independent Disks）技術を適用することもできる。

【0065】

計算ノード４００，５００，６００，・・・は、管理ノード２００からの指示に従って、心臓の冠循環シミュレーションを並列処理で実行する。計算ノード４００，５００，６００，・・・は、マルチコアＣＰＵを有しており、個々の計算ノード内部でも、処理を並列に実行することができる。

【0066】

図９は、ＣＡＤシステムのハードウェア構造の一例を示す図である。ＣＡＤシステム１００は、ＣＰＵ１０１によって装置全体が制御されている。ＣＰＵ１０１には、バス１０９を介してＲＡＭ１０２と複数の周辺機器が接続されている。なおＣＡＤシステム１００が有するＣＰＵ数は１つに限定されず、複数であってもよい。ＣＡＤシステム１００が複数のＣＰＵを有する場合、複数のＣＰＵが連係して動作し、装置全体を制御する。

【0067】

ＲＡＭ１０２は、ＣＡＤシステム１００の主記憶装置として使用される。ＲＡＭ１０２には、ＣＰＵ１０１に実行させるＯＳ（Operating System）のプログラムやアプリケーションプログラムの少なくとも一部が格納される。また、ＲＡＭ１０２には、ＣＰＵ１０１による処理に必要な各種データが格納される。

【0068】

バス１０９に接続されている周辺機器としては、ＨＤＤ１０３、グラフィック処理装置１０４、入力インタフェース１０５、光学ドライブ装置１０６、機器接続インタフェース１０７およびネットワークインタフェース１０８がある。

【0069】

ＨＤＤ１０３は、内蔵したディスクに対して、磁気的にデータの書き込みおよび読み出しを行う。ＨＤＤ１０３は、ＣＡＤシステム１００の補助記憶装置として使用される。ＨＤＤ１０３には、ＯＳのプログラム、アプリケーションプログラム、および各種データが格納される。なお、補助記憶装置としては、フラッシュメモリなどの半導体記憶装置を使用することもできる。

【0070】

グラフィック処理装置１０４には、表示装置としてモニタ６３が接続されている。グラフィック処理装置１０４は、ＣＰＵ１０１からの命令に従って、画像をモニタ６３の画面に表示させる。モニタ６３としては、ＣＲＴ（Cathode Ray Tube）を用いた表示装置や液晶表示装置などがある。

【0071】

入力インタフェース１０５には、入力装置としてキーボード６４とマウス６５とが接続されている。入力インタフェース１０５は、キーボード６４やマウス６５から送られてくる信号をＣＰＵ１０１に送信する。なお、マウス６５は、ポインティングデバイスの一例であり、他のポインティングデバイスを使用することもできる。他のポインティングデバイスとしては、タッチパネル、タブレット、タッチパッド、トラックボールなどがある。

【0072】

光学ドライブ装置１０６は、レーザ光などを利用して、光ディスク６６に記録されたデータの読み取りを行う。光ディスク６６は、光の反射によって読み取り可能なようにデータが記録された可搬型の記録媒体である。光ディスク６６には、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）などがある。

【0073】

機器接続インタフェース１０７は、ＣＡＤシステム１００に周辺機器を接続するための通信インタフェースである。例えば機器接続インタフェース１０７には、メモリ装置６７やメモリリーダライタ６８を接続することができる。メモリ装置６７は、機器接続インタフェース１０７との通信機能を搭載した記録媒体である。メモリリーダライタ６８は、メモリカード６９へのデータの書き込み、またはメモリカード６９からのデータの読み出しを行う装置である。メモリカード６９は、カード型の記録媒体である。

【0074】

ネットワークインタフェース１０８は、ネットワーク６１に接続されている。ネットワークインタフェース１０８は、ネットワーク６１を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。

【0075】

図１０は、計算ノードのハードウェア構成の一例を示す図である。計算ノード４００は、ＣＰＵ４１０、メモリ４２０、およびインターコネクト回路４３０を有している。ＣＰＵ４１０は、複数のコア４１１〜４１８を有する演算処理装置としてのマルチコアＣＰＵである。複数のコア４１１〜４１８により、処理を並列に実行することができる。メモリ４２０は、ＣＰＵ４１０に実行させるプログラムや、処理の実行に用いられるデータを記憶する。例えば、主記憶装置としてのメモリ４２０は、心臓の冠循環シミュレーションを実行するためのプログラムを記憶する。インターコネクト回路４３０は、光による高速通信を行う通信回路である。インターコネクト回路４３０は、ネットワーク６２を介して、他の計算ノード４００，５００，６００，・・・とデータの受け渡しを行う。またインターコネクト回路４３０は、ネットワーク６２を介してストレージ装置３００から心臓のモデルデータを読み出し、ＣＰＵ４１０に転送する。

【0076】

以上のようなハードウェア構成によって、第２の実施の形態の処理機能を実現することができる。
なおＣＡＤシステム１００、管理ノード２００、および計算ノード４００，５００，６００，・・・は、コンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態の処理機能を実現する。ＣＡＤシステム１００、管理ノード２００、および計算ノード４００，５００，６００，・・・それぞれに実行させる処理内容を記述したプログラムは、さまざまな記録媒体に記録しておくことができる。例えば、ＣＡＤシステム１００、管理ノード２００、および計算ノード４００，５００，６００，・・・それぞれに実行させるプログラムをＨＤＤなどに格納しておくことができる。各装置のＣＰＵは、ＨＤＤ内のプログラムの少なくとも一部をＲＡＭにロードし、プログラムを実行する。またＣＡＤシステム１００、管理ノード２００、および計算ノード４００，５００，６００，・・・それぞれに実行させるプログラムを、光ディスク６６、メモリ装置６７、メモリカード６９などの可搬型記録媒体に記録しておくこともできる。なおプログラムを記録する記録媒体には、一時的な伝搬信号自体は含まれない。

【0077】

プログラムを流通させる場合には、例えば、そのプログラムが記録された光ディスク６６、メモリ装置６７、メモリカード６９などの可搬型記録媒体が販売される。また、プログラムをコンピュータの記憶装置に格納しておき、ネットワーク６１を介して、他のサーバコンピュータからＣＡＤシステム１００、管理ノード２００、または計算ノード４００，５００，６００，・・・にそのプログラムを転送することもできる。ＣＡＤシステム１００、管理ノード２００、または計算ノード４００，５００，６００，・・・は、ネットワークを介してプログラムを取得する場合、例えば取得したプログラムをＲＡＭまたはＨＤＤに格納する。

【0078】

なお、第１の実施の形態に示したシミュレーション装置Ｘも、図８〜図１０に示したシステムと同様のハードウェア構成により実現することができる。
次に、第２の実施の形態のシステムにおける各装置の機能について説明する。

【0079】

図１１は、第２の実施の形態の各装置の機能を示すブロック図である。ＣＡＤシステム１００は、有限要素モデル作成部１１０、冠循環モデル作成部１２０、およびシミュレーション条件決定部１３０を有する。有限要素モデル作成部１１０は、ユーザからの指示に従って、心臓の立体形状を示す有限要素モデルを作成する。有限要素モデルは、例えば複数の４面体のメッシュで構成される。４面体のメッシュは、第１の実施の形態における多面体の要素の一例である。冠循環モデル作成部１２０は、ユーザからの指示に従って、冠循環モデルを作成する。冠循環モデルは、心臓の血管ネットワークを表すモデルである。シミュレーション条件決定部１３０は、ユーザからの指示に従って、心臓の冠循環シミュレーションの実行条件を決定する。例えば、シミュレーション条件決定部１３０は、心筋の各圧力節点に、どれくらいの数のミクロモデルを接続させるかを決定する。

【0080】

管理ノード２００は、実行指示部２１０を有する。実行指示部２１０は、ＣＡＤシステム１００で作成された有限要素モデルと冠循環モデルとを用いた心臓の冠循環シミュレーションの実行を、計算ノード４００，５００，６００，・・・に指示する。例えば実行指示部２１０は、ＣＡＤシステム１００から、有限要素モデルと冠循環モデルとを取得し、ストレージ装置３００に格納する。次に実行指示部２１０は、冠循環シミュレーションのジョブを投入する計算ノードを選択する。例えば実行指示部２１０は、ジョブを実行していないすべての計算ノードを、ジョブの投入先として選択することができる。また、予め冠循環シミュレーションの並列度が指定されている場合、実行指示部２１０は、指定された並列度応じた数の計算ノードを、ジョブの投入先として選択する。さらに実行指示部２１０は、選択した計算ノードの１つをマクロ解析用とし、他の計算ノードをミクロ解析用とする。そして実行指示部２１０は、マクロ解析用の計算ノード（図１１の例では計算ノード４００）に、シミュレーション全体の制御とマクロ解析を行うジョブを投入する。また実行指示部２１０は、ミクロ解析用の計算ノード（図１１の例では計算ノード５００，６００，・・・）に、ミクロ解析を行うジョブを投入する。

【0081】

ストレージ装置３００には、有限要素モデル３１０、マクロモデル３２０、およびミクロモデル３３０が格納される。有限要素モデル３１０は、図１に示した第１の実施の形態の形状モデル１ａの一例である。マクロモデル３２０は、図１に示した第１の実施の形態の第１の血管ネットワークモデル１ｂの一例である。ミクロモデル３３０は、図１に示した第１の実施の形態の第２の血管ネットワークモデル１ｃの一例である。なおマクロモデル３２０とミクロモデル３３０とを合わせたものが、冠循環モデルである。

【0082】

計算ノード４００は、シミュレーション制御部４４０とマクロ解析部４５０とを有する。シミュレーション制御部４４０は、心臓の冠循環シミュレーション全体の進行を制御する。マクロ解析部４５０は、冠循環シミュレーションのマクロ部分の解析を行う。例えばマクロ解析部４５０は、心筋の収縮・弛緩による心臓の拍動の解析と、太い冠動脈・静脈の血流の解析とを行う。シミュレーション制御部４４０とマクロ解析部４５０とは、計算ノード４００にマクロ解析用のジョブが投入されることにより実現される機能である。また、マクロ解析部４５０は、図１に示した第１の実施の形態の第１の解析部２ａの一例である。

【0083】

計算ノード５００は、複数のミクロ解析部５１０，５２０，５３０，・・・を有する。ミクロ解析部５１０，５２０，５３０，・・・は、中間スケールの血管、細動脈・細静脈、および毛細血管の血流の解析を行う。複数のミクロ解析部５１０，５２０，５３０，・・・は、ＣＰＵ内の個別のコアによって実行される。ミクロ解析部５１０，５２０，５３０，・・・は、計算ノード５００にミクロ解析用のジョブが投入されることで実現される機能である。

【0084】

計算ノード６００は、複数のミクロ解析部６１０，６２０，６３０，・・・を有する。ミクロ解析部６１０，６２０，６３０，・・・は、中間スケールの血管、細動脈・細静脈、および毛細血管の血流の解析を行う。複数のミクロ解析部６１０，６２０，６３０，・・・は、ＣＰＵ内の個別のコアによって実行される。ミクロ解析部６１０，６２０，６３０，・・・は、計算ノード６００にミクロ解析用のジョブが投入されることで実現される機能である。

【0085】

なお計算ノード５００，６００，・・・のミクロ解析部は、図１に示した第１の実施の形態の第２の解析部３ａ，４ａ，５ａ，・・・の一例である。また図１１に示した各要素間を接続する線は通信経路の一部を示すものであり、図示した通信経路以外の通信経路も設定可能である。

【0086】

次に、ストレージ装置３００が記憶するモデルについて詳細に説明する。
図１２は、有限要素モデルのデータ構造の一例を示す図である。有限要素モデル３１０には、節点情報テーブル３１１とメッシュ情報テーブル３１２とが含まれる。なおＮ_nodesは節点数であり、Ｎ_meshはメッシュ数である。

【0087】

節点情報テーブル３１１は、４面体のメッシュの頂点となる節点の位置が登録されたデータテーブルである節点情報テーブル３１１には、節点番号に対応付けて、節点の座標が登録されている。節点番号は、節点を一意に識別するための識別番号である。座標は、節点の仮想３次元空間上でのｘ軸、ｙ軸、ｚ軸の座標である。

【0088】

メッシュ情報テーブル３１２は、メッシュの形状を定義する情報が登録されたデータテーブルである。メッシュ情報テーブル３１２には、メッシュ番号に対応付けて、４面体のメッシュを構成する節点の節点番号が設定されている。メッシュ番号は、メッシュを一意に識別するための識別番号である。

【0089】

図１３は、マクロモデルのデータ構造の一例を示す図である。マクロモデル３２０は、心臓の血管ネットワークのうち、太い冠動脈・静脈の構造を示す情報である。マクロモデル３２０には、節点情報テーブル３２１と血管要素情報テーブル３２２とが含まれる。節点情報テーブル３２１の構造は、図１２に示す有限要素モデル３１０の節点情報テーブル３１１と同様である。なおＮ_elemは血管の要素数である。

【0090】

血管要素情報テーブル３２２には、血管の要素の要素番号に対応付けて、節点、レベル、半径、およびコンダクタンスが設定されている。要素番号は、血管の要素を一意に識別するための識別番号である。節点は、血管の要素の両端の位置の節点の識別番号である。レベルは、血管のレベルである。レベルは、毛細血管を最も小さい値とし、太くなるほどに値が大きくなる。半径は、血管の半径である。コンダクタンスは、血管内の血液の流れやすさを示す数値である。コンダクタンスが大きいほど、血管内の抵抗か少なく、血液が流れやすいことを示す。

【0091】

図１４は、ミクロモデルのデータ構造の一例を示す図である。ミクロモデル３３０は、心臓の血管ネットワークのうち、中間ネットワーク、細動脈・細静脈、および毛細血管の構造を示す情報である。ミクロモデル３３０には、節点情報テーブル３３１と血管要素情報テーブル３３２とが含まれる。節点情報テーブル３３１の構造は、図１２に示す有限要素モデル３１０の節点情報テーブル３１１と同様である。なおＮ_elemは血管の要素数である。

【0092】

血管要素情報テーブル３３２には、血管の要素の要素番号に対応付けて、節点、レベル、長さ、直径、血管壁弾性率、基準圧力差、本数が設定されている。節点は、血管の要素の両端の位置の節点の識別番号である。レベルは、血管のレベルである。長さは、血管の要素の長さである。直径は、血管の要素の直径である。血管壁弾性率は、血管の要素の壁の弾性率である。基準圧力差は、血管内の血圧と心筋内圧の圧力差の基準値を表す。この基準値とシミュレーションにおける各時刻での圧力差を基に血管径を計算する。本数は、対称性があると仮定したときの同じレベルの要素の本数である。

【0093】

図１１に示した機能が連携して動作し、図１２〜図１４に示した情報を利用して、冠循環シミュレーションを実行する。
図１５は、冠循環シミュレーションの手順を示すフローチャートの一例である。以下、図１５に示す処理をステップ番号に沿って説明する。

【0094】

［ステップＳ１０１］ＣＡＤシステム１００の有限要素モデル作成部１１０は、ユーザからのＣＡＤに対する入力に従って、有限要素モデル３１０を作成する。そしてＣＡＤシステム１００は、作成した有限要素モデル３１０を管理ノード２００に送信する。管理ノード２００は、受信した有限要素モデル３１０をストレージ装置３００に格納する。

【0095】

［ステップＳ１０２］ＣＡＤシステム１００の冠循環モデル作成部１２０は、ユーザからのＣＡＤに対する入力に従って、マクロモデル３２０を作成する。そしてＣＡＤシステム１００は、作成したマクロモデル３２０を管理ノード２００に送信する。管理ノード２００は、受信したマクロモデル３２０をストレージ装置３００に格納する。

【0096】

［ステップＳ１０３］ＣＡＤシステム１００の冠循環モデル作成部１２０は、ユーザからのＣＡＤに対する入力に従って、ミクロモデル３３０を作成する。そしてＣＡＤシステム１００は、作成したミクロモデル３３０を管理ノード２００に送信する。管理ノード２００は、受信したミクロモデル３３０をストレージ装置３００に格納する。

【0097】

［ステップＳ１０４］管理ノード２００の実行指示部２１０は、複数の計算ノードに対して、冠循環シミュレーションを実行するジョブを投入し、ジョブを投入した計算ノードに対して冠循環シミュレーションを指示する。すると複数の計算ノードの複数のコアが、並列処理により冠循環シミュレーションを実行する。

【0098】

以下、図１５に示した手順に沿って、各ステップの処理を詳細に説明する。
＜有限要素モデル作成＞
例えばＣＡＤシステム１００の有限要素モデル作成部１１０は、医療画像（Computed Tomography, Magnetic Resonance Image, Echocardiogramなど）から、心臓の形状を抽出する。そして有限要素モデル作成部１１０は、心臓の形状を示すデータを、有限要素法の実行に必要なメッシュ（例えば４面体または６面体）に分割する。

【0099】

＜冠循環ネットワークモデル作成＞
冠循環ネットワークモデルの作成は、マクロモデルの作成と、ミクロモデルの作成とに分かれる。マクロモデルとミクロモデルの分類には、例えば、Kassabが計測し、血管径に基づいてレベル付けしたデータを利用することができる（非特許文献１−３参照）。
・マクロモデルに含まれる動脈のレベルを１１〜６、静脈のレベルを−６〜−１２とする。
・ミクロモデルのうち、対称性を仮定する中間スケールの血管の動脈側のレベルを５〜３、静脈側のレベルを−３〜−５とし、微小循環のレベルを３以下、−３以上とする。

【0100】

ここで、血管径に基づいたレベル付けについて詳述する。
図１６は、冠循環へのレベル付けの一例を示す図である。例えば冠循環モデル作成部１２０は、血管の分岐と分岐の間を１セグメントとする。そして冠循環モデル作成部１２０は、すべてのセグメントに、そのセグメントの径を基にしたレベル付けを行う。図１６の例では、同レベルのセグメントの連なりを１つのエレメントとしている。

【0101】

冠循環モデル作成部１２０は、毛細血管から順に、太い血管に向かって、各セグメントへのレベル付けを行う。例えば冠循環モデル作成部１２０は、径が８μｍ以下のセグメントは毛細血管として扱う。毛細血管のレベルは０とする。冠循環モデル作成部１２０は、毛細血管以外のセグメントのレベル付け、樹状構造の幹へと向かって行う。例えば冠循環モデル作成部１２０は、レベルＮ（Ｎは、０以上の整数）同士の血管が合流した場合、上流の血管のレベルは１大きくし、レベルＮ＋１とする。また冠循環モデル作成部１２０は、異なるレベルのセグメントが合流した場合、合流位置から上流のセグメントのレベルは、合流した下流のセグメントのうちの大きい方のレベルに合わせる。例えば、レベルＮ＋１のセグメントとレベルＮのセグメントとが合流した場合、合流位置から上流のセグメントのレベルは、Ｎ＋１となる。また冠循環モデル作成部１２０は、以上の規則を原則として、セグメントの径に応じて修正を加える。

【0102】

さらに冠循環モデル作成部１２０は、計測結果に基づいて、以下の項目についてレベルごとに平均と標準偏差を求める。
・血管径
・血管長
・segment-to-element ratio（ＳＮ）
・コネクティビティマトリクス（ＥＲ）
ＳＮは同じレベルのセグメントが連続してつながる数を表す。コネクティビティマトリクスは行列の形で与えられ、ＥＲ（Ｍ、Ｎ）はレベルＭ（Ｍは、０以上の整数）の血管節点からレベルＮの血管が分岐する頻度を表している。毛細血管にはレベル０が割り当てられているが、冠循環モデル作成部１２０は、毛細血管をさらにC0a、C00、C0v、Cccの４種に分類する。C0aは、動脈末端と接続している毛細血管である。C0vは、静脈末端と接続している毛細血管である。C00は、心筋と平行に走行している毛細血管である。Cccは、心筋の長軸方向に対して縦方向に走行している毛細血管である。

【0103】

＜＜マクロモデル作成＞＞
冠循環モデル作成部１２０は、血管の分類が完了すると、レベル１１〜６の動脈と、レベル−６〜−１２の静脈とに基づいて、マクロモデルを作成する。

【0104】

図１７は、マクロモデル作成処理の手順の一例を示すフローチャートである。以下、図１７に示す処理をステップ番号に沿って説明する。
［ステップＳ１１１］冠循環モデル作成部１２０は、初期の心表面冠血管モデルを作成する。例えば冠循環モデル作成部１２０は、人間の心臓を撮影した断層写真を基にして、初期の心表面冠血管データを作成する。次に冠循環モデル作成部１２０は、心表面冠血管データを、心臓有限要素モデルの左心室外壁、右心室外壁に射影する。これにより、初期の心表面冠血管モデルが作成できる。

【0105】

［ステップＳ１１２］冠循環モデル作成部１２０は、初期の心表面冠血管モデルに対して、径、レベル、各節点での分岐の有無といった情報を付与する。冠循環モデル作成部１２０は、作成した心表面冠血管モデルを画面に表示できる。

【0106】

図１８は、心表面冠血管モデルの表示例を示す図である。図１８の例では、心表面冠血管モデル表示画面７０内の左側には、動脈・静脈別画像７１が表示されている。動脈・静脈別画像７１は、動脈と静脈とを色分けして表示したものである。心表面冠血管モデル表示画面７０内の右側には、血管径別画像７２が表示されている。血管径別画像７２は、血管の径に応じて色分けして表示したものである。

【0107】

図１７の説明に戻る。
［ステップＳ１１３］冠循環モデル作成部１２０は、心表面冠血管モデルに中隔枝を追加する。中隔枝とは、心臓の血管の一部に付けられた名称であり、顕微鏡写真などで観察可能である。

【0108】

図１９は、中隔枝を付与した心表面冠循環モデルの一例を示す図である。図１９に示すように、顕微鏡写真や解剖図などには、動脈や静脈の中隔枝がどのように出ているのかが表されている。冠循環モデル作成部１２０は、このような解剖図などに基づいて、中隔枝を付与した心表面冠循環モデル７３を生成する。

【0109】

図１７の説明に戻る。
［ステップＳ１１４］冠循環モデル作成部１２０は、計算データと乱数とを用い、レベル６以上、−６以下の血管ネットワークを生成する。

【0110】

まず、冠循環モデル作成部１２０は、初期の心表面冠血管の分岐節点から派生する血管のレベルを決定する。例えば、冠循環モデル作成部１２０は、分岐元の血管レベルをＮとしたときコネクティビティマトリクスＥＲ（Ｎ、ｘ）と０から１の間の乱数とを用いて、派生する下位の血管のレベルを決定する。冠循環モデル作成部１２０は、分岐元の血管が末端部の場合は、２つの下位の血管のレベルを決定する。次に冠循環モデル作成部１２０は、各レベルのＳＮの分布が、計測結果の平均値と標準偏差が規定する正規分布に従うと仮定し、仮定に従った確率で乱数を生成する。そして冠循環モデル作成部１２０は、得られた乱数の値を、小数点以下を四捨五入することで下位の血管エレメントを構成するセグメント数を決定する。

【0111】

また冠循環モデル作成部１２０は、血管径と長さも、ＳＮ同様、レベルごとの分布が計測結果の平均値と標準偏差が規定する正規分布に従うと仮定する。そして冠循環モデル作成部１２０は、血管径について、１エレメント内では一定として、仮定に従った確率で生成した乱数を使って決定する。また冠循環モデル作成部１２０は、血管長について、セグメントごとにそれぞれ、仮定に従った確率で生成した乱数を使って求める。このようにして、血管ネットワークが生成される。

【0112】

［ステップＳ１１５］最後に冠循環モデル作成部１２０は、生成した血管ネットワークを、心臓有限要素メッシュへマッピングする。その結果、図２に示したようなマクロモデル１０が生成される。

【0113】

＜＜ミクロモデル作成＞＞
次に冠循環モデル作成部１２０は、ミクロモデルを作成する。
図２０は、ミクロモデル作成処理の手順の一例を示すフローチャートである。以下、図２０に示す処理をステップ番号に沿って説明する。

【0114】

まず冠循環モデル作成部１２０は、ステップＳ１２１〜Ｓ１２４において、微小循環モデルの作成処理を行う。微小循環モデルは、動脈側のレベル３〜１の細血管、毛細血管（レベル０）、および静脈側のレベル−３〜−１の細血管を、合わせて表現したものである。太い冠血管を示すマクロモデルが２分岐を基本としたツリー状の分岐構造となっているのに対して、微小循環モデルはネットワーク構造となっていることが大きな特徴である。微小循環系は、例えば細動脈１本と細静脈２本の間に形成され、長軸方向はおおよそ５００μｍ程度のサイズである。そこで冠循環モデル作成部１２０は、ツリー状の血管分岐作成方法に加え、 以下の３つの拘束条件の基に微小循環モデルのネットワークを作成する。これにより、実際の血管形態が再現される。
・条件１：１本の細動脈（レベル３）から派生した毛細血管はやがて２本の細静脈（レベル−３）へと合流していく。
・条件２：毛細血管の長軸方向の長さは５００μｍ 程度である。
・条件３：毛細血管節点は横のつながりを持つ場合がある。

【0115】

冠循環モデル作成部１２０は、微小循環モデルを、マクロモデルの生成で述べた計測データと乱数を用いた方法により作成する。具体的には冠循環モデル作成部１２０は、ステップＳ１２１〜Ｓ１２３の３ステップで微小循環モデルを作成する。

【0116】

［ステップＳ１２１］冠循環モデル作成部１２０は、動脈側のレベル３と、静脈側のレベル−３との血管ネットワークを作成する。例えば冠循環モデル作成部１２０は、動脈側について、Kassabらの計測データに従い、１本のレベル３の細動脈からレベル１の末端の細動脈までツリー状の分岐モデルを作成する。また冠循環モデル作成部１２０は、動脈側について、動脈側同様、２本のレベル−３の細静脈からレベル−１までツリー状の分岐モデルを作成する。

【0117】

［ステップＳ１２２］冠循環モデル作成部１２０は、心筋と平行に走行する毛細血管を作成する。例えば冠循環モデル作成部１２０は、動脈側の末端と静脈側の末端の距離が約５００μｍ程度になるまで動脈側の末端より毛細血管要素をつなげていく。なお冠循環モデル作成部１２０は、毛細血管要素の径、長さは、Kassabの計測データに乱数を与えることで要素ごとに決定する。また一般に動脈側の末端数と静脈側の末端数は等しくないため、冠循環モデル作成部１２０は、両者を接続するためには毛細血管節点において分岐を作成することで、末端数を合わせる。このような毛細血管レベルの分岐は、実際に観察結果に沿ったものである。なお冠循環モデル作成部１２０は、分岐の数については、動脈側の末端数と静脈側の末端数の差によって決定する。また冠循環モデル作成部１２０は、分岐を作成する毛細血管節点は、確率的に選択する。

【0118】

［ステップＳ１２３］冠循環モデル作成部１２０は、毛細血管節点に横のつながりを作成する。例えば冠循環モデル作成部１２０は、心筋細胞の長軸方向とおおよそ垂直に走行する、心筋節点同士を横方向につなげる血管要素を作成する。なお冠循環モデル作成部１２０は、このとき作成する血管要素の血管径、長さについて、心筋と平行に走行する血管の作成時と同様、乱数により決定する。また冠循環モデル作成部１２０は、横のつながりを作成する毛細血管節点については、確率的に選択する。

【0119】

このようして、微小循環モデルが生成される。
図２１は、微小循環モデルの作成過程の一例を示す図である。まずステップＳ１２１で、レベル３の動脈側の血管ネットワーク８１と、レベル−３の静脈側の血管ネットワーク８２とが作成される。

【0120】

次にステップＳ１２２で、血管ネットワーク８１と血管ネットワーク８２との間に、５００μｍ程度の、心筋と平行に走行する毛細血管が生成される。毛細血管は、毛細血管節点で分岐している。

【0121】

次にステップＳ１２３で、心筋細胞の長軸方向とおおよそ垂直に走行する、心筋節点同士を横方向につなげる血管要素が追加される。図２１に示しているように、心筋細胞は縦方向に約１００μｍ、横方向に約１０μｍの縦長の形状をしている。このような心筋細胞の１００μｍの長さを有する方向が、長軸方向である。なお、毛細血管の径は、６μｍ前後である。

【0122】

このようにして、３次元構造の微小循環モデルが完成する。図２０に戻り、ステップＳ１２４以降の処理を説明する。
［ステップＳ１２４］冠循環モデル作成部１２０は、対称性モデルを作成する。例えば冠循環モデル作成部１２０は、マクロモデルの末端から７回分岐した後、微小循環モデルと接続するような対称性モデルを作成する。動脈側の分岐の内訳を示すと、レベル５（１本）→レベル５（２本）→レベル５（４本）→レベル４（８本）→レベル４（１６本）→レベル４（３２本）→レベル３（６４本）→レベル３（１２８本）→レベル３（微小循環モデルの上端）となる。静脈側は、動脈側のレベルに負号をつけたものと同じである。

【0123】

［ステップＳ１２５］冠循環モデル作成部１２０は、微小循環モデルと対称性モデルとを接続する。ステップＳ１２４で作成した対称性モデルは、動脈側・静脈側、それぞれ１２８の末端を有している。つまりマクロモデルの末端に、対称性モデルを介して、微小循環モデルが１２８個並列に接続される。

【0124】

対称性モデルにおいて同じレベルの血管要素が３つ連続して連なっていること。そこでまた冠循環モデル作成部１２０は、血管径を補正する。例えばまた冠循環モデル作成部１２０は、各血管要素のレベルの平均値と標準偏差を用いて、各血管の血管径を、上流より平均値＋σ、平均値、平均値−σとする。なお冠循環モデル作成部１２０は、血管要素の長さは平均値とする。

【0125】

図２２は、対称性モデルの末端に微小循環モデルを１２８個並列接続したミクロモデルの一例を示す図である。図２２に示すように、ミクロモデル７４の対称性モデルの部分は、動脈系・静脈系、それぞれ７回ずつ分岐している。そして対称性モデルの末端に、１２８個の微小循環モデルが接続されている。図２２の右側には、対称性モデルの血管のレベル、長さ、直径、および対称な血管の数を示している。なお長さと直径の単位は、μｍである。

【0126】

以上のようにして、冠循環ネットワークモデルが作成される。
＜シミュレーション実行＞
次に、シミュレーションの実行方法について詳細に説明する。以下、まず心筋の圧力節点ごとのミクロモデルの数の設定方法、心筋と微小循環の連成手法を順に説明し、最後にマルチスケール解析の計算アルゴリズムを説明する。

【0127】

＜＜ミクロモデルの数の設定＞＞
シミュレーション条件決定部１３０は、心筋容積に対する細血管の密度の貫壁性による違いを考慮するため、心外膜側で１、心内膜側で２．２となる血管密度比を、心筋節点に与える。

【0128】

図２３は、血管密度比の設定結果を表した図である。図２３は、心臓の断面を、血管密度比に応じて色分けして表示したものである。シミュレーション条件決定部１３０は、このような血管密度比を用いて、各圧力節点のミクロモデルの個数を決定する。

【0129】

例えば実行指示部２１０は、ミクロモデルの血管の総容積を心臓容積の８％とし、心筋節点に単位体積当たりに含まれるミクロモデル血管容積ρ_kを求める。そしてシミュレーション条件決定部１３０は、第１の実施の形態と同様に、式（１）を用いて、圧力節点ｋに接続されるミクロモデルの数ｎ_kを計算する。

【0130】

さらに、シミュレーション条件決定部１３０は、ユーザからの指示に従い、シミュレーション時間の進行間隔、シミュレーションの終了条件Newton-Raphson法による解析の収束条件などを決定する。そしてＣＡＤシステム１００から管理ノード２００へ、マクロモデル、ミクロモデル、およびシミュレーション条件を送信する。すると管理ノード２００から複数の計算ノードに冠循環シミュレーションのジョブが投入され、複数の計算ノードの並列処理により冠循環シミュレーションが実行される。

【0131】

＜＜ミクロモデルの血管に関する式＞＞
ミクロモデルの血管は心筋に、直接、酸素や栄養を供給する役割を担っている。そこで冠循環シミュレーションでは、血液の輸送に加え血管体積変化による血液の貯蓄や放出を考慮する。一般に血管要素の２端点で圧力値が異なり、要素内でも径が異なる。

【0132】

図２４は、血管径の違いを示す図である。図２４に示すように、コンダクタンスｋ_ijの血管要素を、コンダクタンスｋ_i（ｔ）の分割要素とコンダクタンスｋ_j（ｔ）の分割要素との２つに分けて考える。分割要素の内圧はμ_i（ｔ），μ_j（ｔ）で一定と考えると、血管径Ｄ_i（ｔ）は内圧（冠血圧）と外圧ｐ_m（ｔ）の差により式（１１）のように決定される。なお、外圧とは心筋圧であり、心筋有限要素の圧力節点の値に相当する。αは血管の弾性を示し、＊は基準状態を示す。

【0133】

【数11】

【0134】

分割要素のコンダクタンスｋは便宜的にポアズイユ流れより得られる値で代用する。よって血管要素のコンダクタンスｋ_ij（ｔ）は分割要素のコンダクタンスｋ_i（ｔ），ｋ_j（ｔ）の調和平均として式（１２）のように求められる。

【0135】

【数12】

【0136】

また、分割要素が溜め込む血液は体積の時間微分で表すことができ（式（１３）第二項）、結局マクロモデルの血管節点ｉにおいては次式の流量保存則が成立する。

【0137】

【数13】

【0138】

＜＜心臓有限要素モデルと血液からなる混合体の基礎式＞＞
心筋と微小循環系に含まれる血液からなる混合体の運動方程式は、次の式（１４）で与えられる。

【0139】

【数14】

【0140】

ｐは心筋内圧である。ρは混合体密度である。ａは加速度である。σ^sはCauchy応力テンソルである。体積変化率Ｊを用いた混合体の体積保存則は、次の式（１５）となる。

【0141】

【数15】

【0142】

式（１５）の第二項がミクロモデルの体積変化であり、ミクロモデルに蓄えられる血液が増加した分だけ混合体の容積が増加することを表している。ΔＶ_Microはミクロモデル１個の初期状態からの容積増加量である。ｎは単位心筋体積当たりに含まれるミクロモデルの数である。

【0143】

結局、未知数は心筋の変位３成分、心筋圧力１成分、血管圧力１成分の計５種となり、以下の式（１６）で表される。

【0144】

【数16】

【0145】

Ａ_mic：ミクロモデルのマトリクス
Ｂ_mac：マクロモデルのマトリクス
Ｃ：ミクロモデルとマクロモデルの相互作用
Ｄ：心筋および流体節点の変位に関連するマトリクス
Ｆ：:心筋および流体節点の圧力に関連するマトリクス
Ｋ：ミクロモデルと心筋節点の圧力の相互作用
Ｈ：マクロモデルと流体節点の圧力の相互作用
Δμ_mic：ミクロモデルの節点の圧力に関連するベクトル
Δμ_mac：マクロモデルの節点の圧力に関連するベクトル
Δｕ：心筋および流体節点の変位に関するベクトル
Δｐ：心筋および流体節点の圧力に関するベクトル
ｒ：残差ベクトル
未知数は、以上の式（１６）を解くことで求められる。

【0146】

＜＜シミュレーション手順＞＞
計算ノード４００，５００，６００，・・・は、管理ノード２００の実行指示部２１０の指示に応じ、心臓の冠循環シミュレーションを並列処理で実行する。例えば、ＭＰＩ（Message-Passing Interface）を用いた並列処理が行われる。その際、実行指示部２１０は、マクロモデルの計算を担当するプロセス（macro procs）を、コアに割り当てる。マクロモデルの計算を担当するプロセスが割り当てられコアは、マクロ解析部４５０として機能する。

【0147】

また実行指示部２１０は、ミクロモデルの計算を担当するプロセス（micor procs）を、複数のコアに割り当てる。このときミクロモデルの計算を担当するプロセスを割り当てるコアには、例えば、そのコアで計算するミクロモデルの番号や、ミクロモデルが接続される圧力節点の位置が指定される。ミクロモデルの計算を担当するプロセスを割り当てられたコアは、ミクロ解析部５１０，・・・，として機能する。

【0148】

また実行指示部２１０は、シミュレーション全体を制御するプロセスを、コアに割り当てる。シミュレーション全体を制御するプロセスが割り当てられたコアは、シミュレーション制御部４４０として機能する。なお、シミュレーション制御部４４０は、例えばマクロプロセスが割り当てられたマクロ解析部４５０と同じコアで実行することもできる。

【0149】

図２５は、心臓の冠循環シミュレーションの手順の一例を示すフローチャートである。図２５に示す処理のうち、ステップＳ１３１〜Ｓ１３５は、シミュレーション制御部４４０が実行する全体制御処理である。またステップＳ１４１〜Ｓ１４５は、ミクロ解析部５１０・・・が実行するミクロプロセスの処理である。ステップＳ１５１〜Ｓ１５５は、マクロ解析部４５０が実行するマクロプロセスの処理である。なお、マクロプロセスの処理には、心筋の収縮・弛緩、内腔血液の拍出の解析も含まれている。

【0150】

以下、全体制御、ミクロプロセス、マクロプロセスに分けて、図２５に示す処理をステップ番号に沿って説明する。まず、全体制御の処理について説明する。
［ステップＳ１３１］シミュレーション制御部４４０は、心臓の有限要素モデル３１０と冠循環モデル（マクロモデル３２０とミクロモデル３３０）とを、ストレージ装置３００から各計算ノード４００，５００，６００，・・・に読み込ませる。各計算ノード４００，５００，６００，・・・では、読み込んだ情報を、メモリに記憶する。

【0151】

［ステップＳ１３２］シミュレーション制御部４４０は、シミュレーション上の時刻を示すタイムステップを、１ステップ分進める。
［ステップＳ１３３］シミュレーション制御部４４０は、Newton−Raphson法による非線形解析の反復処理を進める。例えばシミュレーション制御部４４０は、マクロ解析部４５０と複数のミクロ解析部５１０，・・・に、Newton−Raphson法による非線形解析の次の段階の処理の実行を指示する。これにより複数のミクロ解析部５１０，・・・それぞれにおいて、ステップＳ１４１〜Ｓ１４５の処理が並列で実行される。また、マクロ解析部４５０において、ステップＳ１５１〜Ｓ１５５の処理が実行される。

【0152】

［ステップＳ１３４］シミュレーション制御部４４０は、Newton−Raphson法による非線形解析の結果が収束したか否かを判断する。例えばシミュレーション制御部４４０は、前回の解析結果と今回の解析結果との差が所定の閾値以下であれば、結果が収束したと判断する。シミュレーション制御部４４０は、結果が収束した場合、処理をステップＳ１３５に進める。またシミュレーション制御部４４０は、結果が収束していなければ、処理をステップＳ１３３に進める。

【0153】

［ステップＳ１３４］シミュレーション制御部４４０は、シミュレーションを実行するすべてのタイムステップについて計算が完了したか否かを判断する。シミュレーション制御部４４０は、すべてタイムステップの計算が完了していれば、冠循環シミュレーションを終了する。またシミュレーション制御部４４０は、すべてタイムステップ計算が完了していなければ、処理をステップＳ１３２に進める。

【0154】

このようにして、冠循環シミュレーション全体が制御される。
次に、複数のミクロ解析部５１０，・・・それぞれが実行するミクロプロセスのステップＳ１４１〜Ｓ１４５について説明する。

【0155】

［ステップＳ１４１］複数のミクロ解析部５１０，・・・それぞれは、自身が解析するマクロモデルが接続されている圧力節点に応じた、ミクロモデルのマトリクスＡ_Micを組み立てる。

【0156】

［ステップＳ１４２］複数のミクロ解析部５１０，・・・それぞれは、ステップＳ１４１で組み立てたマトリクスＡ_Micに基づいて、数（５）に従ってΔμ_Mac（μはチルダ付き）を算出する。

【0157】

［ステップＳ１４３］複数のミクロ解析部５１０，・・・それぞれは、ステップＳ１４３で計算したΔμ_Mac（μはチルダ付き）を、マクロ解析部４５０に送信する。
［ステップＳ１４４］複数のミクロ解析部５１０，・・・それぞれは、マクロ解析部４５０から、自身が解析するミクロモデルが接続されている圧力節点についての圧力の増分Δμ_macを受信する。

【0158】

［ステップＳ１４５］複数のミクロ解析部５１０，・・・それぞれは、ステップＳ１４４で受信した値を用いて、式（７）の計算を行い、Δμ_Micを算出する。
このようにして、ミクロモデルごとの計算が、複数のミクロ解析部５１０による並列処理で実行される。

【0159】

次に、マクロ解析部４５０が実行するマクロプロセスのステップＳ１５１〜Ｓ１５５について説明する。
［ステップＳ１５１］マクロ解析部４５０は、ミクロモデルのマトリクスＢ_Micを組み立てる。

【0160】

［ステップＳ１５２］マクロ解析部４５０は、各ミクロ解析部５１０，・・・から、Δμ_Mac（μはチルダ付き）の値を受信する。
［ステップＳ１５３］マクロ解析部４５０は、式（４）に示したＢ_Mac（Ｂはチルダ付き）の定義式を組み立てる。

【0161】

［ステップＳ１５４］マクロ解析部４５０は、式（６）に基づいて、各圧力節点におけるΔμ_Macを算出する。
［ステップＳ１５５］マクロ解析部４５０は、複数のミクロ解析部５１０，・・・それぞれに、そのマクロ解析部が解析しているミクロモデルが接続された圧力節点のΔμ_Macを送信する。

【0162】

図２５に示した処理を、概略的に説明すると、以下のようになる。
まずシミュレーション制御部４４０が、入力データとして心臓の有限要素モデルおよび冠循環モデルを読み込み、以降は、タイムステップごとに、マクロ解析部４５０とミクロ解析部５１０，・・・とが、Newton-Raphson法で非線形解析の反復計算を実行する。反復計算において、マクロ解析部４５０は、まずマトリクスＢ_Macを組み立てる。同時に、複数のミクロ解析部５１０，・・・は、マトリクスＡ_Micを組み立て、Ａ_Micが構成する方程式を解く。次に、ミクロ解析部５１０，・・・が、マクロ解析部４５０へΔμ_Mac（μはチルダ付き）が送信する。マクロ解析部４５０は、Ｂ_Mac（Ｂはチルダ付き）を計算し、方程式を解く。マクロ解析部４５０はミクロ解析部５１０，・・・へΔμ_Maを送信する。最後にミクロ解析部５１０，・・・は、Δμ_Micを計算する。

【0163】

このように、第２の実施の形態では、ミクロ解析部５１０によるミクロプロセスを担当するＣＰＵの各コアが独立してミクロモデルの計算をすることができる。またミクロプロセスとマクロプロセスとが、一部において同時に動作すること可能となり、並列計算システムの性能を効果的に利用することが可能となる。さらに、対称性モデルの導入により、導入前に比べてミクロモデルのデータサイズが、約１２８分の１に削減されている。

【0164】

そして、並列計算システムの性能を効率的に利用したことで、毛細血管の血流までを含んだ心臓の冠循環シミュレーションが可能となっている。毛細血管の血流までを含んだ心臓の冠循環シミュレーションが実現したことで、例えば、これまで微小循環の観察には限界があった、心臓の動作の把握が可能となる。また、化合物の輸送と毛細管・筋組織間のやりとりを追加することで、医薬品開発における薬剤等の効果・副作用予測に適用可能なシミュレーションを実現することも可能となる。

【0165】

〔第３の実施の形態〕
次に、第３の実施形態について説明する。第３の実施の形態は、第１または第２の実施形態に示した血流のシミュレーションに、心筋細胞の代謝機能を解析するシミュレーションを組み合わせたものである。代謝とは心筋細胞が毛細血管内の血流から得た酸素・栄養からアデノシン三リン酸を産生する一連の化学反応である。代謝機能の数理モデルについては非特許文献５などがある。血流のシミュレーションと心筋細胞の代謝機能を解析するシミュレーションとの組み合わせは、毛細管までの血流のシミュレートが可能となったことにより、よりはじめて可能となるものである。

【0166】

第３の実施の形態では、心筋細胞の代謝機能を解析するシミュレーションと組み合わせるために、第１または第２実施形態において求解する方程式に、以下の２つの式（１７）、（１８）を追加する。

【0167】

【数17】

【0168】

式（１７）は、冠循環の血液で運ばれる酸素・栄養の濃度の方程式である。式（１７）において、Ｃは、酵素・栄養素の濃度である。ｖは、血液の流速である。Ｄは、拡散係数である。

【0169】

【数18】

【0170】

式（１８）は、血液と心筋細胞との間の酸素・栄養の濃度差に従った、酸素・栄養の移動の方程式である。式（１８）において、Ｃは、酸素・栄養素の濃度である。Ｋは、拡散係数である。ｑは、心筋細胞の代謝による消費を示す値である。

【0171】

式（１７）、（１８）を含めて、第１または第２実施形態と同様の手法で解を求める。なお式（１７）は、マクロモデル・ミクロモデルの双方の解析に適用される。式（１８）は、ミクロモデルおよびそのミクロモデルが接続された心筋圧力節点を持つ心筋有限要素に適用される。

【0172】

このようにして、血流と心筋細胞の代謝機能とを合わせて解析するシミュレーションが可能となる。
〔第４の実施の形態〕
次に、第４の実施の形態について説明する。第４の実施形態は、第３の実施形態における計算量削減手法である。

【0173】

冠循環の血液流量の影響による代謝機能の変化は、心臓１０拍動以上を経過してはじめて現れる。ところが、並列計算機の性能によっては１０拍動分の計算時間が膨大になる可能性がある。そこで、第４の実施の形態では、少ない計算量でシミュレーションを行う。なお、第４の実施の形態のシステム構成は、第２の実施の形態と同様である。そこで、以下の説明では、図８〜図１１に示した第２の実施の形態の各要素を用いて、第４の実施の形態の処理について説明する。

【0174】

図２６は、第４の実施の形態におけるシミュレーション手順の一例を示すフローチャートである。以下、図２６に示す処理をステップ番号に沿って説明する。
［ステップＳ２０１］シミュレーション制御部４４０は、心臓の有限要素モデル３１０、冠循環モデル（マクロモデル３２０とミクロモデル３３０）、および心筋細胞の代謝機能の解析に必要なデータを、ストレージ装置３００から各計算ノード４００，５００，６００，・・・に読み込ませる。各計算ノード４００，５００，６００，・・・では、読み込んだ情報を、メモリに記憶する。

【0175】

［ステップＳ２０２］シミュレーション制御部４４０は、シミュレーション上の時刻を示すタイムステップを、１ステップ分進める。
［ステップＳ２０３］マクロ解析部４５０と複数のミクロ解析部５１０，・・・とは、連携して、１拍動分の計算を行う。このステップＳ２０３では、心筋の運動、内腔血液の拍出、冠循環の流量、酸素・栄養の輸送、および代謝機能が計算される。そしてマクロ解析部４５０と複数のミクロ解析部５１０，・・・とは、冠循環流量のデータをメモリなどに保存する。

【0176】

［ステップＳ２０４］マクロ解析部４５０と複数のミクロ解析部５１０，・・・とは、連携して、１０拍動分の計算を行う。このステップＳ２０４では、酸素・栄養の輸送と代謝機能とが計算される。

【0177】

［ステップＳ２０５］シミュレーション制御部４４０は、シミュレーションを実行するすべてのタイムステップについて計算が完了したか否かを判断する。シミュレーション制御部４４０は、すべてタイムステップの計算が完了していれば、冠循環シミュレーションを終了する。またシミュレーション制御部４４０は、すべてタイムステップ計算が完了していなければ、処理をステップＳ２０２に進める。

【0178】

このように第４の実施の形態では、１拍動分のシミュレーションをまず行い、心筋の運動、内腔血液の拍出、冠循環の流量、酸素・栄養の輸送、代謝機能を解析する。このとき１拍動における各時刻の冠循環の流量データは保存しておく。次に、保存した冠循環の流量データを用いて酸素・栄養素の輸送と代謝機能のみの計算を１０拍動分行う。これにより代謝機能の変化を解析するのに必要な１０拍動以上のシミュレーションの計算量を削減することができる。その結果、例えば、並列計算機の性能が制限された環境下でも、実用に耐え得る時間でシミュレーションを行うことができる。

【0179】

〔第５の実施の形態〕
次に、第５の実施の形態について説明する。第５の実施形態は、第１〜４の実施形態を基にした冠動脈狭窄時の心臓に対する影響を評価するものである。

【0180】

図２７は、第５の実施の形態におけるシステム構成の一例を示す図である。図２７において、第２の実施の形態と同じ要素には同じ符号を付し、説明を省略する。第５の実施の形態では、図８に示した第２の実施の形態のシステムに対して、ストレージ装置７００、ネットワーク９１、および端末装置８００が追加されている。

【0181】

ストレージ装置７００は、ネットワーク６２に接続されている。そしてストレージ装置７００は、計算ノード４００，５００，６００，・・・により、解析結果保存用に使用される。端末装置８００は、ネットワーク９１を介してストレージ装置７００に接続されている。端末装置８００は、ストレージ装置７００に格納された解析結果に基づいて冠動脈狭窄時の心臓に対する影響を評価する。

【0182】

例えばＣＡＤシステム１００を用い、ＣＴ（Computed Tomography）／ＭＲＩ（Magnetic Resonance Imaging）といった医療画像から得られる、狭窄した冠動脈を反映した冠循環モデルのデータを作成する。或いはＣＡＤシステム１００を用い、ある患者の、将来狭窄が起きる可能性のある冠動脈の部位に対して、適切な狭窄を設定・反映した冠循環モデルのデータを作成する。このようなデータにより、計算ノード４００，５００，６００，・・・が、第１〜４の実施形態で説明したシミュレーションを行い、結果保存用のストレージ装置７００にシミュレーション結果を保存する。ユーザは、評価用の端末装置８００により、解析結果を示すデータを取得する。そして、端末装置８００が、心筋の運動、内腔血液の拍出、冠循環の流量、酸素・栄養の輸送、代謝機能について、冠動脈狭窄の影響を評価する。

【0183】

〔その他の実施の形態〕
上記の第１〜第５の実施の形態では、心臓の冠循環シミュレーションの例を用いているが、同様の血行動態シミュレーションを、別の器官に関して行うこともできる。

【0184】

また第２〜第５の実施の形態のシステムは、複数の装置で構成されているが、任意の２つ以上の装置を、１つの装置にまとめることもできる。例えば、ＣＡＤシステム１００と管理ノード２００との機能を、１つのコンピュータにまとめることができる。

【0185】

さらに第２〜第５の実施の形態では、ＣＰＵのコアがプログラムを実行することによって処理を実現するものとしたが、プログラムで記述された処理の一部を、電子回路に置き換えることが可能である。例えば、上記の処理機能の少なくとも一部を、ＤＳＰ（Digital Signal Processor）、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。

【0186】

以上、実施の形態を例示したが、実施の形態で示した各部の構成は同様の機能を有する他のものに置換することができる。また、他の任意の構成物や工程が付加されてもよい。さらに、前述した実施の形態のうちの任意の２以上の構成（特徴）を組み合わせたものであってもよい。

【符号の説明】

【0187】

１ 記憶手段
１ａ 形状モデル
１ｂ 第１の血管ネットワークモデル
１ｃ 第２の血管ネットワークモデル
２，３，４，５，・・・ 演算装置
２ａ 第１の解析部
３ａ，４ａ，５ａ，・・・ 第２の解析部
Ｘ シミュレーション装置

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】

【図20】

【図21】

【図22】

【図23】

【図24】

【図25】

【図26】

【図27】