JP6091913B2 - 生体シミュレーションプログラム、生体シミュレーション方法及び生体シミュレーション装置 - Google Patents

Description

本発明は、生体シミュレーションプログラム、生体シミュレーション方法及び生体シミュレーション装置に関する。

医用計測技術の発展に伴い、生体の臓器、細胞、細胞に含まれる物質のそれぞれのスケールに関しては、膨大な知見が得られている。ところが、臓器、細胞、細胞に含まれる物質のスケール間の関連、例えば、細胞と臓器との挙動が互いにどのように影響しているのかについては、専門家にとっても理解し予測することが困難である。例えば、心臓においては、家族性肥大型心筋症など、網羅的遺伝子解析によって原因遺伝子が同定されているにも関わらず、病態形成のメカニズム解明に至っていない心疾患がある。

そこで、細胞と臓器との挙動が互いにどのように影響しているのかをシミュレーションする技術がある。かかる技術の一例としては、心筋細胞のモデルと、心臓のモデルとの挙動をシミュレーションする技術が挙げられる。

A. Hosoi, T. Washio, J. Okada, Y. Kadooka, K. Nakajima and T. Hisada, "A multi-scale heart simulation on massively parallel computers," ACM/IEEE Supercomputing Conference (SC10), 2010.

しかしながら、上述した心筋細胞のモデルと心臓のモデルとの挙動をシミュレーションする技術は、心筋細胞に含まれる構造、例えば、サルコメアについては、平均モデルを用いてシミュレーションを行う。このため、かかる技術は、サルコメアなどの心筋細胞に含まれる小さい構造については、平均モデルを用いてシミュレーションを行うので、シミュレーションの結果が精緻さに欠けるという問題がある。

１つの側面では、精緻なシミュレーションの結果を得ることを目的とする。

本願の開示する生体シミュレーションプログラムは、コンピュータに、生体の臓器の細胞に含まれる物質の挙動を計算させる。生体シミュレーションプログラムは、コンピュータに、計算された物質の挙動に基づいて、細胞の挙動を計算させる。生体シミュレーションプログラムは、コンピュータに、計算された細胞の挙動に基づいて、臓器の挙動を計算させる。生体シミュレーションプログラムは、コンピュータに、計算された臓器の挙動を、細胞の挙動に反映させる。生体シミュレーションプログラムは、コンピュータに、臓器の挙動が反映された細胞の挙動を、物質の挙動に反映させる。

精緻なシミュレーションの結果を得ることができる。

図１は、実施例１に係る生体シミュレーション装置の機能的な構成の一例を示す図である。 図２は、マクロモデル、メゾモデル及びミクロモデルの一例を示す図である。 図３は、条件データの一例を示す図である。 図４は、実施例１に係る生体シミュレーション処理の手順を示すフローチャートである。 図５は、実施例１に係るシミュレーション処理の手順を示すフローチャートである。 図６Ａは、表示部に表示された、あるタイムステップにおける各モデルのシミュレーション結果の一例である。 図６Ｂは、表示部に表示された、あるタイムステップにおける各モデルのシミュレーション結果の一例である。 図６Ｃは、表示部に表示された、あるタイムステップにおける各モデルのシミュレーション結果の一例である。 図７は、実施例１に係るマクロ処理部及びメゾ・ミクロ処理部の処理の実行タイミングの一例を示すタイムチャートである。 図８は、実施例２に係る生体シミュレーション装置の機能的な構成の一例を示す図である。 図９は、実施例２に係るシミュレーション処理の手順を示すフローチャートである。 図１０は、実施例２に係るマクロ処理部及びメゾ・ミクロ処理部の処理の実行タイミングの一例を示すタイムチャートである。 図１１は、生体シミュレーションプログラムを実行するコンピュータを示す図である。

以下に、本願の開示する生体シミュレーションプログラム、生体シミュレーション方法及び生体シミュレーション装置の各実施例を図面に基づいて詳細に説明する。なお、各実施例は開示の技術を限定するものではない。

実施例１に係る生体シミュレーション装置について説明する。図１は、実施例１に係る生体シミュレーション装置の機能的な構成の一例を示す図である。生体シミュレーション装置１０は、入力部１１、表示部１２、記憶部１３及び制御部１４を有する。

入力部１１は、制御部１４に各種の情報を入力する。例えば、入力部１１は、ユーザから後述の生体シミュレーション処理を実行する指示を受け付けた場合には、受け付けた指示を制御部１４に入力する。入力部１１のデバイスの一例としては、キーボードやマウスなどが挙げられる。

表示部１２は、各種の情報を表示する。例えば、表示部１２は、後述の表示制御部１４ｃの制御によりシミュレーション結果を表示する。

記憶部１３は、制御部１４で実行される各種プログラムを記憶する。また、記憶部１３は、マクロモデル１３ａ、メゾモデル１３ｂ、ミクロモデル１３ｃ及び条件データ１３ｄを記憶する。

マクロモデル１３ａは、心臓全体のモデルを複数個（例えば１０個）に分割したモデルである。なお、心臓は、生体の臓器の一例である。また、メゾモデル１３ｂは、心筋細胞のモデルである。なお、心筋細胞は、細胞の一例である。また、ミクロモデル１３ｃは、心筋細胞に含まれるサルコメアのモデルである。なお、サルコメアは、細胞に含まれる物質の一例である。図２は、マクロモデル、メゾモデル及びミクロモデルの一例を示す図である。図２の例は、心臓全体のモデル２０を複数のマクロモデル１３ａに分割した場合を示す。複数のマクロモデル１３ａのそれぞれが示す心臓の部分は、複数の心筋細胞２１を含む。メゾモデル１３ｂは、マクロモデル１３ａに含まれる各心筋細胞２１のモデルである。図２の例では、メゾモデル１３ｂは、形状が正四面体のモデルであり、４つの節点２２を有する。メゾモデル１３ｂが示す心筋細胞は、複数のサルコメア２３を含む。ミクロモデル１３ｃは、メゾモデル１３ｂに含まれる各サルコメア２３のモデルである。図２の例では、ミクロモデル１３ｃは、形状が六面体のモデルであり、８つの節点２４を有する。

