JP7460907B2 - シミュレーションプログラム、シミュレーション方法およびシミュレーションシステム - Google Patents

Description

本発明はシミュレーションプログラム、シミュレーション方法およびシミュレーションシステムに関する。

近年、コンピュータの計算能力の向上に伴い、精密なモデルを利用したコンピュータシミュレーションが行われるようになっている。精密なモデルとして、対象物が存在する空間を多数の小領域に分割した有限要素モデルを用いることがある。シミュレーションとして、流体の挙動を解析する流体シミュレーションが行われることがある。また、心臓の拍動のように、物体の周期的な運動のシミュレーションが行われることがある。周期的な運動のシミュレーションでは、物体の状態を時系列に並べた場合に、同様の状態変化のパターンが時間軸上に繰り返し出現することになる。

なお、トランジスタの電気的特性のシミュレーションを行う設計支援装置が提案されている。提案の設計支援装置は、モデルのパラメータとして各節点に仮の不純物濃度を設定し、次元縮退によって各節点の表面ポテンシャルを算出し、表面ポテンシャルから推定される電気的特性が所定条件を満たす場合に、パラメータの値を確定する。

また、血管の形状を示す三次元画像データと、血液の流量や圧力損失を示す流体データとを用いて流体解析を行い、血管の血行状態を示す指標値を算出する画像処理装置が提案されている。また、有限要素法の非定常解析により、金属材料の加工における材料流れのシミュレーションを行う解析システムが提案されている。

特開２０１０－２８７６１４号公報 特開２０１７－７０７４２号公報 特開２０１９－１２８６２１号公報

シミュレーションの１つの方法として、主要な部位を有限要素モデルなどの精密なモデルで表現し、その周囲の部位を他のモデルで表現し、主要モデルと周囲モデルとを結合してシミュレーションを行う方法が考えられる。例えば、周囲部位を主要部位よりも簡易なモデルで表現することが考えられる。

その場合に、主要モデルと周囲モデルとの間の相互作用を計算できるようにするため、シミュレーションの開始時に、周囲モデルの変数に初期値を代入することがある。ここで、変数の初期値をどの様に決定すればよいかが問題となる。

周期的な運動のシミュレーションでは、周囲モデルの変数は、与えられた初期値から、相互作用の計算の繰り返しを通じて、一定の数値となる収束値へ収束する。しかし、変数の初期値が収束値から外れているほど、収束に時間を要する。このため、シミュレーションを行う周期数が多くなり、計算量が増大するおそれがある。これに対し、理想的な初期値をユーザが事前に見積もることは容易でない。

１つの側面では、本発明は、周期的な運動のシミュレーションの計算量を削減できるシミュレーションプログラム、シミュレーション方法およびシミュレーションシステムを提供することを目的とする。

１つの態様では、コンピュータに以下の処理を実行させるシミュレーションプログラムが提供される。周期的な運動を行う第１の部位を示す第１のモデルと、第１の部位に接続されており周期的な運動の影響を受ける第２の部位を示す第２のモデルとを取得する。第２のモデルに含まれる変数の初期値として第１の状態値を割り当て、第１のモデルおよび第２のモデルと第１の状態値とを用いて、周期的な運動のシミュレーションを所定周期数分実行して、所定周期数分の特徴量を算出する。所定周期数分の特徴量に基づいて、第１のモデルより変数が少なく第１のモデルと置換可能な第３のモデルを生成する。第３のモデルおよび第２のモデルを用いて、周期的な運動のシミュレーションを所定の収束条件が満たされるまで継続して実行し、所定の収束条件が満たされた周期において第２のモデルに含まれる変数がもつ第２の状態値を抽出する。

また、１つの態様では、コンピュータが実行するシミュレーション方法が提供される。また、１つの態様では、シミュレーションシステムが提供される。

１つの側面では、周期的な運動のシミュレーションの計算量を削減できる。

第１の実施の形態のシミュレーションシステムを説明するための図である。 第２の実施の形態のシミュレーション装置の例を示すブロック図である。 分散処理システムの例を示す図である。 血行動態シミュレーションに用いるモデルの例を示す図である。 有限要素モデルを説明するための図である。 心室容積と心室圧のシミュレーション例を示すグラフである。 血液量のシミュレーションの収束例を示すグラフである。 心室容積の仮のシミュレーション例を示すグラフである。 心室圧の仮のシミュレーション例を示すグラフである。 三次元モデルの簡略化例を示す図である。 フーリエ係数テーブルの例を示す図である。 シミュレーション装置の機能例を示すブロック図である。 節点テーブルとメッシュテーブルの例を示す図である。 パラメータテーブルの例を示す図である。 時系列データテーブルの例を示す図である。 シミュレーションの手順例を示すフローチャートである。 シミュレーションの手順例を示すフローチャート（続き）である。

以下、本実施の形態を図面を参照して説明する。
［第１の実施の形態］
第１の実施の形態を説明する。

図１は、第１の実施の形態のシミュレーションシステムを説明するための図である。
第１の実施の形態のシミュレーションシステム１０は、心臓の拍動のような周期的な運動のシミュレーションを行う。シミュレーションシステム１０は、単一の筐体のクライアントコンピュータまたはサーバコンピュータであってもよいし、複数の筐体のサーバコンピュータを含んでもよい。シミュレーションシステム１０は、以下に説明する処理の一部または全部を並列的に実行する分散処理システムであってもよい。また、シミュレーションシステム１０は、クラウドシステムやデータセンタの計算資源であってもよい。シミュレーションシステム１０を、情報処理システムと言うこともできる。

シミュレーションシステム１０は、記憶装置１１および演算装置１２を有する。記憶装置１１は、記憶装置１１は、ＲＡＭ（Random Access Memory）やＧＰＵ（Graphics Processing Unit）メモリなどの揮発性半導体メモリでもよいし、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性ストレージでもよい。

演算装置１２は、例えば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）、ＤＳＰ（Digital Signal Processor）などのプロセッサである。ただし、演算装置１２は、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などの特定用途の電子回路を含んでもよい。プロセッサは、メモリ（記憶装置１１でもよい）に記憶されたプログラムを実行する。演算装置１２が、同一種類または異なる種類のプロセッサの集合（例えば、ＣＰＵの集合、ＧＰＵの集合、ＣＰＵとＧＰＵの集合など）であってもよい。

記憶装置１１は、モデル１３（第１のモデル）およびモデル１４（第２のモデル）を記憶する。モデル１３は、周期的な運動を行う部位（第１の部位）を示す。例えば、モデル１３は、一定周期で拍動する心臓を示す。モデル１３は、複数の節点とそれら複数の節点の間を接続する複数のエッジとを含み、各節点に変数が割り当てられる有限要素モデルであってもよい。また、モデル１３は、変数が多い複雑なモデルでもよい。

モデル１４は、モデル１３が示す部位に接続されており周期的な運動の影響を受ける他の部位（第２の部位）を示す。例えば、モデル１４は、大動脈、大静脈、肺動脈、肺静脈などの心臓の外部の血管を示す。第１の実施の形態で行うシミュレーションが、モデル１３とモデル１４との間の流体の循環のシミュレーションであってもよい。モデル１４は、抵抗やキャパシタ（コンデンサ）などの電気回路の素子を含む電気回路モデルであってもよい。また、モデル１４は、モデル１３より変数が少ない簡易なモデルでもよい。モデル１４には、予めパラメータが設定されていることがある。例えば、抵抗の抵抗値（レジスタンス）やキャパシタのキャパシタンス（静電容量）が、パラメータとして設定される。

演算装置１２は、モデル１３とモデル１４を組み合わせて、周期的な運動のシミュレーションを行う。シミュレーションを開始する際、演算装置１２は、モデル１４に含まれる変数に初期値を割り当てる。モデル１４の変数の初期値は、例えば、キャパシタンスの電荷量の初期値である。ここで、変数の初期値が不適当であっても、周期的な運動のシミュレーションを継続して実行すれば、変数の値が適切な値に収束していき、最終的には妥当なシミュレーション結果が得られる。しかし、変数の初期値が不適当であると、収束までの時間が長くなる。よって、シミュレーションを行う周期数（サイクル数）が多くなり、シミュレーションに長時間を要するおそれがある。そこで、演算装置１２は、以下のようにしてモデル１４の変数の初期値を決定する。

まず、演算装置１２は、モデル１４の変数の初期値として、状態値１６ａ（第１の状態値）を仮に割り当てる。状態値１６ａは、例えば、シミュレーション開始時点でキャパシタがもつ電荷量を示す。状態値１６ａは、ランダムな値でもよいし、ユーザが指定した値でもよいし、所定の固定値でもよい。演算装置１２は、モデル１３とモデル１４とを結合したモデルに状態値１６ａを適用して、周期的な運動のシミュレーションを実行する。

このとき、演算装置１２は、十分に少ない所定の周期数（例えば、２～３周期）分のシミュレーションを実行すればよく、所定の周期数でシミュレーションを打ち切ってよい。モデル１４の変数の初期値が不適当である場合、所定の周期数では、まだ変数の値が収束しておらず、シミュレーション結果が安定していないことがある。

演算装置１２は、この所定の周期数分のシミュレーションの結果から、モデル１３の状態に関する指標の特徴量１７を算出する。特徴量１７は、例えば、モデル１３が示す部位の容積の時間変化を示す時系列データや、当該部位の内部圧力の時間変化を示す時系列データなどである。特徴量１７は、周期的な運動のシミュレーションによって直接的に算出される場合もあるし、シミュレーション結果を変換することで得られることもある。

特徴量１７が得られると、演算装置１２は、特徴量１７に基づいて、モデル１３を近似するモデル１５を生成する。モデル１５は、モデル１３より変数が少ない簡易モデルであり、例えば、電気回路モデルである。モデル１５の生成では、例えば、演算装置１２は、モデル１５とモデル１４とを結合したモデルを用いて、周期的な運動のシミュレーションを所定周期分実行する。そして、演算装置１２は、そのシミュレーションの結果が特徴量１７に近付くように、モデル１５に含まれるパラメータの値を修正する。これは、特徴量１７を教師データとして用いて、モデル１５のパラメータを最適化していると言える。

モデル１５が生成されると、演算装置１２は、モデル１５とモデル１４とを結合したモデルを用いて、周期的な運動のシミュレーションを実行する。モデル１４の変数の初期値として、状態値１６ａを用いてもよいし他の状態値を用いてもよい。ここで、演算装置１２は、所定の収束条件を満たすまで当該シミュレーションを継続して実行する。このシミュレーションの周期数は、通常、特徴量１７を算出した際の周期数よりも多くなる。収束条件は、例えば、最新の周期の特徴量と１つ前の周期の特徴量との間の誤差が、閾値より小さいことである。また、収束条件は、例えば、最新の周期の特徴量と１つ前の周期の特徴量との間の相関係数が、閾値より大きいことである。

そして、演算装置１２は、収束条件が満たされた周期において、モデル１４に含まれる変数がもつ状態値１６ｂ（第２の状態値）を抽出する。演算装置１２は、状態値１６ｂを記憶装置１１に保存してもよいし、表示装置に表示してもよいし、他のシステムに転送してもよい。この状態値１６ｂは、シミュレーション結果が収束して安定した状態におけるモデル１４の変数の値であり、変数の理想的な初期値と推定される。例えば、演算装置１２は、モデル１３とモデル１４とを結合したモデルに状態値１６ｂを適用し、周期的な運動のシミュレーションを収束条件が満たされるまで継続して実行する。このときの所要周期数は、状態値１６ａを用いた場合よりも少なくなると期待される。

第１の実施の形態のシミュレーションシステム１０によれば、モデル１４の変数の初期値として状態値１６ａが割り当てられ、モデル１３，１４を用いて所定周期数分のシミュレーションが実行され、特徴量１７が算出される。特徴量１７に基づいて、モデル１３より変数が少ないモデル１５が生成される。そして、モデル１４，１５を用いて所定の収束条件が満たされるまでシミュレーションが実行され、収束条件が満たされた周期におけるモデル１４の変数から状態値１６ｂが抽出される。

これにより、モデル１３に接続されるモデル１４の変数に、収束後の値に近似する適切な初期値を代入することが可能となる。よって、モデル１３，１４を用いた周期的な運動のシミュレーションが収束するまでの所要周期数を少なくすることができる。そのため、シミュレーションの計算量を削減でき、所要時間を短縮することができる。

［第２の実施の形態］
次に、第２の実施の形態を説明する。
第２の実施の形態のシミュレーション装置は、人体の血行動態のシミュレーションを行う。血行動態シミュレーションは、心臓を中心とする血流のシミュレーションである。血行動態シミュレーションを、循環動態シミュレーションなどと言うこともある。

図２は、第２の実施の形態のシミュレーション装置の例を示すブロック図である。
シミュレーション装置１００は、ＣＰＵ１０１、ＲＡＭ１０２、ＨＤＤ１０３、画像インタフェース１０４、入力インタフェース１０５、媒体リーダ１０６および通信インタフェース１０７を有する。シミュレーション装置１００が有するこれらのユニットは、バスに接続されている。ＣＰＵ１０１は、第１の実施の形態の演算装置１２に対応する。ＲＡＭ１０２またはＨＤＤ１０３は、第１の実施の形態の記憶装置１１に対応する。

ＣＰＵ１０１は、プログラムの命令を実行するプロセッサである。ＣＰＵ１０１は、ＨＤＤ１０３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ１０２にロードし、プログラムを実行する。ＣＰＵ１０１は複数のプロセッサコアを備えてもよく、シミュレーション装置１００は複数のプロセッサを備えてもよい。複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」と言うことがある。

ＲＡＭ１０２は、ＣＰＵ１０１が実行するプログラムやＣＰＵ１０１が演算に使用するデータを一時的に記憶する揮発性半導体メモリである。シミュレーション装置１００は、ＲＡＭ以外の種類のメモリを備えてもよく、複数のメモリを備えてもよい。

ＨＤＤ１０３は、ＯＳ（Operating System）やミドルウェアやアプリケーションソフトウェアなどのソフトウェアのプログラム、および、データを記憶する不揮発性ストレージである。シミュレーション装置１００は、フラッシュメモリやＳＳＤ（Solid State Drive）など他の種類のストレージを備えてもよく、複数のストレージを備えてもよい。

画像インタフェース１０４は、ＣＰＵ１０１からの命令に従って、シミュレーション装置１００に接続された表示装置１１１に画像を出力する。表示装置１１１として、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイ、プロジェクタなど、任意の種類の表示装置を使用することができる。シミュレーション装置１００に、プリンタなど表示装置１１１以外の出力デバイスが接続されてもよい。

入力インタフェース１０５は、シミュレーション装置１００に接続された入力デバイス１１２から入力信号を受け付ける。入力デバイス１１２として、マウス、タッチパネル、タッチパッド、キーボードなど、任意の種類の入力デバイスを使用することができる。シミュレーション装置１００に複数種類の入力デバイスが接続されてもよい。

媒体リーダ１０６は、記録媒体１１３に記録されたプログラムやデータを読み取る読み取り装置である。記録媒体１１３として、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤなどの磁気ディスク、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）などの光ディスク、半導体メモリなど、任意の種類の記録媒体を使用することができる。媒体リーダ１０６は、例えば、記録媒体１１３から読み取ったプログラムやデータを、ＲＡＭ１０２やＨＤＤ１０３などの他の記録媒体にコピーする。読み取られたプログラムは、例えば、ＣＰＵ１０１によって実行される。なお、記録媒体１１３は可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体１１３やＨＤＤ１０３を、コンピュータ読み取り可能な記録媒体と言うことがある。

通信インタフェース１０７は、ネットワーク１１４に接続され、ネットワーク１１４を介して他の情報処理装置と通信する。通信インタフェース１０７は、スイッチやルータなどの有線通信装置に接続される有線通信インタフェースでもよいし、基地局やアクセスポイントなどの無線通信装置に接続される無線通信インタフェースでもよい。

ここで、第２の実施の形態では説明を簡単にするため、シミュレーション装置１００の処理をＣＰＵ１０１が実行することとしている。ただし、後述するシミュレーション装置１００の処理を分散処理システムによって実行することも可能である。

図３は、分散処理システムの例を示す図である。
分散処理システムは、計算ノード３１～３４を含む複数の計算ノード、管理ノード３５およびストレージノード３６を有する。計算ノード３１～３４は、複数のプロセスまたは複数のスレッドを並列に実行する。例えば、計算ノード３１～３４は、大規模連立方程式の解を反復法によって求める求解演算が複数のスレッドに分割された場合における、当該複数のスレッドを並列に実行する。管理ノード３５は、計算ノード３１～３４にプロセスやスレッドを割り当て、並列処理を制御する。ストレージノード３６は、計算ノード３１～３４によって使用されるデータを記憶する。計算ノード３１～３４、管理ノード３５およびストレージノード３６は、例えば、ネットワークに接続されている。

計算ノード３１～３４は、同一筐体内の複数のＣＰＵでもよいし、同一筐体内の複数のＧＰＵでもよい。計算ノード３１～３４それぞれがＲＡＭまたはＧＰＵメモリを備えてもよいし、計算ノード３１～３４が同一のＲＡＭまたはＧＰＵメモリを共有してもよい。また、計算ノード３１～３４が、異なる筐体のサーバコンピュータであってもよい。管理ノード３５は、計算ノード３１～３４と同一筐体内のＣＰＵであってもよく、ＲＡＭを備えていてもよい。また、管理ノード３５が、計算ノード３１～３４と異なる筐体のサーバコンピュータまたはクライアントコンピュータでもよい。ストレージノード３６は、ＨＤＤやＳＳＤなどの不揮発性ストレージでもよい。また、ストレージノード３６が、計算ノード３１～３４と異なる筐体のストレージサーバであってもよい。