条件データ１３ｄは、生体シミュレーションを行う場合の条件を示すデータである。図３は、条件データの一例を示す図である。図３の例では、条件データ１３ｄは、「タイムステップ」及び「Ｃａ^２＋濃度」の各項目を有する。「タイムステップ」の項目には、生体シミュレーションのタイムステップが登録される。「Ｃａ^２＋濃度」の項目には、「タイムステップ」の項目に登録されたタイムステップでの各心筋細胞におけるカルシウムイオン（Ｃａ^２＋）濃度が登録される。すなわち、図３の例に示す条件データ１３ｄは、各タイムステップの各心筋細胞におけるカルシウムイオン濃度を示す。

記憶部１３は、例えば、フラッシュメモリなどの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。なお、記憶部１３は、上記の種類の記憶装置に限定されるものではなく、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）等であってもよい。

制御部１４は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。図１に示すように、複数の第２の演算処理部であるマクロ処理部１４ａと、複数の第１の演算処理部であるメゾ・ミクロ処理部１４ｂと、表示制御部１４ｃとを有する。

複数のマクロ処理部１４ａのそれぞれは、複数のマクロモデル１３ａのそれぞれに対応し、各マクロ処理部１４ａは、対応するマクロモデル１３ａの挙動を計算する。また、複数のメゾ・ミクロ処理部１４ｂのそれぞれは、複数のメゾモデル１３ｂのそれぞれに対応し、各メゾ・ミクロ処理部１４ｂは、対応するメゾモデル１３ｂの挙動を計算する。さらに、各メゾ・ミクロ処理部１４ｂは、対応するメゾモデルが示す心筋細胞２１に含まれる各サルコメア２３を示す各ミクロモデル１３ｃの挙動を計算する。すなわち、本実施例では、１つのマクロモデル１３ａに対して、複数のメゾモデル１３ｂが対応するため、１つのマクロ処理部１４ａに対して複数のメゾ・ミクロ処理部１４ｂが対応することとなる。また、制御部１４は、演算処理部であるコアを複数有するマルチコアのＣＰＵ（Central Processing Unit）を１つ又は複数有する。若しくは、演算処理部であるコアを１つ有するシングルコアのＣＰＵを複数有してもよい。マクロ処理部１４ａ、メゾ・ミクロ処理部１４ｂ及び表示制御部１４ｃの各部は、後述の生体シミュレーションプログラムを各コアが実行することにより、各コア上で実現される。

第２の演算処理部であるマクロ処理部１４ａは、計算部１５ａを有する。第１の演算処理部であるメゾ・ミクロ処理部１４ｂは、第１の計算部１５ｂ、第２の計算部１５ｃ、第１の更新部１５ｄ、第２の更新部１５ｅを有する。

ここで、制御部１４で実行される処理の一例について図４を参照して説明する。図４は、実施例１に係る生体シミュレーション処理の手順を示すフローチャートである。生体シミュレーション処理は、例えば、入力部１１から制御部１４に、生体シミュレーション処理を実行する指示が入力された場合に実行される。図４に示すように、マクロ処理部１４ａ及びメゾ・ミクロ処理部１４ｂは、シミュレーション処理を実行する（Ｓ１０１）。そして、表示制御部１４ｃは、実行されたシミュレーション処理におけるタイムステップｔが、生体シミュレーションにおける最後のタイムステップであるか否かを判定する（Ｓ１０２）。最後のタイムステップでない場合（Ｓ１０２：Ｎｏ）には、表示制御部１４ｃは、タイムステップｔの値を１つインクリメントして、タイムステップを１つ進める（Ｓ１０３）。そして、マクロ処理部１４ａ及びメゾ・ミクロ処理部１４ｂは、Ｓ１０１に戻り、再び、シミュレーション処理を実行する。

一方、最後のタイムステップである場合（Ｓ１０２：Ｙｅｓ）には、表示制御部１４ｃは、次の処理を行う。表示制御部１４ｃは、シミュレーション処理において計算されたマクロモデル１３ａ、メゾモデル１３ｂ及びミクロモデル１３ｃの変位に応じた形状を１つめのタイムステップから最後のタイムステップまで順に表示するように表示部１２を制御する（Ｓ１０４）。

ここで、本実施例に係るシミュレーション処理の一例について説明する。本実施例では、マクロモデル１３ａとメゾモデル１３ｂとの間のマルチスケール解析は、各タイムステップでの非線形方程式を線形化し、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法を用いて反復的に解く問題（下記の式（１）〜（４））に帰着させる。

ここで、Ａ_Ｆ，ｅは、メゾモデル１３ｂの要素ｅに対応する剛性マトリクス（メゾ要素剛性マトリクス）である。Ａ_Ｃ，ｅは、マクロモデル１３ａの要素ｅに対応する剛性マトリクス（マクロ要素剛性マトリクス）である。ｆ_Ｆ，ｅは、Ａ_Ｆ，ｅに対応するメゾ等価節点力である。ｆ_Ｃ，ｅは、Ａ_Ｃ，ｅに対応するメゾ等価節点力である。ｘ_Ｆ，ｅは、メゾモデル１３ｂの要素ｅの解（メゾ解）、例えば、メゾモデル１３ｂの要素ｅの変位である。ｘ_Ｃは、マクロモデル１３ａの解（マクロ解）、例えば、マクロモデル１３ａの変位である。

は、マクロモデル１３ａの歪みからメゾ座標系における応力を求める作用素である。Ｂ_Ｃ，ｅは、マクロモデル１３ａの要素ｅにおける変位から歪みを求める作用素である。

は、メゾ周期変位への埋め込みを示す作用素である。Ｃ_Ｃは、マクロモデル１３ａに対応する剛性マトリクス（マクロ剛性マトリクス）である。ｇ_Ｃは、Ｃ_Ｃに対応するマクロ等価節点力である。

Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法を用いて反復して解く際には、下記の式（５）に従って、特性変位モードと称されるベクトルχ_Ｃ，ｅを計算する。

そして、メゾ解ｘ_Ｆ，ｅが互いに直接影響することがない性質を利用し、下記の式（６）にしたがって、シュアコンプリメント（Schur complement）Ｓ_Ｃ，ｅを計算する。

続いて、Ｓ_Ｃ，ｅを用いて、下記の式（７）に示すように、ブロックＬＵ（Lower triangle matrix Upper triangle matrix）分解を行う。このブロックＬＵ分解は、 正方行列を部分行列（ブロック）の集まりと捉え、下三角ブロック行列と上三角ブロック行列の積に分解する手法である。また、下三角ブロック行列は、対角ブロックより右上のブロックがすべてゼロ行列である行列をいい、上三角ブロック行列は、対角ブロックより左下のブロックがすべてゼロ行列である行列をいう。なお、ＬＵ分解は、すべてのブロック行列のサイズを1x1とした場合のブロックＬＵ分解であると見なすこともできる。

ここで、式（７）を用いて、式（４）を変形すると、下記の式（８）〜式（１１）を用いて、メゾ解ｘ_Ｆ，ｅを計算し、マクロ解ｘ_Ｃを計算し、メゾ解ｘ_Ｆ，ｅを更新することができる。すなわち、メゾ解ｘ_Ｆ，ｅの計算結果に基づいてマクロ解ｘ_Ｃを計算し、マクロ解ｘ_Ｃの計算結果をメゾ解ｘ_Ｆ，ｅに反映することができる。ここで、式（８）を用いて、メゾ解ｘ_Ｆ，ｅを計算する場合には、メゾモデルの要素ｅごとにメゾ解ｘ_Ｆ，ｅを独立して計算することができる。したがって、本実施例のように、各メゾ・ミクロ処理部１４ｂに各コアを割り当てて、各コアが、並列してメゾ解ｘ_Ｆ，ｅを計算する場合には、並列してメゾ解ｘ_Ｆ，ｅを計算するので、高速にメゾ解ｘ_Ｆ，ｅを計算することができる。よって、本実施例では、高効率なメゾ解ｘ_Ｆ，ｅの計算を実現することができる。

また、本実施例では、メゾモデル１３ｂとミクロモデル１３ｃとの間のマルチスケール解析においても、各タイムステップでの非線形方程式を線形化し、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法を用いて反復的に解く問題（下記の式（１２）〜（１５））に帰着させる。

ここで、Ａ_Ｍ，ｅは、ミクロモデル１３ｃの要素ｅに対応する剛性マトリクス（ミクロ要素剛性マトリクス）である。ｆ_Ｍ，ｅは、ミクロ等価節点力である。ｘ_Ｍ，ｅは、ミクロモデル１３ｃの要素ｅの解（ミクロ解）、例えば、メゾモデル１３ｂの要素ｅの変位である。Ｇは、メゾモデル１３ｂの歪みからミクロ座標系における応力を求める作用素である。Ｂ_Ｆ，ｅは、メゾモデル１３ｂの要素ｅにおける変位から歪みを求める作用素である。Ｐは、ミクロ周期変位への埋め込みを示す作用素である。Ｃ_Ｆは、メゾモデル１３ｂの剛性マトリクス（メゾ剛性マトリクス）である。ｇ_Ｆは、Ｃ_Ｆに対応するメゾ等価節点力である。

Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ法を用いて反復して解く際には、下記の式（１６）に従って、特性変位モードと称されるベクトルχ_Ｆ，ｅを計算する。

そして、ミクロ解ｘ_Ｍ，ｅが互いに直接影響することがない性質を利用し、下記の式（１７）にしたがって、シュアコンプリメントＳ_Ｆ，ｅを計算する。

続いて、Ｓ_Ｆ，ｅを用いて、下記の式（１８）に示すように、ＬＵ分解を行う。

ここで、式（１８）を用いて、式（１５）を変形すると、下記の式（１９）〜式（２２）を用いて、ミクロ解_Ｍ，ｅを計算し、メゾ解ｘ_Ｆ，ｅを計算し、ミクロ解_Ｍ，ｅを更新することができる。すなわち、ミクロ解_Ｍ，ｅの計算結果に基づいてメゾ解ｘ_Ｆ，ｅを計算し、メゾ解ｘ_Ｆ，ｅの計算結果をミクロ解_Ｍ，ｅに反映することができる。

ここで、マクロ処理部１４ａ及びメゾ・ミクロ処理部１４ｂが実行する先のシミュレーション処理の一例について説明する。図５は、実施例１に係るシミュレーション処理の手順を示すフローチャートである。なお、本実施例では、シミュレーション処理において、複数のマクロ処理部１４ａのそれぞれと、各マクロ処理部１４ａに対応する複数のメゾ・ミクロ処理部１４ｂとの間で情報のやりとりが発生する。しかしながら、説明の便宜上、図５の例では、１つのマクロ処理部１４ａが実行する処理と、１つのメゾ・ミクロ処理部１４ｂが実行する処理とを示す。

図５に示すように、マクロ処理部１４ａの計算部１５ａは、現在のタイムステップｔに対応する各心筋細胞のカルシウムイオン濃度を条件データ１３ｄから取得し、取得したカルシウムイオン濃度を各心筋細胞に対応する各メゾ・ミクロ処理部１４ｂに送信する（Ｓ２０１）。

そして、各メゾ・ミクロ処理部１４ｂの各第１の計算部１５ｂは、計算部１５ａから送信されたカルシウムイオン濃度を受信する（Ｓ２０２）。そして、各第１の計算部１５ｂは、モンテカルロ法を用いてミクロモデル１３ｃの収縮力などを計算する（Ｓ２０３）。一方、計算部１５ａは、各節点の圧力などを含む節点情報を各メゾ・ミクロ処理部１４ｂに送信する（Ｓ２０４）。