次に、血行動態シミュレーションに用いるモデルについて説明する。
図４は、血行動態シミュレーションに用いるモデルの例を示す図である。
シミュレーション装置１００は、モデル１４０を記憶する。モデル１４０は、三次元モデル１４１および周囲モデル１４２を含む。三次元モデル１４１は、人体の心臓の左心室および右心室を表す心室モデルである。三次元モデル１４１は、１分間に６０回といった一定周期の拍動を表現するように作成される。周囲モデル１４２は、心臓の心室の外側であって、血液を循環させるための人体の主な部位を表す循環系モデルである。具体的には、周囲モデル１４２は、大動脈、大静脈、肺動脈、肺静脈、左心房および右心房を表す。心室から送出された血液が心室の外部を経由して再び心室に戻るという血液循環が表現されるように、三次元モデル１４１と周囲モデル１４２とが結合して使用される。

三次元モデル１４１は、有限要素法に用いられる有限要素モデルである。三次元モデル１４１では、左心室および右心室の三次元形状が、小さい四面体に分割される。四面体の頂点が節点であり、四面体の辺がエッジである。節点に変数が割り当てられる。有限要素モデルの構造例については後述する。周囲モデル１４２は、大動脈、大静脈、肺動脈、肺静脈、左心房および右心房の機能を電気回路の素子で表現した電気回路モデルである。周囲モデル１４２は、三次元モデル１４１より変数が少なく複雑度が低い。

周囲モデル１４２は、抵抗１４２－１～１４２－８およびキャパシタ１４２－９～１４２－１４を含む。抵抗１４２－１は、大動脈に対応し、抵抗１４２－２，１４２－７と隣接する。抵抗１４２－１は、抵抗値Ｒ_Ａをもつ。抵抗１４２－２は、大静脈に対応し、抵抗１４２－１，１４２－５と隣接する。抵抗１４２－２は、抵抗値Ｒ_Ｖをもつ。抵抗１４２－３は、肺動脈に対応し、抵抗１４２－４，１４２－６と隣接する。抵抗１４２－３は、抵抗値Ｒ_ＰＡをもつ。抵抗１４２－４は、肺静脈に対応し、抵抗１４２－３，１４２－８と隣接する。抵抗１４２－４は、抵抗値Ｒ_ＰＶをもつ。

抵抗１４２－５は、右心室の入口に対応し、抵抗値Ｒ_ＴＲをもつ。抵抗１４２－６は、右心室の出口に対応し、抵抗値Ｒ_ＰＵをもつ。抵抗１４２－７は、左心室の出口に対応し、抵抗値Ｒ_Ｃをもつ。抵抗１４２－８は、左心室の入口に対応し、抵抗値Ｒ_ＭＩをもつ。

キャパシタ１４２－９は、大動脈に対応し、抵抗１４２－７と抵抗１４２－１の間にある。キャパシタ１４２－９は、キャパシタンスＣ_Ａをもつ。キャパシタ１４２－９の電荷量Ｑ_Ａは大動脈の血液量に相当し、その圧力Ｐ_Ａに依存する。キャパシタ１４２－１０は、大静脈に対応し、抵抗１４２－１と抵抗１４２－２の間にある。キャパシタ１４２－１０は、キャパシタンスＣ_Ｖをもつ。キャパシタ１４２－１０の電荷量Ｑ_Ｖは大静脈の血液量に相当し、その圧力Ｐ_Ｖに依存する。

キャパシタ１４２－１１は、肺動脈に対応し、抵抗１４２－６と抵抗１４２－３の間にある。キャパシタ１４２－１１は、キャパシタンスＣ_ＰＡをもつ。キャパシタ１４２－１１の電荷量Ｑ_ＰＡは肺動脈の血液量に相当し、その圧力Ｐ_ＰＡに依存する。キャパシタ１４２－１２は、肺静脈に対応し、抵抗１４２－３と抵抗１４２－４の間にある。キャパシタ１４２－１２は、キャパシタンスＣ_ＰＶをもつ。キャパシタ１４２－１２の電荷量Ｑ_ＰＶは肺静脈の血液量に相当し、その圧力Ｐ_ＰＶに依存する。

キャパシタ１４２－１３は、右心房に対応し、抵抗１４２－２と抵抗１４２－５の間にある可変キャパシタ（可変コンデンサ）である。キャパシタ１４２－１３は、キャパシタンスの逆数に相当するエラスタンスＥ_ＲＡをもつ。エラスタンスＥ_ＲＡは、所定周期で規則的に変化する。キャパシタ１４２－１３の電荷量Ｑ_ＲＡは右心房の血液量に相当し、その圧力Ｐ_ＲＡに依存する。キャパシタ１４２－１４は、左心房に対応し、抵抗１４２－４と抵抗１４２－８の間にある可変キャパシタである。キャパシタ１４２－１４は、エラスタンスＥ_ＬＡをもつ。エラスタンスＥ_ＬＡは、所定周期で規則的に変化する。キャパシタ１４２－１４の電荷量Ｑ_ＬＡは左心房の血液量に相当し、その圧力Ｐ_ＬＡに依存する。

三次元モデル１４１の右心室の入力は圧力Ｐ_ＴＲをもち、右心室の出力は圧力Ｐ_ＲＶをもち、左心室の出力は圧力Ｐ_ＬＶをもち、左心室の入力は圧力Ｐ_ＭＩをもつ。左心室の出力は、抵抗１４２－７，１４２－１，１４２－２，１４２－５を通って、右心室の入力へと伝播する。また、右心室の出力は、抵抗１４２－６，１４２－３，１４２－４，１４２－８を通って、左心室の入力へと伝播する。このように、血液循環が電流で表現される。

抵抗値Ｒ_Ａ，Ｒ_Ｖ，Ｒ_ＰＡ，Ｒ_ＰＶ，Ｒ_ＴＲ，Ｒ_ＰＵ，Ｒ_Ｃ，Ｒ_ＭＩ、キャパシタンスＣ_Ａ，Ｃ_Ｖ，Ｃ_ＰＡ，Ｃ_ＰＶおよびエラスタンスＥ_ＲＡ，Ｅ_ＬＡは、シミュレーション前に予め値が与えられるパラメータである。圧力Ｐ_Ａ，Ｐ_Ｖ，Ｐ_ＰＡ，Ｐ_ＰＶ，Ｐ_ＲＡ，Ｐ_ＬＡは、シミュレーションの中で値が変化する変数である。電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡは、対応する圧力に依存して値が変化する変数である。電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡの初期値が決まれば、圧力Ｐ_Ａ，Ｐ_Ｖ，Ｐ_ＰＡ，Ｐ_ＰＶ，Ｐ_ＲＡ，Ｐ_ＬＡの初期値が決まる。

モデル１４０を用いた血行動態シミュレーションは、臨床医学に応用することが可能である。例えば、患者の現在の心臓を示す三次元モデル１４１を作成して、現在の血液循環を把握することが考えられる。また、手術などの治療行為を行った後の心臓を示す三次元モデル１４１を作成して、治療後の血液循環を予測することが考えられる。なお、周囲モデル１４２を患者に合わせるために、第２の実施の形態で説明するシミュレーションを、上記のパラメータの値を変えながら繰り返すようにしてもよい。

ここで、周囲モデル１４２に含まれるパラメータおよび変数の間の関係を説明する。シミュレーションでは、以下の関係式に従って連立方程式が生成される。
キャパシタ１４２－９および抵抗１４２－１に着目すると、この区間の体積流量Ｆ_ＡＯが数式（１）のように算出される。体積流量は、ある面を単位時間当たりに通過する流体の体積である。変数名の上に付したドットは、時間についての一階微分を表す。抵抗１４２－７に着目すると、体積流量Ｆ_ＡＯは数式（２）のようにも算出される。よって、数式（１）と数式（２）から、これらのパラメータおよび変数の間に成立する関係が規定される。また、抵抗１４２－１、キャパシタ１４２－１０および抵抗１４２－２に着目すると、この区間の体積流量の保存から数式（３）の関係が成立する。

また、抵抗１４２－２およびキャパシタ１４２－１３に着目すると、右心室の入口の体積流量Ｆ_ＴＲについて数式（４）の関係が成立する。抵抗１４２－５に着目すると、体積流量Ｆ_ＴＲは数式（５）のようにも算出される。よって、数式（４）と数式（５）から、これらのパラメータおよび変数の間に成立する関係が規定される。また、キャパシタ１４２－１１および抵抗１４２－３に着目すると、この区間の体積流量Ｆ_ＰＵが数式（６）のように算出される。抵抗１４２－６に着目すると、体積流量Ｆ_ＰＵは数式（７）のようにも算出される。よって、数式（６）と数式（７）から、これらのパラメータおよび変数の間に成立する関係が規定される。