続いて、各第１の計算部１５ｂは、計算部１５ａから送信された節点情報を受信する（Ｓ２０５）。続いて、各第１の計算部１５ｂは、式（１６）を用いて、ベクトルχ_Ｆ，ｅを計算する（Ｓ２０６）。そして、Ｓ２０６において、各第１の計算部１５ｂは、式（１７）にしたがって、シュアコンプリメントＳ_Ｆ，ｅを計算する。続いて、Ｓ２０６において、各第１の計算部１５ｂは、Ｓ_Ｆ，ｅを用いて、ＬＵ分解を行って、式（１８）を得る。続いて、Ｓ２０６において、各第１の計算部１５ｂは、式（１８）を用いて、式（１５）を変形して、式（１９）〜式（２２）を得る。そして、Ｓ２０６において、各第１の計算部１５ｂは、式（２０）に従って、ｂ_Ｆ，ｅを計算する。

続いて、各第１の計算部１５ｂは、式（１９）に従って、ミクロ解ｘ_Ｍ，ｅを計算する（Ｓ２０７）。

そして、各メゾ・ミクロ処理部１４ｂの各第２の計算部１５ｃは、式（５）を用いて、ベクトルχ_Ｃ，ｅを計算する（Ｓ２０８）。そして、Ｓ２０８において、各第２の計算部１５ｃは、式（６）にしたがって、シュアコンプリメントＳ_Ｃ，ｅを計算する。続いて、Ｓ２０８において、各第２の計算部１５ｃは、Ｓ_Ｃ，ｅを用いて、ＬＵ分解を行って、式（７）を得る。続いて、Ｓ２０８において、各第２の計算部１５ｃは、式（７）を用いて、式（４）を変形して、式（８）〜式（９）、式（１１）を得る。そして、Ｓ２０８において、各第２の計算部１５ｃは、式（９）に従って、ｂ_Ｃ，ｅを計算する。

続いて、各第２の計算部１５ｃは、式（８）に従って、メゾモデル１３ｂの各要素のメゾ解ｘ_Ｆ，ｅを計算し、式（２１）に従って、メゾモデル１３ｂの各要素のメゾ解ｘ_Ｆ，ｅからメゾモデル１３ｂのメゾ解ｘ_Ｆを計算する（Ｓ２０９）。そして、各第２の計算部１５ｃは、メゾ解ｘ_Ｆが収束したか否かを判定する（Ｓ２１０）。収束しない場合（Ｓ２１０：Ｎｏ）には、Ｓ２０６に戻る。一方、収束した場合（Ｓ２１０：Ｙｅｓ）には、各第２の計算部１５ｃは、式（２）を計算した結果得られたＡ_Ｃ，ｅ、式（３）を計算した結果得られたｆ_Ｃ，ｅ、シュアコンプリメントＳ_Ｃ，ｅ、及び、ｂ_Ｃ，ｅをマクロ処理部１４ａに送信する（Ｓ２１１）。

マクロ処理部１４ａの計算部１５ａは、マクロ処理部１４ａに対応する全てのメゾ・ミクロ処理部１４ｂから、Ａ_Ｃ，ｅ、ｆ_Ｃ，ｅ、Ｓ_Ｃ，ｅ及びｂ_Ｃ，ｅを受信した場合に（Ｓ２１２）、式（７）を用いて、式（４）を変形して、式（１０）を得る（Ｓ２１３）。そして、Ｓ２１３において、計算部１５ａは、式（１０）に従って、マクロ解ｘ_Ｃを計算する。そして、計算部１５ａは、マクロ解ｘ_Ｃが収束したか否かを判定する（Ｓ２１４）。収束しない場合（Ｓ２１４：Ｎｏ）には、計算部１５ａは、マクロ解ｘ_Ｃを各メゾ・ミクロ処理部１４ｂに送信し（Ｓ２１５）、Ｓ２０４に戻って、各節点の圧力などを含む節点情報を更新し、更新後の節点情報を各メゾ・ミクロ処理部１４ｂに送信する。一方、収束した場合（Ｓ２１４：Ｙｅｓ）には、計算部１５ａは、収束したことを示すメッセージを各メゾ・ミクロ処理部１４ｂに送信し（Ｓ２１６）、処理結果を内部メモリに格納し、リターンする。

また、Ｓ２１１での処理を実行すると、各メゾ・ミクロ処理部１４ｂの各第１の更新部１５ｄは、収束したことを示すメッセージを受信したか否かを判定することで、マクロ解ｘ_Ｃが収束したか否かを判定する（Ｓ２１７）。ここで、第１の更新部１５ｄは、収束したことを示すメッセージを受信した場合には、マクロ解ｘ_Ｃが収束したと判定し、収束したことを示すメッセージを受信していない場合には、マクロ解ｘ_Ｃが収束していないと判定する。

収束していない場合（Ｓ２１７：Ｎｏ）には、第１の更新部１５ｄは、計算部１５ａからマクロ解ｘ_Ｃを受信すると（Ｓ２１８）、受信したマクロ解ｘ_Ｃを式（１１）に代入して、メゾ解ｘ_Ｆ，ｅを更新する（Ｓ２１９）。

続いて、各メゾ・ミクロ処理部１４ｂの各第２の更新部１５ｅは、メゾ解ｘ_Ｆを式（２２）に代入して、ミクロ解ｘ_Ｍ，ｅを更新し（Ｓ２２０）、処理結果を内部メモリに格納し、リターンする。

上述してきたように、本実施例に係る生体シミュレーション装置１０は、ミクロモデル１３ｃが示す、心臓の心筋細胞に含まれるサルコメアの挙動を計算する。そして、シミュレーション装置１０は、計算されたサルコメアの挙動に基づいて、メゾモデル１３ｂが示す心筋細胞の挙動を計算する。続いて、シミュレーション装置１０は、計算された心筋細胞の挙動に基づいて、マクロモデル１３ａが示す心臓の挙動を計算する。その後、シミュレーション装置１０は、計算された心臓の挙動に基づいて、メゾモデル１３ｂが示す心筋細胞の挙動を更新する。そして、シミュレーション装置１０は、更新された心筋細胞の挙動に基づいて、ミクロモデル１３ｃが示すサルコメアの挙動を更新する。このように、シミュレーション装置１０によれば、心臓全体を示すマクロモデル１３ａ、心筋細胞を示すメゾモデル１３ｂ、及び、サルコメアを示すミクロモデル１３ｃの３レベルのモデルを用いて、挙動のシミュレーションを行う。そのため、シミュレーション装置１０によれば、３レベルのマルチスケール解析という精緻なシミュレーションを行うことができる。したがって、シミュレーション装置１０によれば、精緻なシミュレーションの結果を得ることができる。

図６Ａ〜図６Ｃは、表示部に表示された、あるタイムステップにおける各モデルのシミュレーション結果の一例である。図６Ａは、表示部１２に表示された、あるタイムステップにおけるマクロモデル１３ａが示す心臓３０の形状を示す。また、図６Ｂは、表示部１２に表示された、あるタイムステップにおけるメゾモデル１３ｂが示す心筋細胞３１の形状を示す。また、図６Ｃは、表示部１２に表示された、あるタイムステップにおけるミクロモデル１３ｃが示すサルコメア３２の形状を示す。図６Ｃに示すように、本実施例に係るシミュレーション装置１０によれば、アクチン・フィラメント３２ａと、ミオシン・ヘッド３２ｂとのクロスブリッジ形成の様子が微細な部分の構造が可視化できる。

ここで、実施例１に係るマクロ処理部１４ａ及びメゾ・ミクロ処理部１４ｂの処理の実行タイミングの一例について説明する。図７は、実施例１に係るマクロ処理部及びメゾ・ミクロ処理部の処理の実行タイミングの一例を示すタイムチャートである。図７の例は、マクロ処理部１４ａが実行する各ステップＳでの処理時間及び処理が実行されるタイミング、及び、メゾ・ミクロ処理部１４ｂが実行する各ステップＳでの処理時間及び処理が実行されるタイミングを示す。図７の例に示すタイムチャートは、メゾ・ミクロ処理部１４ｂがＳ２０６〜Ｓ２１０での処理を実行している間、マクロ処理部１４ａは、処理を実行せず、待機中であることを示す。

そこで、実施例２では、マクロ処理部が待機中となる時間を抑制することができる生体シミュレーション装置について説明する。なお、実施例１と同一の構成については、同一の符号を付し、説明を省略する場合がある。図８は、実施例２に係る生体シミュレーション装置の機能的な構成の一例を示す図である。生体シミュレーション装置２０は、入力部１１、表示部１２、記憶部１３及び制御部２４を有する。

入力部１１、表示部１２及び記憶部１３については、実施例１と同様の構成であるので、説明を省略する。

制御部２４は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。図８に示すように、複数のマクロ処理部２４ａと、複数のメゾ・ミクロ処理部２４ｂと、表示制御部１４ｃとを有する。

複数のマクロ処理部２４ａのそれぞれは、複数のマクロモデル１３ａのそれぞれに対応し、各マクロ処理部２４ａは、対応するマクロモデル１３ａの挙動を計算する。また、複数のメゾ・ミクロ処理部２４ｂのそれぞれは、複数のメゾモデル１３ｂのそれぞれに対応し、各メゾ・ミクロ処理部２４ｂは、対応するメゾモデル１３ｂの挙動を計算する。さらに、各メゾ・ミクロ処理部２４ｂは、対応するメゾモデルが示す心筋細胞２１に含まれる各サルコメア２３を示す各ミクロモデル１３ｃの挙動を計算する。すなわち、本実施例では、１つのマクロモデル１３ａに対して、複数のメゾモデル１３ｂが対応するため、１つのマクロ処理部２４ａに対して複数のメゾ・ミクロ処理部２４ｂが対応することとなる。また、制御部２４は、マルチコアの複数のＣＰＵを有する。マクロ処理部２４ａ、メゾ・ミクロ処理部２４ｂ及び表示制御部１４ｃの各部は、後述の生体シミュレーションプログラムを各コアが実行することにより、各コア上で実現される。

マクロ処理部２４ａは、計算部２５ａを有する。メゾ・ミクロ処理部２４ｂは、第１の計算部２５ｂ、第２の計算部２５ｃ、第３の計算部２５ｄ、第１の更新部２５ｅ、第２の更新部２５ｆを有する。

ここで、実施例２に係る生体シミュレーション装置２０は、実施例１と同様に、先の図４の例のフローチャートが示す生体シミュレーション処理を実行するため、生体シミュレーション処理についての説明は省略する。実施例２に係るシミュレーション処理は、実施例１に係るシミュレーション処理と異なる。実施例２に係るシミュレーション処理の一例について説明する。実施例２に係るシミュレーション処理と、実施例１に係るシミュレーション処理との異なる点について説明する。実施例２では、実施例１の式（７）に代えて、以下の式（２３）を用いる。式（２３）は、近似的にＬＵ分解を行った結果の式である。

また、実施例１では、Ｓ２０８において、式（９）を用いてｂ_Ｃ，ｅを計算したが、実施例２では、ｂ_Ｃ，ｅを計算しない。また、実施例１では、Ｓ２１３において、式（１０）を用いて、マクロ解ｘ_Ｃを計算したが、実施例２では、式（１０）に代えて、以下の式（２４）を用いて、マクロ解ｘ_Ｃを計算する。

ここで、第２の演算処理部であるマクロ処理部２４ａ及び第１の演算処理部であるｓメゾ・ミクロ処理部２４ｂが実行する先のシミュレーション処理の一例について説明する。図９は、実施例２に係るシミュレーション処理の手順を示すフローチャートである。なお、本実施例においても、シミュレーション処理で、複数のマクロ処理部２４ａのそれぞれと、各マクロ処理部２４ａに対応する複数のメゾ・ミクロ処理部２４ｂとの間で情報のやりとりが発生する。しかしながら、説明の便宜上、図９の例では、１つのマクロ処理部２４ａが実行する処理と、１つのメゾ・ミクロ処理部２４ｂが実行する処理とを示す。