また、抵抗１４２－３、キャパシタ１４２－１２および抵抗１４２－４に着目すると、この区間の体積流量の保存から数式（８）の関係が成立する。また、抵抗１４２－４およびキャパシタ１４２－１４に着目すると、左心室の入口の体積流量Ｆ_ＭＩについて数式（９）の関係が成立する。抵抗１４２－８に着目すると、体積流量Ｆ_ＭＩは数式（１０）のようにも算出される。よって、数式（９）と数式（１０）から、これらのパラメータおよび変数の間に成立する関係が規定される。

次に、有限要素モデルについて説明する。
図５は、有限要素モデルを説明するための図である。
三次元モデル１４１として、図５に示すような有限要素モデル１４４が使用される。なお、三次元モデル１４１の要素（メッシュ）は三次元の四面体であるものの、図５では説明を簡単にするため、要素を二次元の三角形として記載している。

有限要素モデル１４４は、節点１４４－１～１４４－９（節点ｎ１～ｎ９）を含む複数の節点（ノード）と、それら複数の節点の間を接続する複数のエッジ（リンク）とを含む。四面体の１つの要素は、４つの節点および６つのエッジによって形成され、三角形の面を４つもつ。有限要素モデル１４４の要素は、正四面体でなくてもよい。隣接する要素の間では、一部の節点および一部のエッジが共有される。

図５の例では、節点１４４－１，１４４－２，１４４－４によって１つの面が形成されている。節点１４４－２，１４４－４，１４４－５によって１つの面が形成されている。節点１４４－２，１４４－３，１４４－５によって１つの面が形成されている。節点１４４－３，１４４－５，１４４－６によって１つの面が形成されている。また、節点１４４－４，１４４－７，１４４－８によって１つの面が形成されている。節点１４４－４，１４４－５，１４４－８によって１つの面が形成されている。節点１４４－５，１４４－８，１４４－９によって１つの面が形成されている。節点１４４－５，１４４－６，１４４－９によって１つの面が形成されている。

複数の節点それぞれに、１以上の変数が割り当てられる。例えば、複数の節点それぞれに、その節点の位置における圧力を示す変数、その節点の当初の位置からの変位量を示す変数、その節点の位置における流速（流体速度）を示す変数などが割り当てられる。変位量の変数を割り当てるのは、心臓の拍動では心室が膨張と収縮を繰り返すためである。また、流速の変数は、血液が流れる領域の節点に割り当てられる。

例えば、節点１４４－１に、圧力を示す変数ｐ１、変位量を示す変数ｕ１、流速を示す変数ｆ１が割り当てられる。また、節点１４４－９に、圧力を示す変数ｐ９、変位量を示す変数ｕ９、流速を示す変数ｆ９が割り当てられる。

右心室の入口の圧力Ｐ_ＴＲは、右心室の入口の断面に相当する節点から、圧力を示す変数の値を抽出することで特定できる。同様に、右心室の出口の圧力Ｐ_ＲＶは、右心室の出口の断面に相当する節点から、圧力を示す変数の値を抽出することで特定できる。左心室の入口の圧力Ｐ_ＭＩは、左心室の入口の断面に相当する節点から、圧力を示す変数の値を抽出することで特定できる。左心室の出口の圧力Ｐ_ＬＶは、左心室の出口の断面に相当する節点から、圧力を示す変数の値を抽出することで特定できる。

シミュレーションにあたり、同一断面に属する複数の節点の間では圧力が均一であるという制約条件を設けるようにしてもよい。これらの圧力Ｐ_ＴＲ，Ｐ_ＲＶ，Ｐ_ＬＶ，Ｐ_ＭＩを通じて、三次元モデル１４１と周囲モデル１４２とが連携する。また、右心室の容積は、右心室を示す複数の要素の体積を合計することで算出できる。左心室の容積は、左心室を示す複数の要素の体積を合計することで算出できる。各節点の現在の位置座標は、当初の位置座標と変位量から算出できる。よって、ある要素の体積は、変位量を考慮して、その要素を形成する４つの節点それぞれの現在の位置座標から算出できる。

ここで、異なる変数の間に成立する関係を規定するための基礎方程式（支配方程式）について説明する。以下に説明する基礎方程式から連立方程式が生成される。前述の周囲モデル１４２の連立方程式と三次元モデル１４１の連立方程式とが結合されて、一組の連立方程式が形成され、係数行列が生成される。係数行列を用いた行列演算を繰り返す反復法によって、連立方程式の解が算出される。

三次元モデル１４１は、右心室の心筋部分の領域、右心室の流体部分の領域、左心室の心筋部分の領域および左心室の流体部分の領域を含む。心筋部分の物理的特性を示す基礎方程式は、数式（１１）に示す通りである。また、流体部分の物理的特性を示す基礎方程式は、数式（１２）に示す通りである。

数式（１１）において、ρは心筋の重量密度、ｐ_ｍは心筋圧、Ｊは体積の変形率である。また、ｕは変位量ベクトル、τ_ｆは心室内面Γ_ｆｓの上の表面力（Traction Force）ベクトルである。また、Ｅはグリーン・ラグランジェ（Green-Lagrange）ひずみテンソル、Ｓは第２ピオラ・キルヒホッフ（Piola-Kirchhoff）応力テンソル、Ｃは右コーシー・グリーン（Cauchy-Green）変形テンソルである。数式（１２）において、ρ_ｆは血液の重量密度、μ_ｆは血液の粘度、ｐは流体圧力である。また、ｖは速度ベクトル、ｃはメッシュに対する相対速度ベクトル、τ_ｓは心室内面Γ_ｆｓの上の表面力ベクトル、ｔ_ｃは血液の流入面Γ_ｃの上の表面力ベクトルである。また、Ｄは、変形速テンソルである。

次に、シミュレーション結果の収束について説明する。心臓の拍動を伴う血行動態シミュレーションでは、時刻ｔを微少量ずつ増加させていき、左心室および右心室の容積（心室容積）の時間変化と、左心室および右心室の圧力（心室圧）の時間変化を求める。心臓は同じ運動を一定周期で繰り返すため、安定状態における一周期分の心室容積および心室圧の時間変化を求めればよい。ただし、シミュレーションでは最初から安定状態の心室容積および心室圧を計算できるわけではなく、安定状態に到達するまでに時間を要する。そのため、収束するまで、複数周期分の心室容積および心室圧を計算することになる。

図６は、心室容積と心室圧のシミュレーション例を示すグラフである。
グラフ１５１は、心室容積と心室圧の組の時間変化を示すＰＶ（Pressure-Volume）ループを表したグラフである。曲線１５１－１は、開ループ（オープンループ）のモデルから算出された左心室の心室容積および心室圧を示す。曲線１５１－２は、開ループのモデルから算出された右心室の心室容積および心室圧を示す。曲線１５１－３は、閉ループ（クローズドループ）のモデルから算出された左心室の心室容積および心室圧を示す。曲線１５１－４は、閉ループのモデルから算出された右心室の心室容積および心室圧を示す。

開ループのモデルは、周辺モデルにおいて血流を途中でアースすることで、全血液量を保存しないモデルである。閉ループのモデルは、周辺モデルにおいて血流を途中でアースせず、全血液量を保存するモデルである。前述のモデル１４０は、閉ループのモデルである。グラフ１５１に示すように、開ループのモデルを用いた場合、ＰＶループの収束が比較的速い。一方、閉ループのモデルを用いた場合、ＰＶループの収束が遅い。ただし、閉ループのモデルの方が現実の血行動態との適合性が高く、シミュレーション結果の精度が高くなりやすい。このように、シミュレーションの前提条件によっては、収束が遅くなり、多くの周期数分の心室容積および心室圧を計算することになる場合がある。

収束が遅くなる原因の１つとして、変数の初期値が収束後の値（収束値）と乖離していることが挙げられる。第２の実施の形態の血行動態シミュレーションでは、シミュレーションの開始時点で、大動脈、大静脈、肺動脈、肺静脈、右心房および左心房に血液が存在する。よって、周囲モデル１４２のキャパシタ１４２－９～１４２－１４に対して、電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡの初期値を与える。この初期値が収束値と乖離していると、電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡが収束値に到達するまでに長時間を要し、結果として心室容積および心室圧が収束するまでに長時間を要することになる。

図７は、血液量のシミュレーションの収束例を示すグラフである。
グラフ１５２は、血液量の時間変化を表す時系列データである。曲線１５２－１は、大静脈の血液量を表しており、キャパシタ１４２－１０の電荷量Ｑ_Ｖに対応する。曲線１５２－２は、肺静脈の血液量を表しており、キャパシタ１４２－１２の電荷量Ｑ_ＰＶに対応する。曲線１５２－１，１５２－２に示すように、電荷量Ｑ_ｖ，Ｑ_ＰＶの初期値が１０秒後の収束値から乖離しているため、電荷量Ｑ_ｖ，Ｑ_ＰＶの収束が遅くなっている。