図９に示すように、マクロ処理部２４ａの計算部２５ａは、現在のタイムステップｔに対応する各心筋細胞のカルシウムイオン濃度を条件データ１３ｄから取得し、取得したカルシウムイオン濃度を各心筋細胞に対応する各メゾ・ミクロ処理部２４ｂに送信する（Ｓ３０１）。

そして、各メゾ・ミクロ処理部２４ｂの第１の計算部２５ｂは、計算部２５ａから送信されたカルシウムイオン濃度を受信する（Ｓ３０２）。そして、各第１の計算部２５ｂは、モンテカルロ法を用いてミクロモデル１３ｃの収縮力などを計算する（Ｓ３０３）。一方、計算部２５ａは、各節点の圧力などを含む節点情報を各メゾ・ミクロ処理部２４ｂに送信する（Ｓ３０４）。

続いて、各第１の計算部２５ｂは、計算部２５ａから送信された節点情報を受信する（Ｓ３０５）。続いて、各第１の計算部２５ｂは、式（１６）を用いて、ベクトルχ_Ｆ，ｅを計算する（Ｓ３０６）。そして、Ｓ３０６において、各第１の計算部２５ｂは、式（１７）にしたがって、シュアコンプリメントＳ_Ｆ，ｅを計算する。続いて、Ｓ３０６において、各第１の計算部２５ｂは、Ｓ_Ｆ，ｅを用いて、ＬＵ分解を行って、式（１８）を得る。続いて、Ｓ３０６において、各第１の計算部２５ｂは、式（１８）を用いて、式（１５）を変形して、式（１９）〜式（２２）を得る。そして、Ｓ３０６において、各第１の計算部２５ｂは、式（２０）に従って、ｂ_Ｆ，ｅを計算する。

続いて、各第１の計算部２５ｂは、式（１９）に従って、ミクロ解ｘ_Ｍ，ｅを計算する（Ｓ３０７）。そして、Ｓ３０７において、各メゾ・ミクロ処理部２４ｂの各第２の計算部２５ｃは、式（５）を用いて、ベクトルχ_Ｃ，ｅを計算する。そして、Ｓ３０７において、各第２の計算部２５ｃは、式（６）にしたがって、シュアコンプリメントＳ_Ｃ，ｅを計算する。次に、Ｓ３０７において、各第２の計算部２５ｃは、Ｓ_Ｃ，ｅを用いて、ＬＵ分解を行って、式（２３）を得る。続いて、Ｓ３０７において、各メゾ・ミクロ処理部２４ｂの各第３の計算部２５ｄは、式（３）を用いて、ｆ_Ｃ，ｅを計算する。

そして、各第３の計算部２５ｄは、式（２）を計算した結果得られたＡ_Ｃ，ｅ、式（３）を計算した結果得られたｆ_Ｃ，ｅ、及び、シュアコンプリメントＳ_Ｃ，ｅをマクロ処理部２４ａに送信する（Ｓ３０８）。

マクロ処理部２４ａの計算部２５ａは、マクロ処理部２４ａに対応する全てのメゾ・ミクロ処理部２４ｂから、Ａ_Ｃ，ｅ、ｆ_Ｃ，ｅ及びＳ_Ｃ，ｅを受信した場合に（Ｓ３０９）、式（２３）を用いて、式（４）を変形して、式（２４）を得る（Ｓ３１０）。そして、Ｓ３１０において、計算部２５ａは、式（２４）に従って、マクロ解ｘ_Ｃを計算する。そして、計算部２５ａは、マクロ解ｘ_Ｃが収束したか否かを判定する（Ｓ３１１）。収束しない場合（Ｓ３１１：Ｎｏ）には、計算部２５ａは、マクロ解ｘ_Ｃを各メゾ・ミクロ処理部２４ｂに送信し（Ｓ３１２）、Ｓ３０４に戻って、各節点の圧力などを含む節点情報を更新し、更新後の節点情報を各メゾ・ミクロ処理部２４ｂに送信する。一方、収束した場合（Ｓ３１１：Ｙｅｓ）には、計算部２５ａは、処理結果を内部メモリに格納し、リターンする。

また、Ｓ３０８での処理を実行すると、各第２の計算部２５ｃは、式（２３）を用いて、式（４）を変形して、式（８）、式（１１）を得る（Ｓ３１３）。

続いて、各第２の計算部２５ｃは、式（８）に従って、メゾモデル１３ｂの各要素のメゾ解ｘ_Ｆ，ｅを計算し、式（２１）に従って、メゾモデル１３ｂの各要素のメゾ解ｘ_Ｆ，ｅからメゾモデル１３ｂのメゾ解ｘ_Ｆを計算する（Ｓ３１４）。そして、各第２の計算部２５ｃは、メゾ解ｘ_Ｆが収束したか否かを判定する（Ｓ３１５）。収束しない場合（Ｓ３１５：Ｎｏ）には、各メゾ・ミクロ処理部２４ｂの各第１の更新部２５ｅは、計算部２５ａからマクロ解ｘ_Ｃを受信すると（Ｓ３１６）、受信したマクロ解ｘ_Ｃを式（１１）に代入して、メゾ解ｘ_Ｆ，ｅを更新する（Ｓ３１７）。

続いて、各メゾ・ミクロ処理部２４ｂの各第２の更新部２５ｆは、メゾ解ｘ_Ｆを式（２２）に代入して、ミクロ解ｘ_Ｍ，ｅを更新し（Ｓ３１８）、処理結果を内部メモリに格納し、リターンする。