電荷量Ｑ_ｖ，Ｑ_ＰＶの初期値を収束値に近付けることができれば、電荷量Ｑ_ｖ，Ｑ_ＰＶの収束が速くなり、結果としてシミュレーションの所要周期数が少なくて済む。ただし、電荷量Ｑ_ｖ，Ｑ_ＰＶの収束値をユーザがシミュレーション前に見積もることは容易でない。そこで、シミュレーション装置１００は、以下のようにして適切な初期値を探索する。

まず、シミュレーション装置１００は、前述のモデル１４０を用いて、シミュレーションを所定周期数分だけ仮に実行する。モデル１４０を用いたシミュレーションでは、一周期分の心室容積および心室圧の計算に１時間程度要することがある。この点、仮のシミュレーションでは２～３周期分の計算を行えばよく、心室容積および心室圧が収束する前にシミュレーションを打ち切ってよい。三次元モデル１４１および周囲モデル１４２のパラメータには、予め用意したパラメータ値を代入する。周囲モデル１４２の変数（電荷量）には、仮の初期値を代入する。仮の初期値は、ランダムな数値でもよいし、ユーザが指定した数値でもよいし、所定の固定値でもよい。

図８は、心室容積の仮のシミュレーション例を示すグラフである。
グラフ１５３は、仮のシミュレーションによって算出される心室容積の時間変化を表す時系列データである。曲線１５３－１は、左心室の心室容積を表す。曲線１５３－２は、右心室の心室容積を表す。図８の例ではグラフ１５３に９拍分の心室容積が記載されているが、最初の所定拍数（例えば、３拍）分の心室容積を抽出できればよい。

図９は、心室圧の仮のシミュレーション例を示すグラフである。
グラフ１５４は、仮のシミュレーションによって算出される心室圧の時間変化を表す時系列データである。曲線１５４－１は、左心室の心室圧を表す。曲線１５４－２は、右心室の心室圧を表す。図９の例ではグラフ１５４に９拍分の心室圧が記載されているが、最初の所定拍数（例えば、３拍）分の心室圧を抽出できればよい。なお、仮のシミュレーションでは、前述のグラフ１５２が示す電荷量Ｑ_ｖ，Ｑ_ＰＶのように、所定拍数分の電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡも抽出しておく。

次に、シミュレーション装置１００は、所定周期数分の仮のシミュレーション結果を教師データとして用いて、三次元モデル１４１と置換される簡易モデルを生成する。簡易モデルは、仮のシミュレーション結果をできる限り高精度に再現可能なモデルであって、三次元モデル１４１よりも変数が少ないものである。第２の実施の形態では簡易モデルとして、有限要素モデルではなく電気回路モデルを使用する。

図１０は、三次元モデルの簡略化例を示す図である。
三次元モデル１４１を簡略化することでモデル１４０ａが得られる。モデル１４０ａは、三次元モデル１４１と置換された簡易モデル１４３と、周囲モデル１４２とを含む。周囲モデル１４２は、モデル１４０ａでもそのまま使用される。簡易モデル１４３は、有限要素モデルを近似する電気回路モデルであり、キャパシタ１４３－１，１４３－２を含む。キャパシタ１４３－１は、右心室に対応する可変キャパシタである。キャパシタ１４３－２は、左心室に対応する可変キャパシタである。

ここで、簡易モデル１４３に含まれるパラメータおよび変数の間に成立する関係を説明する。第２の実施の形態で使用する簡易モデル１４３は、心室の動作を簡易的に表現したモデルの一例である。他の簡易モデルが存在する場合、簡易モデル１４３に代えて他の簡易モデルを使用することも可能である。心室の動作を近似的に表現する方法は、医学論文などの専門的文献で提案される可能性もある。

モデル１４０ａを用いたシミュレーションでは、以下の関係式に従って連立方程式が生成される。前述の周囲モデル１４２の連立方程式と簡易モデル１４３の連立方程式とが結合されて、一組の連立方程式が形成され、係数行列が生成される。係数行列を用いた行列演算を繰り返す反復法によって、連立方程式の解が算出される。

時刻ｔにおける左心室の圧力Ｐ_ＬＶは、数式（１３）のように算出される。数式（１３）において、Ｅ_ＬＶはキャパシタ１４３－２のエラスタンスである。Ｖ_ＬＶは時刻ｔにおける左心室の容積、Ｖ_０，ＬＶは左心室の容積の最小値である。エラスタンスＥ_ＬＶは、数式（１４）に従って算出される。数式（１４）において、Ｅ_{ｍａｘ，ＬＶ}はキャパシタ１４３－２のエラスタンスの最大値である。Ｅ_Ｎは正規化された可変エラスタンスであり、後述するようにフーリエ係数として定義される。

時刻ｔにおける右心室の圧力Ｐ_ＲＶは、数式（１５）のように算出される。数式（１５）において、Ｅ_ＲＶはキャパシタ１４３－１のエラスタンスである。Ｖ_ＲＶは時刻ｔにおける右心室の容積、Ｖ_０，ＲＶは右心室の容積の最小値である。エラスタンスＥ_ＲＶは、数式（１６）に従って算出される。数式（１６）において、Ｅ_{ｍａｘ，ＲＶ}はキャパシタ１４３－１のエラスタンスの最大値である。Ｅ_Ｎは数式（１４）のものと同じである。正規化エラスタンスＥ_Ｎは、数式（１７）のように１２次のフーリエ係数として定義される。後述するように、振幅Ａ_０，Ａ_ｎおよび位相ρ_ｎがフーリエ係数である。

圧力Ｐ_ＬＶが左心室の心室圧であり、容積Ｖ_ＬＶが左心室の心室容積であり、圧力Ｐ_ＲＶが右心室の心室圧であり、容積Ｖ_ＲＶが右心室の心室容積である。圧力Ｐ_ＬＶ，Ｐ_ＲＶおよび容積Ｖ_ＬＶ，Ｖ_ＲＶは、シミュレーションの中で値が変化する変数である。三次元モデル１４１と異なり、簡易モデル１４３では心室容積および心室圧を直接表現する変数が存在する。最大エラスタンスＥ_{ｍａｘ，ＬＶ}，Ｅ_{ｍａｘ，ＲＶ}および最小心室容積Ｖ_０，ＬＶ，Ｖ_０，ＲＶは、シミュレーションの中で値が変化しないパラメータである。なお、左右心室の同期不全を表現するために、数式（１４）に代えて数式（１８）を用いてもよい。数式（１８）において、位相φは右心室に対する左心室の拍動の遅延時間を示すパラメータである。

図１１は、フーリエ係数テーブルの例を示す図である。
フーリエ係数テーブル１３５は、数式（１７）に出現するフーリエ係数の数値を示す。フーリエ係数テーブル１３５は、次数ｎ、振幅Ａ_ｎおよび位相ρ_ｎを対応付けている。次数ｎは、０次から１２次までの１３通りである。振幅Ａ_ｎとして、０次の振幅Ａ_０から１２次の振幅Ａ_１２までの１３個の振幅が登録されている。位相ρ_ｎとして、１次の位相ρ_１から１２次の位相ρ_１２までの１２個の位相が登録されている。

簡易モデル１４３の生成では、シミュレーション装置１００は、簡易モデル１４３のパラメータである最大エラスタンスＥ_{ｍａｘ，ＬＶ}，Ｅ_{ｍａｘ，ＲＶ}および最小心室容積Ｖ_０，ＬＶ，Ｖ_０，ＲＶに、仮のパラメータ値を代入する。仮のパラメータ値は、ランダムな値でもよいし、ユーザが指定した値でもよいし、所定の固定値でもよい。また、シミュレーション装置１００は、簡易モデル１４３の変数である圧力Ｐ_ＬＶ，Ｐ_ＲＶおよび容積Ｖ_ＬＶ，Ｖ_ＲＶに、モデル１４０を用いた前述のシミュレーションによって算出された心室圧と心室容積の初期値を代入する。この初期値は、例えば、グラフ１５３，１５４の先頭値である。

また、シミュレーション装置１００は、周囲モデル１４２のパラメータである抵抗値Ｒ_Ａ，Ｒ_Ｖ，Ｒ_ＰＡ，Ｒ_ＰＶ，Ｒ_ＴＲ，Ｒ_ＰＵ，Ｒ_Ｃ，Ｒ_ＭＩ、キャパシタンスＣ_Ａ，Ｃ_Ｖ，Ｃ_ＰＡ，Ｃ_ＰＶおよびエラスタンスＥ_ＲＡ，Ｅ_ＬＡに、前述のシミュレーションと同じ値を代入する。また、シミュレーション装置１００は、周囲モデル１４２の変数である電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡに、前述のシミュレーションと同じ初期値を代入する。

そして、シミュレーション装置１００は、簡易モデル１４３および周囲モデル１４２を含むモデル１４０ａを用いて、所定周期分のシミュレーションを実行する。ここで実行する周期数は、前述のシミュレーションと同じく２～３拍程度である。また、シミュレーション装置１００は、圧力Ｐ_ＬＶ，Ｐ_ＲＶ、容積Ｖ_ＬＶ，Ｖ_ＲＶおよび電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡのうち、１以上の変数を評価指標として選択する。シミュレーション装置１００は、選択した評価指標について、モデル１４０のシミュレーション結果とモデル１４０ａのシミュレーション結果との間の類似度を示す評価値を算出する。評価値は、例えば、２つの時系列データの間の誤差または相関係数である。

シミュレーション装置１００は、この評価値が改善するように簡易モデル１４３のパラメータ値を更新して、シミュレーションを繰り返す。評価値が誤差である場合は評価値が小さくなるようにパラメータ値を更新し、評価値が相関係数である場合は評価値が大きくなるようにパラメータ値を更新する。モデル１４０ａを用いた所定周期分のシミュレーションは、変数が少ないため１～２分程度で終わることがある。これにより、モデル１４０のシミュレーション結果を教師データとして簡易モデル１４３のパラメータが最適化され、モデル１４０を近似するようなモデル１４０ａを得ることができる。

簡易モデル１４３のパラメータが最適化されると、シミュレーション装置１００は、圧力Ｐ_ＬＶ，Ｐ_ＲＶおよび容積Ｖ_ＬＶ，Ｖ_ＲＶの周期的変化が収束するまで、モデル１４０ａのシミュレーションを実行する。収束条件は、最新の一周期分の時系列データと直前の一周期分の時系列データとが十分に類似していることである。類似度の評価値は、例えば、誤差または相関係数である。評価値が誤差である場合は評価値が閾値未満に低下したことが条件であり、評価値が相関係数である場合は評価値が閾値を超えたことが条件である。

このシミュレーションでは、簡易モデル１４３のパラメータ値としては、上記の最適化された値を用いる。簡易モデル１４３の変数の初期値、周囲モデル１４２のパラメータ値および周囲モデル１４２の変数の初期値としては、簡易モデル１４３を生成したときに用いたものと同じ値を用いる。シミュレーション装置１００は、モデル１４０ａのシミュレーションを収束するまで実行すると、収束条件を満たした周期の先頭の電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡを周囲モデル１４２から抽出する。これら電荷量が、周囲モデル１４２の変数に対する好ましい初期値と推定される。

周囲モデル１４２の好ましい初期値が決定されると、シミュレーション装置１００は、モデル１４０に戻って本来のシミュレーションを実行する。この本来のシミュレーションは、収束条件が満たされるまで継続して実行される。収束条件は、モデル１４０ａのシミュレーションと同じく、最新の一周期分の心室圧および心室容積の時系列データと直前の一周期分の時系列データとが十分に類似していることである。本来のシミュレーションでは、三次元モデル１４１のパラメータ値および周囲モデル１４２のパラメータ値として、前述の所定周期分のシミュレーションで用いたものと同じ値を用いる。また、周囲モデル１４２の変数である電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡの初期値として、モデル１４０ａを用いて算出された好ましい初期値を用いる。

これにより、任意の初期値を周囲モデル１４２の変数に代入して、モデル１４０のシミュレーションを収束するまで実行する場合と比べて、収束に要する周期数（拍数）が少なくなる。よって、シミュレーションの計算量を削減でき所要時間を短縮できる。なお、簡易モデル１４３は三次元モデル１４１よりも著しく変数が少ないため、モデル１４０ａを用いたシミュレーションの実行時間は短い。よって、モデル１４０ａのシミュレーションを行ってもトータルの所要時間は短縮される。

次に、シミュレーション装置１００の機能について説明する。
図１２は、シミュレーション装置の機能例を示すブロック図である。
シミュレーション装置１００は、メッシュデータ記憶部１２１およびパラメータ記憶部１２２を有する。これらの記憶部は、例えば、ＲＡＭ１０２またはＨＤＤ１０３の記憶領域を用いて実現される。また、シミュレーション装置１００は、三次元シミュレーション部１２３、簡易シミュレーション部１２４、指標算出部１２５、パラメータ推定部１２６、初期値決定部１２７および可視化部１２８を有する。これらの処理部は、例えば、ＣＰＵ１０１が実行するプログラムを用いて実現される。

メッシュデータ記憶部１２１は、有限要素モデルである三次元モデル１４１に相当する三次元メッシュデータを記憶する。三次元メッシュデータは、心室の形状を四面体の要素（メッシュ）に細分化したものである。例えば、ＣＡＤアプリケーションによって、心室の形状を示すＣＡＤ（Computer Aided Design）データが生成され、メッシュ生成アプリケーションによって、ＣＡＤデータから三次元メッシュデータが生成される。ＣＡＤアプリケーションやメッシュ生成アプリケーションは、シミュレーション装置１００で実行されてもよいし、他の情報処理装置で実行されてもよい。

パラメータ記憶部１２２は、三次元モデル１４１に含まれる心筋の重力密度ρや血液の重力密度ρ_ｆなどのパラメータの値を記憶する。また、パラメータ記憶部１２２は、周囲モデル１４２に含まれるパラメータである抵抗値Ｒ_Ａ，Ｒ_Ｖ，Ｒ_ＰＡ，Ｒ_ＰＶ，Ｒ_ＴＲ，Ｒ_ＰＵ，Ｒ_Ｃ，Ｒ_ＭＩ、キャパシタンスＣ_Ａ，Ｃ_Ｖ，Ｃ_ＰＡ，Ｃ_ＰＶおよびエラスタンスＥ_ＲＡ，Ｅ_ＬＡの値を記憶する。これらのパラメータ値は、例えば、ユーザにより設定される。ただし、前述のように、シミュレーション装置１００を用いてパラメータ値の探索を行ってもよい。

三次元シミュレーション部１２３は、三次元モデル１４１および周囲モデル１４２を用いたシミュレーションを行う。三次元シミュレーション部１２３は、三次元モデル１４１に含まれる変数およびパラメータの間の関係を示す連立方程式を、前述の基礎方程式に基づいて生成する。また、三次元シミュレーション部１２３は、周囲モデル１４２に含まれる変数およびパラメータの間の関係を示す連立方程式を、前述の関係式に基づいて生成する。三次元シミュレーション部１２３は、これらの連立方程式の係数行列を生成し、反復法による行列演算によって連立方程式の解を求める。三次元シミュレーション部１２３は、複数のＣＰＵや複数のＧＰＵを用いて行列演算を並列処理してもよい。

簡易シミュレーション部１２４は、簡易モデル１４３および周囲モデル１４２を用いたシミュレーションを行う。簡易シミュレーション部１２４は、簡易モデル１４３に含まれる変数およびパラメータの間の関係を示す連立方程式を、前述の関係式に基づいて生成する。また、簡易シミュレーション部１２４は、周囲モデル１４２に含まれる変数およびパラメータの間の関係を示す連立方程式を、前述の関係式に基づいて生成する。簡易シミュレーション部１２４は、これらの連立方程式の係数行列を生成し、反復法による行列演算によって連立方程式の解を求める。簡易シミュレーション部１２４は、複数のＣＰＵや複数のＧＰＵを用いて行列演算を並列処理してもよい。

指標算出部１２５は、三次元シミュレーション部１２３に、所定周期分の三次元シミュレーションを実行させる。指標算出部１２５は、三次元シミュレーション部１２３の実行結果から、左心室の心室圧、右心室の心室圧、左心室の心室容積および右心室の心室容積の所定周期分の時系列データを生成する。また、指標算出部１２５は、電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡの所定周期分の時系列データを抽出する。

パラメータ推定部１２６は、指標算出部１２５が生成した１０個の指標の時系列データのうちの一部または全部を、教師データとして選択する。パラメータ推定部１２６は、選択した教師データに基づいて、簡易モデル１４３のパラメータ値を修正しながら、簡易シミュレーション部１２４に所定周期分の簡易シミュレーションを繰り返し実行させる。これにより、パラメータ推定部１２６は、簡易モデル１４３のパラメータ値を決定する。

初期値決定部１２７は、パラメータ推定部１２６が決定したパラメータ値を指定して、簡易シミュレーション部１２４に簡易シミュレーションを収束するまで実行させる。初期値決定部１２７は、簡易シミュレーション部１２４の実行結果から、最後の周期の先頭時刻における電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡを抽出する。

可視化部１２８は、初期値決定部１２７が抽出した電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡを周囲モデル１４２の変数の初期値に指定して、三次元シミュレーション部１２３に三次元シミュレーションを収束するまで実行させる。可視化部１２８は、三次元シミュレーション部１２３の実行結果から、左心室の心室圧、右心室の心室圧、左心室の心室容積および右心室の心室容積の時系列データを生成する。可視化部１２８は、生成した時系列データを出力する。出力する時系列データは、収束した一周期分のものでもよい。例えば、可視化部１２８は、時系列データをグラフとして可視化して表示装置１１１に表示させる。ただし、可視化部１２８は、時系列データを不揮発性ストレージに保存してもよく、ネットワーク１１４を介して他の情報処理装置に転送してもよい。