次に、実施例２に係るマクロ処理部２４ａ及びメゾ・ミクロ処理部２４ｂの処理の実行タイミングの一例について説明する。図１０は、実施例２に係るマクロ処理部及びメゾ・ミクロ処理部の処理の実行タイミングの一例を示すタイムチャートである。図１０の例は、マクロ処理部２４ａが実行する各ステップＳでの処理時間及び処理が実行されるタイミング、及び、メゾ・ミクロ処理部２４ｂが実行する各ステップＳでの処理時間及び処理が実行されるタイミングを示す。本実施例では、実施例１の式（７）に代えて、以下の式（２３）を用いて式（４）を変形して、式（２４）を得る。そして、本実施例では、式（２４）を用いて、マクロ解ｘ_Ｃを計算する。そのため、図１０の例に示すように、メゾ・ミクロ処理部２４ｂがＳ２０６〜Ｓ２１０での処理を実行している間、マクロ処理部２４ａは、Ｓ３１０及びＳ３１１での処理を実行することができる。すなわち、本実施例では、メゾ・ミクロ処理部２４ｂによるメゾ解ｘ_Ｆを計算する処理と、マクロ処理部２４ａによるマクロ解ｘ_Ｃを計算する処理とを並行して実行することができる。したがって、本実施例に係るシミュレーション装置２０によれば、計算時間の増大を抑制することができる。また、本実施例に係るシミュレーション装置２０によれば、ＣＰＵのコアの利用率の低下を抑制することができる。

上述してきたように、本実施例に係る生体シミュレーション装置２０は、ミクロモデル１３ｃが示す、心臓の心筋細胞に含まれるサルコメアの挙動を計算する。続いて、シミュレーション装置２０は、マクロモデル１３ａが示す心臓に加えられる力を計算する。そして、シミュレーション装置２０は、計算された力に基づいて、心臓の挙動を計算する。続いて、シミュレーション装置２０は、計算されたサルコメアの挙動に基づいて、メゾモデル１３ｂが示す心筋細胞の挙動を計算する。その後、シミュレーション装置２０は、計算された心臓の挙動に基づいて、メゾモデル１３ｂが示す心筋細胞の挙動を更新する。そして、シミュレーション装置２０は、更新された心筋細胞の挙動に基づいて、ミクロモデル１３ｃが示すサルコメアの挙動を更新する。このように、シミュレーション装置２０によれば、心臓全体を示すマクロモデル１３ａ、心筋細胞を示すメゾモデル１３ｂ、及び、サルコメアを示すミクロモデル１３ｃの３レベルのモデルを用いて、挙動のシミュレーションを行う。そのため、シミュレーション装置２０によれば、３レベルのマルチスケール解析という精緻なシミュレーションを行うことができる。したがって、シミュレーション装置２０によれば、精緻なシミュレーションの結果を得ることができる。

また、上述したように、本実施例に係るシミュレーション装置２０によれば、メゾ・ミクロ処理部２４ｂがＳ２０６〜Ｓ２１０での処理を実行している間、マクロ処理部２４ａは、Ｓ３１０及びＳ３１１での処理を実行することができる。すなわち、本実施例では、メゾ・ミクロ処理部２４ｂによるメゾ解ｘ_Ｆを計算する処理と、マクロ処理部２４ａによるマクロ解ｘ_Ｃを計算する処理とを並行して実行することができる。したがって、本実施例に係るシミュレーション装置２０によれば、計算時間の増大を抑制することができる。また、本実施例に係るシミュレーション装置２０によれば、ＣＰＵのコアの利用率の低下を抑制することができる。

さて、これまで開示の装置に関する実施例について説明したが、本発明は上述した実施例以外にも、種々の異なる形態にて実施されてよいものである。そこで、以下では、本発明に含まれる他の実施例を説明する。

例えば、上述した実施例１及び実施例２では、１つのメゾ・ミクロ処理部１４ａ、２４ａが１つのコア上で実現され、１つのマクロ処理部１４ｂ、２４ｂが１つのコア上で実現される場合について説明した。しかしながら、開示の装置は、これに限られない。例えば、複数のメゾ・ミクロ処理部１４ａ、２４ａが１つのコア上で実現されてもよい。また、複数のマクロ処理部１４ｂ、２４ｂが１つのコア上で実現されてもよい。更には、メゾ・ミクロ処理部１４ａ、２４ａ及びマクロ処理部１４ｂ、２４ｂが１つのコア上で実現されてもよい。

また、実施例において説明した各処理のうち、自動的に行われるものとして説明した処理の全部または一部を手動的に行うこともできる。また、本実施例において説明した各処理のうち、手動的に行われるものとして説明した処理の全部または一部を公知の方法で自動的に行うこともできる。

また、各種の負荷や使用状況などに応じて、実施例において説明した各処理の各ステップでの処理を任意に細かくわけたり、あるいはまとめたりすることができる。また、ステップを省略することもできる。

また、各種の負荷や使用状況などに応じて、実施例において説明した各処理の各ステップでの処理の順番を変更できる。

また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的状態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。

［生体シミュレーションプログラム］
また、上記の実施例で説明した生体シミュレーション装置（１０、２０）の各種の処理は、あらかじめ用意されたプログラムをパーソナルコンピュータやワークステーションなどのコンピュータシステムで実行することによって実現することもできる。そこで、以下では、図１１を用いて、上記の実施例で説明した生体シミュレーション装置と同様の機能を有するプログラムを実行するコンピュータの一例を説明する。図１１は、生体シミュレーションプログラムを実行するコンピュータを示す図である。

図１１に示すように、コンピュータ３００は、マルチコアの複数のＣＰＵ３１０、ＲＯＭ３２０、Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ（ＨＤＤ）３３０、ＲＡＭ３４０を有する。これら３１０〜３４０は、バス３５０を介して接続される。

ＲＯＭ３２０には、ＯＳなどの基本プログラムが記憶されている。また、ＨＤＤ３３０には、上記の実施例１で示す計算部１５ａ、第１の計算部１５ｂ、第２の計算部１５ｃ、第１の更新部１５ｄ、第２の更新部１５ｅと同様の機能を発揮する生体シミュレーションプログラム３３０ａが予め記憶される。なお、生体シミュレーションプログラム３３０ａについては、適宜分離しても良い。また、生体シミュレーションプログラム３３０ａは、上記の実施例２で示す計算部２５ａ、第１の計算部２５ｂ、第２の計算部２５ｃ、第３の計算部２５ｄ、第１の更新部２５ｅ、第２の更新部２５ｆと同様の機能を発揮するプログラムであってもよい。