図１３は、節点テーブルとメッシュテーブルの例を示す図である。
節点テーブル１３１およびメッシュテーブル１３２は、心臓の心室の形状を示すメッシュデータであり、メッシュデータ記憶部１２１に記憶される。節点テーブル１３１は、節点番号と座標とを対応付けた複数のレコードを記憶する。節点番号は、節点を識別する識別子である。座標は、節点の位置を示すＸ座標、Ｙ座標およびＺ座標の組である。メッシュテーブル１３２は、メッシュ番号と節点番号とを対応付けた複数のレコードを記憶する。メッシュ番号は、四面体の要素（メッシュ）を識別する識別子である。節点番号は、四面体の４つの頂点に相当する４つの節点の節点番号を列挙したものである。

図１４は、パラメータテーブルの例を示す図である。
パラメータテーブル１３３は、パラメータ記憶部１２２に記憶される。パラメータテーブル１３３は、パラメータとその値（パラメータ値）とを対応付けた複数のレコードを含む。パラメータには、周囲モデル１４２の抵抗値Ｒ_Ａ，Ｒ_Ｖ，Ｒ_ＰＡ，Ｒ_ＰＶ，Ｒ_ＴＲ，Ｒ_ＰＵ，Ｒ_Ｃ，Ｒ_ＭＩ、キャパシタンスＣ_Ａ，Ｃ_Ｖ，Ｃ_ＰＡ，Ｃ_ＰＶおよびエラスタンスＥ_ＲＡ，Ｅ_ＬＡが含まれる。また、パラメータには、三次元モデル１４１の重力密度ρ，ρ_ｆなど、前述の基礎方程式に出現する各種のパラメータが含まれる。

図１５は、時系列データテーブルの例を示す図である。
時系列データテーブル１３４は、指標算出部１２５によって生成される。時系列データテーブル１３４は、１０個の指標の所定周期分（例えば、３拍分）の時系列データを記憶する。１０個の指標は、左心室容積Ｖ_ＬＶ、右心室容積Ｖ_ＲＶ、左心室圧Ｐ_ＬＶ、右心室圧Ｐ_ＲＶおよび電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡである。時系列データテーブル１３４は、所定の時間刻みで、これら１０個の指標の数値を列挙する。

次に、シミュレーション装置１００の処理手順について説明する。
図１６は、シミュレーションの手順例を示すフローチャートである。
（Ｓ１０）シミュレーション装置１００は、メッシュデータ、三次元モデル１４１のパラメータ値および周囲モデル１４２のパラメータ値を取得する。

（Ｓ１１）指標算出部１２５は、周囲モデル１４２の電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡの初期値を設定する。初期値は、ランダムでもよく所定の固定値でもよい。

（Ｓ１２）指標算出部１２５は、三次元モデル１４１のパラメータ値、周囲モデル１４２のパラメータ値、ステップＳ１１の周囲モデル１４２の初期値、および、拍数を指定する。拍数は、例えば、２～３拍である。三次元シミュレーション部１２３は、指標算出部１２５からの指定に従って三次元シミュレーションを実行する。

（Ｓ１３）指標算出部１２５は、三次元シミュレーション部１２３の実行結果から、左心室容積Ｖ_ＬＶ、右心室容積Ｖ_ＲＶ、左心室圧Ｐ_ＬＶ、右心室圧Ｐ_ＲＶおよび電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡの時系列データを生成する。各時刻の左心室容積Ｖ_ＬＶおよび右心室容積Ｖ_ＲＶは、節点に対応付けられた変数である変位量の値から各四面体の体積を求め、これを合算することで算出することができる。各時刻の左心室圧Ｐ_ＬＶおよび右心室圧Ｐ_ＲＶは、心室出口の節点に対応付けられた変数である圧力の値に相当する。各時刻の電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡは、圧力Ｐ_Ａ，Ｐ_Ｖ，Ｐ_ＰＡ，Ｐ_ＰＶ，Ｐ_ＲＡ，Ｐ_ＬＡとキャパシタンスの積として算出される。

（Ｓ１４）パラメータ推定部１２６は、左心室容積Ｖ_ＬＶ、右心室容積Ｖ_ＲＶ、左心室圧Ｐ_ＬＶ、右心室圧Ｐ_ＲＶおよび電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡの１０個の指標から、１以上の評価指標を選択する。好ましくは、１０個の指標を全て選択する。

（Ｓ１５）パラメータ推定部１２６は、ステップＳ１３の左心室容積Ｖ_ＬＶ、右心室容積Ｖ_ＲＶ、左心室圧Ｐ_ＬＶおよび右心室圧Ｐ_ＲＶの時系列データから、それぞれ先頭の値を抽出し、抽出した値を簡易モデル１４３の初期値に設定する。

（Ｓ１６）パラメータ推定部１２６は、簡易モデル１４３の最大エラスタンスＥ_{ｍａｘ，ＬＶ}，Ｅ_{ｍａｘ，ＲＶ}および最小心室容積Ｖ_０，ＬＶ，Ｖ_０，ＲＶにパラメータ値を設定する。パラメータ値は、ランダムでもよく所定の固定値でもよい。

（Ｓ１７）パラメータ推定部１２６は、ステップＳ１２と同じ周囲モデル１４２のパラメータ値、周囲モデル１４２の初期値、および、拍数を指定する。また、パラメータ推定部１２６は、ステップＳ１５の簡易モデル１４３の初期値、および、ステップＳ１６の簡易モデル１４３のパラメータ値を指定する。簡易シミュレーション部１２４は、パラメータ推定部１２６からの指定に従って簡易シミュレーションを実行する。

（Ｓ１８）パラメータ推定部１２６は、簡易シミュレーション部１２４の実行結果から、左心室容積Ｖ_ＬＶ、右心室容積Ｖ_ＲＶ、左心室圧Ｐ_ＬＶ、右心室圧Ｐ_ＲＶおよび電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡの時系列データを生成する。パラメータ推定部１２６は、ステップＳ１４で選択した評価指標について、ここで生成した時系列データとステップＳ１３の時系列データとの間の類似度を示す評価値を算出する。２以上の評価指標が選択されている場合、例えば、評価指標毎の評価値の平均を全体の評価値とする。評価値は、例えば、誤差（平均二乗誤差や二乗平均平方根誤差など）または相関係数である。

（Ｓ１９）パラメータ推定部１２６は、ステップＳ１８の評価値が所定条件を満たしているか判断する。所定条件は、例えば、誤差を示す評価値が所定の閾値未満であること、または、相関係数を示す評価値が所定の閾値を超えていることである。評価値が所定条件を満たす場合はステップＳ２１に進み、それ以外の場合はステップＳ２０に進む。

（Ｓ２０）パラメータ推定部１２６は、評価値が改善するように（例えば、誤差が小さくなる、または、相関係数が大きくなるように）、簡易モデル１４３の最大エラスタンスＥ_{ｍａｘ，ＬＶ}，Ｅ_{ｍａｘ，ＲＶ}および最小心室容積Ｖ_０，ＬＶ，Ｖ_０，ＲＶのパラメータ値を修正する。そして、ステップＳ１７に戻る。次回のステップＳ１７では、簡易モデル１４３のパラメータ値以外のパラメータ値、初期値および拍数は、前回と同じものを使用する。

図１７は、シミュレーションの手順例を示すフローチャート（続き）である。
（Ｓ２１）パラメータ推定部１２６は、所定条件を満たす評価値が得られた時点（例えば、誤差が閾値未満になった時点、または、相関係数が閾値を超えた時点）における簡易モデル１４３のパラメータ値を、最適なパラメータ値と決定する。初期値決定部１２７は、ステップＳ１７と同じ周囲モデル１４２のパラメータ値、周囲モデル１４２の初期値、および、簡易モデル１４３の初期値を指定する。また、初期値決定部１２７は、パラメータ推定部１２６が決定した簡易モデル１４３のパラメータ値を指定する。

（Ｓ２２）簡易シミュレーション部１２４は、初期値決定部１２７からの指定に従って、簡易シミュレーションを一周期分だけ実行する。
（Ｓ２３）初期値決定部１２７は、左心室容積Ｖ_ＬＶ、右心室容積Ｖ_ＲＶ、左心室圧Ｐ_ＬＶおよび右心室圧Ｐ_ＲＶの最新の一周期分の時系列データを抽出する。初期値決定部１２７は、ここで抽出した時系列データと直前の一周期分の時系列データとの間の類似度を示す評価値を算出する。例えば、上記の４つの評価指標の評価値の平均を全体の評価値とする。評価値は、例えば、誤差または相関係数である。なお、直前の一周期分の時系列データがまだ無い場合は、ステップＳ２４の判断がＮＯになる。

（Ｓ２４）初期値決定部１２７は、ステップＳ２３の評価値が収束を示しているか判断する。例えば、誤差を示す評価値が所定の閾値未満であるとき、または、相関係数を示す評価値が所定の閾値を超えたとき、収束していると判断する。評価値が収束を示している場合はステップＳ２５に進み、それ以外の場合はステップＳ２２に戻る。

（Ｓ２５）初期値決定部１２７は、最終周期の先頭時刻における周囲モデル１４２の電荷量Ｑ_Ａ，Ｑ_Ｖ，Ｑ_ＰＡ，Ｑ_ＰＶ，Ｑ_ＲＡ，Ｑ_ＬＡを抽出する。
（Ｓ２６）可視化部１２８は、ステップＳ１２の三次元モデル１４１のパラメータ値、周囲モデル１４２のパラメータ値、および、周囲モデル１４２の初期値のうち、周囲モデル１４２の初期値をステップＳ２５で抽出された値に変更する。

（Ｓ２７）三次元シミュレーション部１２３は、可視化部１２８からの指定に従って、三次元シミュレーションを一周期分だけ実行する。
（Ｓ２８）可視化部１２８は、三次元シミュレーション部１２３の実行結果から、左心室容積Ｖ_ＬＶ、右心室容積Ｖ_ＲＶ、左心室圧Ｐ_ＬＶおよび右心室圧Ｐ_ＲＶの最新の一周期分の時系列データを生成する。各時刻の左心室容積Ｖ_ＬＶおよび右心室容積Ｖ_ＲＶは、各四面体の体積を合算することで算出することができる。各時刻の左心室圧Ｐ_ＬＶおよび右心室圧Ｐ_ＲＶは、心室出口の節点に対応付けられた変数である圧力の値に相当する。

（Ｓ２９）可視化部１２８は、ステップＳ２８の時系列データと直前の一周期分の時系列データとの間の類似度を示す評価値を算出する。例えば、上記の４つの評価指標の評価値の平均を全体の評価値とする。評価値は、例えば、誤差または相関係数である。なお、直前の一周期分の時系列データがまだ無い場合は、ステップＳ３０の判断がＮＯになる。

（Ｓ３０）可視化部１２８は、ステップＳ２９の評価値が収束を示しているか判断する。例えば、誤差を示す評価値が所定の閾値未満であるとき、または、相関係数を示す評価値が所定の閾値を超えたとき、収束していると判断する。評価値が収束を示している場合はステップＳ３１に進み、それ以外の場合はステップＳ２７に戻る。

（Ｓ３１）可視化部１２８は、左心室容積Ｖ_ＬＶ、右心室容積Ｖ_ＲＶ、左心室圧Ｐ_ＬＶおよび右心室圧Ｐ_ＲＶの時間変化を示すグラフを生成し、表示装置１１１に表示させる。
第２の実施の形態のシミュレーション装置１００によれば、三次元モデル１４１と接続される周囲モデル１４２の変数に、収束後の値に近似する適切な初期値をシミュレーション開始時点で代入することが可能となる。よって、三次元モデル１４１および周囲モデル１４２を用いた血行動態シミュレーションが収束するまでの所要周期数を少なく抑えることができる。そのため、血行動態シミュレーションの計算量を削減でき、所要時間を短縮することができる。また、適切な初期値を求めるために使用する簡易モデル１４３は、三次元モデル１４１よりも変数が大幅に少ない。このため、簡易モデル１４３および周囲モデル１４２を用いたシミュレーションは短時間で完了することが可能である。そのため、簡易モデル１４３および周囲モデル１４２を用いたシミュレーションを追加しても、トータルの計算量は削減され、トータルの所要時間を短縮できる。

１０ シミュレーションシステム
１１ 記憶装置
１２ 演算装置
１３，１４，１５ モデル
１６ａ，１６ｂ 状態値
１７ 特徴量

Claims (11)

1. コンピュータに、
   周期的な運動を行う第１の部位を示す第１のモデルと、前記第１の部位に接続されており前記周期的な運動の影響を受ける第２の部位を示す第２のモデルとを取得し、
   前記第２のモデルに含まれる変数の初期値として第１の状態値を割り当て、前記第１のモデルおよび前記第２のモデルと前記第１の状態値とを用いて、前記周期的な運動のシミュレーションを所定周期数分実行して、前記所定周期数分の特徴量を算出し、
   前記所定周期数分の特徴量に基づいて、前記第１のモデルより変数が少なく前記第１のモデルと置換可能な第３のモデルを生成し、
   前記第３のモデルおよび前記第２のモデルを用いて、前記周期的な運動のシミュレーションを所定の収束条件が満たされるまで継続して実行し、前記所定の収束条件が満たされた周期において前記第２のモデルに含まれる変数がもつ第２の状態値を抽出する、
   処理を実行させるシミュレーションプログラム。
2. 前記コンピュータに更に、
   前記第２のモデルに含まれる変数の初期値として前記第２の状態値を割り当て、前記第１のモデルおよび前記第２のモデルと前記第２の状態値とを用いて、前記周期的な運動のシミュレーションを実行する、処理を実行させる
   請求項１記載のシミュレーションプログラム。
3. 前記所定の収束条件は、最新の周期の特徴量と前記最新の周期から１つ前の周期の特徴量との差が、閾値未満であることである、
   請求項１記載のシミュレーションプログラム。
4. 前記第３のモデルの生成では、前記所定周期数分の特徴量と、前記第３のモデルおよび前記第２のモデルを用いて前記周期的な運動のシミュレーションを前記所定周期数分実行した結果との間の差に基づいて、前記第３のモデルのパラメータの値を修正する、
   請求項１記載のシミュレーションプログラム。
5. 前記第２のモデルは、キャパシタを含む電気回路を示す電気回路モデルであり、
   前記第１の状態値および前記第２の状態値は、前記キャパシタの電荷量に対応する、
   請求項１記載のシミュレーションプログラム。
6. 前記第１のモデルは、複数の節点と前記複数の節点の間を接続する複数のエッジとを含み、前記複数の節点それぞれに変数が割り当てられる有限要素モデルであり、
   前記第２のモデルおよび前記第３のモデルは、電気回路を示す電気回路モデルである、
   請求項１記載のシミュレーションプログラム。
7. 前記周期的な運動のシミュレーションは、前記第１の部位と前記第２の部位との間で流体を循環させる運動のシミュレーションである、
   請求項１記載のシミュレーションプログラム。
8. 前記所定周期数分の特徴量は、前記第１の部位の容積の時間変化を示す第１の時系列データと、前記第１の部位の圧力の時間変化を示す第２の時系列データとを含む、
   請求項１記載のシミュレーションプログラム。
9. 前記第１の部位は、心臓の少なくとも一部分であり、
   前記第２の部位は、前記心臓の外部の血管を含む、
   請求項１記載のシミュレーションプログラム。
10. コンピュータが、
    周期的な運動を行う第１の部位を示す第１のモデルと、前記第１の部位に接続されており前記周期的な運動の影響を受ける第２の部位を示す第２のモデルとを取得し、
    前記第２のモデルに含まれる変数の初期値として第１の状態値を割り当て、前記第１のモデルおよび前記第２のモデルと前記第１の状態値とを用いて、前記周期的な運動のシミュレーションを所定周期数分実行して、前記所定周期数分の特徴量を算出し、
    前記所定周期数分の特徴量に基づいて、前記第１のモデルより変数が少なく前記第１のモデルと置換可能な第３のモデルを生成し、
    前記第３のモデルおよび前記第２のモデルを用いて、前記周期的な運動のシミュレーションを所定の収束条件が満たされるまで継続して実行し、前記所定の収束条件が満たされた周期において前記第２のモデルに含まれる変数がもつ第２の状態値を抽出する、
    シミュレーション方法。
11. 周期的な運動を行う第１の部位を示す第１のモデルと、前記第１の部位に接続されており前記周期的な運動の影響を受ける第２の部位を示す第２のモデルとを記憶する記憶装置と、
    前記第２のモデルに含まれる変数の初期値として第１の状態値を割り当て、前記第１のモデルおよび前記第２のモデルと前記第１の状態値とを用いて、前記周期的な運動のシミュレーションを所定周期数分実行して、前記所定周期数分の特徴量を算出し、前記所定周期数分の特徴量に基づいて、前記第１のモデルより変数が少なく前記第１のモデルと置換可能な第３のモデルを生成し、前記第３のモデルおよび前記第２のモデルを用いて、前記周期的な運動のシミュレーションを所定の収束条件が満たされるまで継続して実行し、前記所定の収束条件が満たされた周期において前記第２のモデルに含まれる変数がもつ第２の状態値を抽出する演算装置と、
    を有するシミュレーションシステム。

Images