そして、ＣＰＵ３１０が、生体シミュレーションプログラム３３０ａを、ＨＤＤ３３０から読み出して実行する。

なお、上記した生体シミュレーションプログラム３３０ａについては、必ずしも最初からＨＤＤ３３０に記憶させておく必要はない。

例えば、コンピュータ３００に挿入されるフレキシブルディスク（ＦＤ）、ＣＤ−ＲＯＭ、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの「可搬用の物理媒体」に生体シミュレーションプログラム３３０ａを記憶させておく。そして、コンピュータ３００がこれらから生体シミュレーションプログラム３３０ａを読み出して実行するようにしてもよい。

さらには、公衆回線、インターネット、ＬＡＮ、ＷＡＮなどを介してコンピュータ３００に接続される「他のコンピュータ（またはサーバ）」などに生体シミュレーションプログラム３３０ａを記憶させておく。そして、コンピュータ３００がこれらから生体シミュレーションプログラム３３０ａを読み出して実行するようにしてもよい。

１０ 生体シミュレーション装置
１３ａ マクロモデル
１３ｂ メゾモデル
１３ｃ ミクロモデル
１４ａ マクロ処理部
１４ｂ メゾ・ミクロ処理部
１５ａ 計算部
１５ｂ 第１の計算部
１５ｃ 第２の計算部
１５ｄ 第１の更新部
１５ｅ 第２の更新部

Claims (8)

1. コンピュータに、
   生体の心臓の細胞に含まれるサルコメアのモデルの収縮力を算出し、算出した前記収縮力と、前記細胞のモデルにかかる圧力を含む情報とに基づいて、前記サルコメアのモデルの変位を算出し、
   算出した前記サルコメアのモデルの変位に基づいて、前記細胞のモデルの変位を算出するとともに、算出した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新し、
   算出した前記細胞のモデルの変位に基づいて、前記心臓のモデルの変位を算出し、
   算出した前記心臓のモデルの変位に基づいて、前記細胞のモデルの変位を更新し、更新した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新する、
   処理を実行させることを特徴とする生体シミュレーションプログラム。
2. 前記心臓のモデルの変位を算出する処理は、算出した前記心臓のモデルの変位に基づいて、前記細胞のモデルにかかる圧力を含む情報を更新する、
   ことを特徴する請求項１に記載の生体シミュレーションプログラム。
3. コンピュータに、
   生体の心臓の細胞に含まれるサルコメアのモデルの収縮力を算出し、算出した前記収縮力と、前記細胞のモデルにかかる圧力を含む情報とに基づいて、前記サルコメアのモデルの変位を算出し、
   算出した前記サルコメアのモデルの変位に基づいて、前記心臓のモデルの変位を算出し、
   算出した前記サルコメアのモデルの変位に基づいて、前記細胞のモデルの変位を算出し、
   算出した前記心臓のモデルの変位に基づいて、前記細胞のモデルの変位を更新し、更新した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新する、
   処理を実行させることを特徴とする生体シミュレーションプログラム。
4. 前記心臓のモデルの変位を算出する処理と、前記細胞のモデルの変位を算出する処理とは、並列に実行させる、
   ことを特徴する請求項３に記載の生体シミュレーションプログラム。
5. コンピュータが、
   生体の心臓の細胞に含まれるサルコメアのモデルの収縮力を算出し、算出した前記収縮力と、前記細胞のモデルにかかる圧力を含む情報とに基づいて、前記サルコメアのモデルの変位を算出し、
   算出した前記サルコメアのモデルの変位に基づいて、前記細胞のモデルの変位を算出するとともに、算出した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新し、
   算出した前記細胞のモデルの変位に基づいて、前記心臓のモデルの変位を算出し、
   算出した前記心臓のモデルの変位に基づいて、前記細胞のモデルの変位を更新し、更新した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新する、
   処理を実行することを特徴とする生体シミュレーション方法。
6. コンピュータが、
   生体の心臓の細胞に含まれるサルコメアのモデルの収縮力を算出し、算出した前記収縮力と、前記細胞のモデルにかかる圧力を含む情報とに基づいて、前記サルコメアのモデルの変位を算出し、
   算出した前記サルコメアのモデルの変位に基づいて、前記心臓のモデルの変位を算出し、
   算出した前記サルコメアのモデルの変位に基づいて、前記細胞のモデルの変位を算出し、
   算出した前記心臓のモデルの変位に基づいて、前記細胞のモデルの変位を更新し、更新した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新する、
   処理を実行することを特徴とする生体シミュレーション方法。
7. 生体の心臓の細胞に含まれるサルコメアのモデルの収縮力を算出し、算出した前記収縮力と、前記細胞のモデルにかかる圧力を含む情報とに基づいて、前記サルコメアのモデルの変位を算出する第１算出部と、
   算出した前記サルコメアのモデルの変位に基づいて、前記細胞のモデルの変位を算出するとともに、算出した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新する第１更新部と、
   算出した前記細胞のモデルの変位に基づいて、前記心臓のモデルの変位を算出する第２算出部と、
   算出した前記心臓のモデルの変位に基づいて、前記細胞のモデルの変位を更新し、更新した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新する第２更新部と、
   を有することを特徴とする生体シミュレーション装置。
8. 生体の心臓の細胞に含まれるサルコメアのモデルの収縮力を算出し、算出した前記収縮力と、前記細胞のモデルにかかる圧力を含む情報とに基づいて、前記サルコメアのモデルの変位を算出する第１算出部と、
   算出した前記サルコメアのモデルの変位に基づいて、前記心臓のモデルの変位を算出する第２算出部と、
   算出した前記サルコメアのモデルの変位に基づいて、前記細胞のモデルの変位を算出する第３算出部と、
   算出した前記心臓のモデルの変位に基づいて、前記細胞のモデルの変位を更新し、更新した前記細胞のモデルの変位に基づいて、前記サルコメアのモデルの変位を更新する更新部と、
   を有することを特徴とする生体シミュレーション装置。

Images