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

Description

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

分子生物学の分野では、実験事実に基づきサルコメア内のミオシン・アクチン間の架橋運動を記述した多様なサルコメアモデルが提案されている。例えば、サルコメア内のミオシン・アクチン間の結合に寄与する分子に対して複数の代表的な状態が定義され、近傍分子間の相互作用および分子のエネルギーなどが考慮されて状態間の遷移率が定義されたサルコメアモデルがある。かかるサルコメア内の分子の運動と筋肉の連続体モデルを連成させた解析を行う際には、分子モデルに対してはその平均的な挙動を想定し状態濃度を変数とする常微分方程式を導く方法が一般的である。また、サルコメアモデルの計算結果から連続体モデルでの収縮力を算出し、かかる結果をもとに連続体としての筋肉の運動を計算する方法がとられることが一般的である。

Peterson J.N., Hunter W.C., Berman M.R.: Estimated time course of Ca2+bound to troponin C during relaxation in isolated cardiac muscle, American Physiological Society, H1013-H1024(1991) Hunter P.J., McCulloch A.D., ter Keurs H.E.D.J.: Modeling the mechanical properties of cardiac muscle, Progress in Biophysics & Molecular Biology 69, 289-331(1998) Razumova M.V., Bukatina A.E., Campbell K.B.: Stiffness-distortion sarcomere model for muscle simulation. J. Appl. Physiol. 87, 1861-1876(1999) Rice J.J., Wang F., Bers D.M., de Tombe P.P.: Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations. Biophys. J. 95, 2368-2390(2008)

しかしながら、上述したサルコメアモデルは、サルコメアモデル内の分子モデルの平均化した挙動を表すため、かかるサルコメアモデルを用いてサルコメアの挙動のシミュレーションを行う場合には、分子モデルの状態遷移に対するシミュレーションの結果が精緻さに欠け、結果として連続体としての筋肉の運動に対しても精緻さに欠けるという問題点がある。さらに、筋肉の運動からサルコメアモデル内の分子の状態遷移への強いフィードバックもあるので、サルコメアモデルから連続体モデルへの一方向の情報伝達に基づくシミュレーションでは、安定した解析を実行することが困難である。

１つの側面では、確率的な挙動を有する分子モデルと連続体としての筋肉の運動の精緻な連成シミュレーションの結果を得ることを目的とする。他の側面では、筋肉の連続体モデルの運動から分子モデルへのフィードバックを陰的解法に基づき連成シミュレーションに取り込むことで、かかる連成解析を適切な時間刻み幅で安定に遂行することを目的とする。

本願の開示する生体シミュレーションプログラムは、コンピュータに、生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、複数の状態間の遷移率が定義されたモデルを用いて、次の処理を実行させる。すなわち、生体シミュレーションプログラムは、コンピュータに、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態を計算する処理を実行させる。そして、生体シミュレーションプログラムは、コンピュータに、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態に基づいて、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動を計算する処理を実行させる。そして、生体シミュレーションプログラムは、コンピュータに、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動に基づいて、サルコメアの挙動を計算する処理を実行させる。そして、生体シミュレーションプログラムは、コンピュータに、サルコメアの挙動に基づいて、筋肉の挙動を計算する処理を実行させる。

確率的な挙動を有する分子モデルと連続体としての筋肉の運動の精緻な連成シミュレーションの結果を得ることができる。

図１は、実施例に係る生体シミュレーション装置の機能的な構成の一例を示す図である。 図２は、心筋モデルの一例の一部を示す図である。 図３は、サルコメアモデルの一例を示す図である。 図４は、Ｔ／Ｔユニットの状態遷移の一例を示す図である。 図５は、ミオシンヘッドの状態遷移の一例を示す図である。 図６は、ＳＬに応じてアクチンフィラメント２１のオーバラップ状態が変化することを説明するための図である。 図７は、実施例に係る生体シミュレーション処理の手順を示すフローチャートである。 図８は、実施例に係るシミュレーション処理の手順を示すフローチャートである。 図９は、ミオシンアームの伸びについての計算方法の一例について説明するための図である。 図１０は、実施例に係るモンテカルロシミュレーション処理の手順を示すフローチャートである。 図１１は、実施例に係る有限要素解析の手順を示すフローチャートである。 図１２は、実施例に係る受付部が受け付けるＴ／Ｔユニットの遷移情報を説明するための図である。 図１３は、実施例に係る受付部が受け付けるミオシンヘッドの情報を説明するための図である。 図１４は、実施例に係る受付部が受け付けるＴ／Ｔユニットの状態に関わる関数情報を説明するための図である。 図１５は、実施例に係る受付部が受け付けるミオシンヘッドの状態に関わる関数情報を説明するための図である。 図１６は、実施例に係る受付部が受け付けるミオシンヘッドの非結合-結合状態間の遷移を定義する情報を説明するための図である。 図１７は、実施例に係る受付部が受け付けるミオシンヘッドの首振り前後の状態間の遷移を定義する情報を説明するための図である。 図１８は、実施例に係る受付部が受け付ける結合解離の遷移を定義する情報を説明するための図である。 図１９は、ミオシンヘッドの状態と状態間遷移の定義をもとにモンテカルロステップを実行するコード（モンテカルロコード）を自動生成した結果を示す図である。 図２０Ａは、心拍動のシミュレーションを実行するときの心筋の連続体モデルに関わるパラメータ指定の一例を示す図である。 図２０Ｂは、心拍動のシミュレーションを実行するときの線維方向を指定する方法についての一例を示す図である。 図２１Ａは、心拍動性能に関わる緒量のシミュレーション結果の一例を示す図である。 図２１Ｂは、左心室圧容積変化（上）と仕事率及びエネルギー消費率（下）のシミュレーション結果の一例を示す図である。 図２１Ｃは、収縮力の分布に関するシミュレーション結果の一例を示す図である。 図２２は、生体シミュレーションプログラムを実行するコンピュータを示す図である。

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

実施例に係る生体シミュレーション装置について説明する。本実施例に係る生体シミュレーション装置は、生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、複数の状態間の遷移率が定義されたモデルと有限要素メッシュにより表現された筋肉の連続体モデルを用いて、次の処理を行う。すなわち、生体シミュレーション装置は、筋肉の連続体モデルの各有限要素に埋め込まれた複数のアクチン及び複数のミオシンの状態遷移を計算する。また、生体シミュレーション装置は、分子モデルの状態遷移と連続体として筋肉の運動から導かれる個々の分子が生じる収縮力と分子モデルが埋め込まれた有限要素の連続体モデルでの収縮力を整合させた運動方程式を解き、次のような解析を行う。すなわち、生体シミュレーション装置は、分子の状態遷移運動と筋肉の運動を連成させた解析を行う。図１は、実施例に係る生体シミュレーション装置の機能的な構成の一例を示す図である。図１の例に示すように、生体シミュレーション装置１０は、入力部１１、表示部１２、記憶部１３及び制御部１４を有する。

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

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

記憶部１３は、制御部１４で実行される各種プログラムを記憶する。また、記憶部１３は、心筋モデル１３ａ、サルコメアモデル１３ｂを記憶する。

心筋モデル１３ａは、連続体である心臓全体の筋肉のモデルを複数個の要素に分割したモデルである。例えば、心臓全体の筋肉のモデルを、左心房、左心室、右心房、右心室の４個のモデルに分割し、左心房のモデル、左心室のモデル、右心房のモデル、右心室のモデルのそれぞれを心筋モデル１３ａとすることができる。心筋モデル１３ａは、有限要素解析によって、心筋モデル１３ａが示す筋肉のモデルの挙動を計算する際に用いられる。図２は、心筋モデルの一例の一部を示す図である。図２の例では、心筋モデル１３ａは左心室のモデルであり、左心室のモデルの一部が示されている。また、図２の例では、心筋モデル１３ａの有限要素メッシュ内の複数の要素１３ａ＿１のそれぞれの線維方向が矢印で示されている。以下の説明では、かかる線維方向のベクトルをｆで表す場合がある。

また、各要素１３ａ＿１には、複数のサルコメアモデル１３ｂが埋め込まれている。以下、各要素１３ａ＿１に、「ｎｓ」個のサルコメアモデル１３ｂが埋め込まれている場合について説明する。サルコメアモデル１３ｂは、構成要素である分子が確率的な挙動を示すモデルである。サルコメアモデル１３ｂは、いわゆる協調性、及び、いわゆるサルコメア長依存性などを有するように定義される。図３は、サルコメアモデルの一例を示す図である。図３の例に示すように、サルコメアモデル１３ｂは、アクチンフィラメント２１上の複数のＴ／Ｔユニット（トロポニン／トロポミオシンのユニット）２０と、ミオシンフィラメント２３上の複数のミオシンヘッド２２とを有する。なお、後述するが、サルコメアモデル１３ｂには、ミオシンフィラメント２３とミオシンヘッド２２とをつないでいるミオシンアーム２４も含まれる。また、以下、１つのサルコメアモデル１３ｂ内に、「ｎｍ」個のミオシンヘッド２２及びミオシンアーム２４の組が存在する場合について説明する。

サルコメアモデル１３ｂにおいては、Ｔ／Ｔユニット２０及びミオシンヘッド２２に対して複数の状態が定義されており、予め定められた状態間の遷移率に従いモンテカルロシミュレーションが実行される。すなわち、モンテカルロシミュレーションにより、Ｔ／Ｔユニット２０及びミオシンヘッド２２の確率的な挙動が計算される。

Ｔ／Ｔユニット２０に対しては、直下のミオシンヘッド２２の状態に応じて状態の遷移率が制御される。図４は、Ｔ／Ｔユニットの状態遷移の一例を示す図である。図４の例に示すように、本実施例では、Ｔ／Ｔユニット２０にカルシウムイオン（Ｃａ^２＋）が結合している状態（結合状態：Ｃａ−ｏｎ）と、結合していない状態（非結合状態：Ｃａ−ｏｆｆ）との２種類の状態がある。図４の例に示す状態遷移において、例えば、カルシウムイオン濃度が高くなるほど、結合状態に遷移しやすくなる。また、カルシウムイオンの解離確率Ｋ_ｏｆｆが直下のミオシンヘッド２２の状態に応じて定められる。また、Ｋ_ｏｎ、Ｋ_ｏｆｆは、当該Ｔ／Ｔユニット２０の配下にあるアクチン切片にミオシンヘッド２２が結合しているかいないかによって、値が変化する。例えば、Ｔ／Ｔユニット２０の配下にあるアクチン切片にミオシンヘッド２２が結合している場合には、Ｋ_ｏｆｆが小さくなり、Ｔ／Ｔユニット２０からカルシウムイオンが外れにくくなる。すなわち、本実施例のサルコメアモデル２０のＴ／Ｔユニット２０は、当該Ｔ／Ｔユニット２０の配下にあるアクチン切片にミオシンヘッド２２が結合している場合には、Ｔ／Ｔユニット２０からカルシウムイオンが外れにくくなるという協調性を有するモデルである。ここで、Ｔ／Ｔユニット２０にカルシウムイオンが結合すると、ミオシンヘッド２２との結合部位を隠している紐状のタンパク質であるトロポミオシンの位置がずれ、Ｔ／Ｔユニット２０の配下にあるアクチン切片とミオシンヘッド２２との結合がより容易になる。そのため、Ｔ／Ｔユニット２０にカルシウムイオンが結合すると、Ｔ／Ｔユニット２０の配下にあるアクチン切片とミオシンヘッド２２とが結合する確率が高くなる。そして、Ｔ／Ｔユニット２０の配下にあるアクチン切片とミオシンヘッド２２とが結合すると、さらに、トロポミオシンの位置がずれて、当該ミオシンの近傍のミオシンヘッド２２が、直上のアクチン切片と結合する確率が高くなる。なお、Ｔ／Ｔユニット２０の配下にあるアクチン切片とは図３の２０の線で区切られたものを指し、近傍のミオシンヘッド２２とは、あるミオシンヘッド２２から所定範囲に存在するミオシンヘッド２２のことを指す。

また、ミオシンヘッド２２に対しては、近傍のミオシンヘッド２２および直上のアクチン切片をその配下とするＴ／Ｔユニット２０の状態に応じて状態遷移が制御される。図５は、ミオシンヘッドの状態遷移の一例を示す図である。図５の例に示すように、本実施例では、Ｎ_ＸＢと、Ｐ_ＸＢと、ＸＢ_ＰｒｅＲと、ＸＢ_{ＰｏｓｔＲ}との４種類の状態がある。Ｎ_ＸＢは、ミオシンヘッド２２が、Ｔ／Ｔユニット２０の配下にあるアクチン切片に結合していない状態（非結合状態）である。また、Ｐ_ＸＢは、ミオシンヘッド２２が、Ｔ／Ｔユニット２０の配下にあるアクチン切片との結合を開始した状態（弱結合状態）である。また、ＸＢ_ＰｒｅＲは、Ｐ_ＸＢから、ミオシンヘッド２２に化学状態の変化が発生して、より強く、ミオシンヘッド２２が、Ｔ／Ｔユニット２０の配下にあるアクチン切片に結合した状態（強結合（１））である。また、ＸＢ_{ＰｏｓｔＲ}は、ミオシンヘッド２２が、Ｔ／Ｔユニット２０の配下にあるアクチン切片との結合を維持した状態で、ＸＢ_ＰｒｅＲから、ミオシンヘッド２２の回転運動（首振り運動）が発生した状態（強結合（２））である。図５中、ｋ_ｎｐは、Ｎ_ＸＢからＰ_ＸＢへ遷移する際の係数である。サルコメアモデル１３ｂのアクチンフィラメント２１のオーバラップが大きくなると、ミオシンヘッド２２と結合することができるアクチンフィラメントの範囲が短くなる。このため、サルコメアモデル１３ｂのアクチンフィラメント２１のオーバラップが大きくなると、ミオシンヘッド２２がアクチンに結合しにくくなる。本実施例では、サルコメアモデル１３ｂのアクチンフィラメント２１のオーバラップが大きくなると、ｋ_ｎｐの値は小さくなり、ミオシンヘッド２２がアクチンに結合する確率を低くする。

また、ｋ_ｐｎは、Ｐ_ＸＢ、ＸＢ_ＰｒｅＲ、ＸＢ_{ＰｏｓｔＲ}からＮ_ＸＢへ遷移する際の係数である。また、γ^ｎ及びγ^−ｎは、ミオシンヘッド２２の非結合-結合状態間の遷移に関わる協調性を示す。ここで、ｎは、隣接し、かつ、アクチンに結合しているミオシンヘッド２２の数（０〜２）を示す。例えば、γ＝８０の場合、隣接し、かつ、アクチンに結合しているミオシンヘッド２２の数が「１」である場合、γ^ｎは「８０」となる。すなわち、図５の例に示す場合において、８０倍結合しやすくなることを示す。また、γ＝８０の場合、隣接し、かつ、アクチンに結合しているミオシンヘッド２２の数が「２」である場合、γ^ｎは、「６４００」となる。すなわち、図５の例に示す場合において、６４００倍結合しやすくなることを示す。また、γ＝８０の場合、隣接し、かつ、アクチンに結合しているミオシンヘッド２２の数が「１」である場合、γ^−ｎは「１／８０」となる。すなわち、図５の例に示す場合において、解離のしやすさが１／８０倍となることを示す。また、γ＝８０の場合、隣接し、かつ、アクチンに結合しているミオシンヘッド２２の数が「２」である場合、γ^−ｎは「１／６４００」となる。すなわち、図５の例に示す場合において、解離のしやすさが１／６４００倍となることを示す。これにより、あるミオシンヘッド２２について、自身がアクチンに結合していないが、隣接するミオシンヘッド２２がアクチンに結合している場合には、アクチンに結合しやすくなる。また、あるミオシンヘッド２２について、自身がアクチンに結合しており、隣接するミオシンヘッド２２もアクチンに結合している場合には、解離しにくくなる。

また、図５中、ミオシンヘッド２２がＮ_ＸＢ、Ｐ_ＸＢ、ＸＢ_ＰｒｅＲの状態にある場合、ミオシンヘッド２２には、ＡＴＰ（アデノシン三リン酸）が結合している。また、ＡＴＰ→ＡＤＰ（アデノシン二リン酸）＋ｐｉ（リン酸）の加水分解反応でエネルギーが生じる。

また、図５中、ｆ_ａｐｐは、Ｐ_ＸＢからＸＢ_ＰｒｅＲへの遷移によって、状態分布がボルツマン分布則に基づく平衡状態となるように定められる。ｇ_ａｐｐについても、ＸＢ_ＰｒｅＲからＰ_ＸＢへの遷移によって、状態分布がボルツマン分布則に基づく平衡状態となるように定められる。

また、ＸＢ_ＰｒｅＲからＸＢ_{ＰｏｓｔＲ}への遷移、及び、ＸＢ_{ＰｏｓｔＲ}からＸＢ_ＰｒｅＲへの遷移において、ミオシンヘッド２２が首振り運動を行う。図５中、ｈ_ｆ、ｈ_ｂは、ミオシンヘッド２２の首振り運動によるミオシンアームの長さの変化によって、状態分布がミオシンヘッド内の化学エネルギーとミオシンアームの弾性エネルギーの和から定まるボルツマン分布則に基づく平衡状態となるように定められる。また、図５中、ｇ_ｘｂは、前記協調性の効果以外のＸＢ_{ＰｏｓｔＲ}からの解離のしやすさを表す遷移率である。また、図５中、αは、強結合状態ＸＢ_ＰreＲおよびＸＢ_{ＰｏｓｔＲ}において、アクチンからのミオシンヘッド２２の外れにくさを示す１以下の係数である。

ここで、サルコメア長（ＳＬ）に応じてアクチンフィラメント２１のオーバラップ状態が変化することについて説明する。まず、ＳＬは、次のような方法によって計算できる。すなわち、下記の式（１）を用いて、時刻Ｔにおける心筋モデル１３ａの各要素１３ａ＿１の線維方向の伸び（線維方向ストレッチ）λ（Ｔ）を計算する。

ただし、Ｆ（Ｔ）は、時刻Ｔでの心筋モデル１３ａの変形勾配テンソルである。

続いて、下記の式（２）を用いて、各要素１３ａ＿１の線維方向の速度（ストレッチ速度）

を計算する。

そして、下記の式（３）を用いて、時刻Ｔにおけるサルコメアモデル１３ｂのサルコメア長ＳＬ（Ｔ）を計算する。

ただし、ＳＬ_０は、時刻Ｔ＝０のときの無負荷状態におけるサルコメア長である。

そして、下記の式（４）を用いて、サルコメアモデル１３ｂのすべり速度

を計算する。

図６を参照して、サルコメア長（ＳＬ）に応じてアクチンフィラメント２１のオーバラップ状態が変化することについて説明する。図６は、ＳＬに応じてアクチンフィラメント２１のオーバラップ状態が変化することを説明するための図である。図６の例では、時刻Ｔ１、Ｔ２、Ｔ３のそれぞれの場合におけるサルコメア長ＳＬ（Ｔ１）、ＳＬ（Ｔ２）、ＳＬ（Ｔ３）が示されている。図６の上の例に示すように、サルコメア長（ＳＬ）が長すぎると中央部のミオシンヘッドに対して結合できるアクチンがなくなる。また、図６下の例に示すように、サルコメア長（ＳＬ）が短すぎた場合も、アクチンフィラメント２１がオーバラップする部分が中央付近にでき、この部分に位置するミオシンヘッドがアクチンに結合しにくくなる。そのため、本実施例では、アクチンフィラメントのオーバラップの状態をミオシンヘッド２２の状態遷移に反映するための関数χ_ＬＡ及びχ_ＲＢを構成し、ミオシンヘッド２２の結合状態への状態遷移において、関数χ_ＬＡ及びχ_ＲＢを参照して遷移率を決定することもできる。

本実施例では、従来の常微分方程式が用いられたモデルのように平均的な振る舞いを一つの代表ユニットで記述するのではなく、個々のＴ／Ｔユニット２０および個々のミオシンヘッド２２の状態遷移が近傍の状態も参照しながら制御される。それゆえ、本実施例によれば、ミクロのモデルの状態遷移をより現実に近いかたちで平均化などにともなう誤差を抑制してシミュレーションを行い、以下に述べるように心筋の連続体モデルとの正確な連成解析を行うことが可能となる。

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

制御部１４は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。図１に示すように、複数の心筋モデル処理部１４ａと、複数のサルコメアモデル処理部１４ｂと、表示制御部１４ｃと、受付部１４ｄとを有する。

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

サルコメアモデル処理部１４ｂは、第１の計算部１５ａ、第２の計算部１５ｂを有する。

ここで、制御部１４で実行される処理の一例について図７を参照して説明する。図７は、実施例に係る生体シミュレーション処理の手順を示すフローチャートである。生体シミュレーション処理は、例えば、入力部１１から制御部１４に、生体シミュレーション処理を実行する指示が入力された場合に実行される。なお、生体シミュレーション処理では、心筋モデル１３ａの有限要素解析が時間ΔＴ間隔で実行される。また、サルコメアモデルのモンテカルロシミュレーション処理では、サルコメアモデル１３ｂのモンテカルロ法を用いたシミュレーションが、時間区分［Ｔ，Ｔ＋ΔＴ］がｎ個のステップの細かい時間区分に分割された時間Δｔ（ΔＴ／ｎ）間隔で実行される。ΔＴの値としては、例えば、２．５ミリ秒を採用することができる。また、Δｔの値としては、例えば、５マイクロ秒を採用することができる。モンテカルロ法では、遷移率と時間区分Δｔの積が１を超えないことが要請されるので、上記のような細かい時間刻みとなる。一方で、有限要素解析も細かい時間刻みで計算した方が精度の上から好ましいが、陰解法における線形求解に時間がかかるため上記のような時間区分ΔＴが適当である。

図７に示すように、心筋モデル処理部１４ａ及びサルコメアモデル処理部１４ｂは、シミュレーション処理を実行する（Ｓ１０１）。図８は、実施例に係るシミュレーション処理の手順を示すフローチャートである。図８に示すように、第１の計算部１５ａは、時間幅［Ｔ，Ｔ＋ΔＴ］における初期化を行う（Ｓ２０１）。例えば、Ｓ２０１で、第１の計算部１５ａは、変数ＳｕｍＬ_ＩＲに「０」を設定する。また、Ｓ２０１で、第１の計算部１５ａは、変数Ｘ_ＢＤに「０」を設定する。また、Ｓ２０１で、第１の計算部１５ａは、変数Ｘ_ＢＤ０に「０」を設定する。また、Ｓ２０１で、第１の計算部１５ａは、フラグｆｌａｇ_Ｄ０に「１」を設定する。また、Ｓ２０１で、第１の計算部１５ａは、変数Ｌ_Ｓ０に「Ｌ_ｓ（Ｔ）」を設定する。ここで、Ｌ_ｓ（Ｔ）は、時刻Ｔにおけるミオシンアーム２４の伸びのうちフィラメント間のすべり速度による寄与分のみを抜き出したものである。この変数は、処理開始時に「０」に設定され、後述するように毎回の有限要素ステップの終了後に正しい値に更新される。

図９を参照して、ミオシンアーム２４の伸びについての計算方法の一例について説明する。図９は、ミオシンアームの伸びについての計算方法の一例について説明するための図である。図９の例は、時刻ｔｂでミオシンヘッド２２がアクチンフィラメント２１に結合し、時刻Ｔに至るまでの間、結合が維持されている場合を示す。また、図９の例は、ミオシンヘッド２２が時刻ｔｂから時刻Ｔまでの間で首振り運動をした場合を示す。図９の例に示す場合において、第１の計算部１５ａは、下記の式（５）を用いて、Ｌ_ｓ（Ｔ）を計算する。

ただし、Ｌ_ＩＮＩＴは、結合時のミオシンアーム２４の初期伸びである。また、Ｌ_ＲＯＴは、結合後のミオシンヘッド２２の首振り運動による伸びである。

図８に戻り、第１の計算部１５ａは、変数ｋの値を「１」に設定する（Ｓ２０２）。続いて、第１の計算部１５ａ及び第２の計算部１５ｂは、モンテカルロシミュレーション処理を実行する（Ｓ２０３）。なお、Ｓ２０３で実行されるモンテカルロシミュレーション処理は、時刻（Ｔ＋ｋΔｔ）における処理である。図１０は、実施例に係るモンテカルロシミュレーション処理の手順を示すフローチャートである。図１０に示すように、第１の計算部１５ａは、ミオシンヘッド２２の状態を判定し、結合状態であるか否かを判定する（Ｓ３０１）。例えば、Ｓ３０１では、第１の計算部１５ａは、Ｐ_ＸＢ、ＸＢ_ＰｒｅＲ、ＸＢ_{ＰｏｓｔＲ}の何れかである場合には、結合状態であると判定し、Ｎ_ＸＢである場合には、結合状態でないと判定する。

結合状態でない場合（Ｓ３０１；Ｎｏ）には、第１の計算部１５ａは、Ｓ３０３へ進む。一方、結合状態である場合（Ｓ３０１；Ｙｅｓ）には、第１の計算部１５ａは、各種の前処理を行う（Ｓ３０２）。例えば、Ｓ３０２では、第１の計算部１５ａは、変数Ｘ_ＢＤの値を１つインクリメントする。また、Ｓ３０２では、第１の計算部１５ａは、フラグｆｌａｇ_Ｄ０の値が「１」である場合には、変数Ｘ_ＢＤ０の値を１つインクリメントする。そして、Ｓ３０２では、第１の計算部１５ａは、下記の式（６）を用いて、ミオシンアームの伸びＬを計算する。

ただし、Ｌ_ＩＲは、上述したＬ_ＩＮＩＴとＬ_ＲＯＴとの和を示す。ここで、Ｌ_ＩＮＩＴは結合状態への遷移直後のミオシンアームの伸びであり、例えば熱揺らぎを考慮してミオシンアームの弾性エネルギーから定まるボルツマン分布をもとに計算される。

そして、第１の計算部１５ａは、乱数を生成し、生成した乱数を用いて、ミオシンヘッド２２の状態を計算する（Ｓ３０３）。これにより、第１の計算部１５ａは、ミオシンヘッド２２の確率的な挙動を計算することができる。

続いて、第２の計算部１５ｂは、ミオシンヘッド２２の状態が遷移したか否かを判定する（Ｓ３０４）。ミオシンヘッド２２の状態が遷移しない場合（Ｓ３０４：Ｎｏ）には、第２の計算部１５ｂは、Ｓ３０６へ進む。

一方、ミオシンヘッド２２の状態が遷移した場合（Ｓ３０４：Ｙｅｓ）には、第２の計算部１５ｂは、状態遷移処理を行う（Ｓ３０５）。例えば、Ｓ３０５では、第２の計算部１５ｂは、ミオシンヘッド２２の状態が、非結合状態Ｎ_ＸＢから結合状態Ｐ_ＸＢへ遷移した場合には、変数Ｌ_ＩＮＩＴの値を変数Ｌ_ＩＲに設定し、フラグｆｌａｇ_Ｄ０の値を「０」に設定する。また、Ｓ３０５では、第２の計算部１５ｂは、ミオシンヘッド２２の状態が、結合状態Ｐ_ＸＢ、ＸＢ_ＰｒｅＲ、ＸＢ_{ＰｏｓｔＲ}のいずれかから、非結合状態Ｎ_ＸＢへ遷移した場合には、変数Ｌ_ＩＲ、変数Ｌ_Ｓ０及び変数Ｘ_Ｂｄのそれぞれの値を「０」に設定する。また、Ｓ３０５では、第２の計算部１５ｂは、ミオシンヘッド２２の状態が、結合状態Ｐ_ＸＢ、ＸＢ_ＰｒｅＲ、ＸＢ_{ＰｏｓｔＲ}のいずれかの状態のままで、ミオシンヘッド２２が首振り運動を行い、回転した場合には、次の処理を行う。すなわち、第２の計算部１５ｂは、変数Ｌ_ＩＲの値と、ミオシンヘッド２２の回転により伸びたミオシンアーム２４の長さＸ_{ＳＷＩＮＧ}との和（Ｌ_ＩＲ＋Ｘ_{ＳＷＩＮＧ}）を、変数Ｌ_ＩＲに設定する。

そして、第２の計算部１５ｂは、後処理を実行し（Ｓ３０６）、処理結果を内部メモリに格納し、リターンする。例えば、Ｓ３０６では、第２の計算部１５ｂは、変数ＳｕｍＬ_ＩＲの値と、変数Ｌ_ＩＲの値との和（ＳｕｍＬ_ＩＲ＋Ｌ_ＩＲ）を変数ＳｕｍＬ_ＩＲに設定する。また、Ｓ３０６では、第２の計算部１５ｂは、変数ＳｕｍＸＢ_Ｄの値と、変数ＸＢ_Ｄの値との和（ＳｕｍＸＢ_Ｄ＋ＸＢ_Ｄ）を変数ＳｕｍＸＢ_Ｄに設定する。上述した処理が有限要素解析における１ステップの間にｎステップ回行われることで、変数ＳｕｍＸＢ_Ｄには、有限要素解析における１ステップの間における、ミオシンヘッド２２の結合状態の連続回数が設定される。また、変数ＳｕｍＬ_ＩＲには、有限要素解析における１ステップの間におけるＬ_ＩＲの和が設定され、後述するようにこれらの和は有限要素解析におけるアクティブ応力を計算する際に用いられる。

図８の説明に戻り、第２の計算部１５ｂは、変数ｋの値が、ｎを超えたか否かを判定する（Ｓ２０４）。変数ｋの値がｎを超えていない場合（Ｓ２０４；Ｎｏ）には、第２の計算部１５ｂは、変数ｋの値を１つインクリメントし（Ｓ２０５）、Ｓ２０３に戻る。一方、変数ｋの値がｎを超えた場合（Ｓ２０４；Ｙｅｓ）には、心筋モデル処理部１４ａは、有限要素解析を実行する（Ｓ２０６）。図１１は、実施例に係る有限要素解析の手順を示すフローチャートである。図１１に示すように、心筋モデル処理部１４ａは、１ミオシンアームの平均諸量を計算する（Ｓ４０１）。例えば、Ｓ４０１では、心筋モデル処理部１４ａは、下記の式（７）を用いて、有限要素解析の１ステップの間における各ミオシンアーム２４の伸びの平均値Ｌ_ＡＶＲを計算する。

また、Ｓ４０１では、心筋モデル処理部１４ａは、下記の式（８）を用いて、有限要素解析の１ステップの間におけるミオシンアーム２４の力の平均値Ｆ_ＭＨ、及び、剛性の平均値Ｋ_ＭＨを計算する。

ここで、Ｗ_ａｒｍは、ミオシンアーム２４の伸びＬに対する関数として与えられたバネの弾性エネルギーである。

また、下記の式（９）は、与えられた線維方向のストレッチの変分δλに対するミオシンアームの仮想的な仕事δＷ_ＭＨを表したものである。

そして、心筋モデル処理部１４ａは、下記の式（１０）に基づき、心筋の連続体モデル上での任意の歪みの変分δEに対して、有限要素解析の１ステップにおける単位体積あたりの心筋モデル１３ａがした仮想的な仕事とミオシンアームの集合がした仮想的な仕事が等しくなるようにアクティブ応力Ｓ_{ａｃｔｉｖｅ}を計算する。以下、その方法の一例について説明する。

式（１０）において、ｉｍは、ミオシンフィラメント２３上のミオシンヘッド２２のインデックスであり、ｉｓは、要素１３ａ＿１内のサルコメアモデル１３ｂのインデックスである。また、Ｒ_ＳＡは、筋肉に対するサルコメアが占める割合であり、ＳＡ_０は、１サルコメアモデル１３ｂあたりの断面積である。また、Ｓ_{ａｃｔｉｖｅ}は、連続体モデルのアクティブ応力テンソル、δＥはＧｒｅｅｎ−Ｌａｇｒａｎｇｅ歪みテンソルの変分であり、この二つのテンソルの積が連続体の仮想仕事を表す。

心筋モデル処理部１４ａは、式（１０）に基づき下記の式（１１）を用いて、単位断面積あたりの心筋モデル１３ａの収縮力Ｆを計算する（Ｓ４０２）。

また、Ｓ４０２では、心筋モデル処理部１４ａは、下記の式（１２）を用いて、単位断面積あたりの心筋モデル１３の剛性Ｋを計算する。

ここで、下記の式（１３）と（１４）はＧｒｅｅｎ−Ｌａｇｒａｎｇｅ歪みテンソルの変分と線維方向のストレッチλの変分の関係を示すものである。式（１３）より、式（１１）で計算したＦを用いてアクティブ応力テンソルＳ_{ａｃｔｉｖｅ}を下記の式（１５）に示すように計算すれば式（１０）が成立することがわかる。

それゆえ、心筋モデル処理部１４ａは、式（１５）を用いて、Ｓ_{ａｃｔｉｖｅ}を計算する（Ｓ４０３）。そして、有限要素解析をＮｅｗｍａｒｋ−β法のような陰解法で実行する場合には、心筋モデル処理部１４ａは、陰解法で用いられる剛性行列を、下記の式（１６）及び式（１７）を用いて計算する（Ｓ４０４）。

これらの式は、線維方向ストレッチλおよびその時間微分とＧｒｅｅｎ−Ｌａｇｒａｎｇｅ歪みテンソルＥおよびその時間微分の関係式（１３）および式（１４）から導かれたものである。Ｓ４０４では、心筋モデル処理部１４ａは、式（１０）により算出したアクティブ応力Ｓ_{ａｃｔｉｖｅ}から等価節点力を、式（１６）および式（１７）をもとに剛性行列を計算し、これらをすべての要素にわたって重ね合わせることによりＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復における線形方程式の右辺ベクトルと係数行列を作成する。残差ベクトルに対応する右辺ベクトルが十分小さければ（Ｓ４０５肯定）、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復を終了して処理結果を内部メモリに格納しリターンする。そうでなければ（Ｓ４０５否定）、心筋モデル処理部１４ａは、線形方程式を解き、得られた解をもとに連続体モデルの変位ベクトル、速度ベクトル、加速度ベクトルを更新し（Ｓ４０６）、Ｓ４０１へ戻り、次のＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復に移る。ただし、これらベクトルの時間ステップＴ＋ΔＴにおける初期値は、これらベクトルの時刻Ｔにおける値とする。式（７）にみるように各ミオシンアームの時間幅［Ｔ，Ｔ＋ΔＴ］における平均伸びＬ_ＡＶＲを時刻Ｔ＋ΔＴにおける線維方向ストレッチλの時間微分から計算しているので、収縮力の計算において連続体モデルの運動のフィードバックが強連的にかかり、剛性行列に対しても時間幅［Ｔ，Ｔ＋ΔＴ］での架橋運動の影響が適切に反映される結果となり、安定した計算が可能となる。

Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復が収束したら有限要素解析を終了し、図８に戻り、心筋モデル処理部１４ａは、下記の式（１８）を用いて、変数Ｌｓを次の時間ステップでの計算のために更新する（Ｓ２０７）。この更新によりＬ_Ｓには時刻Ｔ＋ΔＴまでにバネの伸びＬの中ですべり速度が寄与した分が正しく設定し直される。

受付部１４ｄは、サルコメアモデル１３ｂに定義されるＴ／Ｔユニット２０の状態と遷移、及び、ミオシンヘッド２２の状態と遷移の定義を受け付ける。そして、受付部１４ｄは、受け付けた内容を各種のモデルに反映させる。図１２〜図１８は、実施例に係る受付部が受け付ける情報を説明するための図である。図１２の例に示すように、受付部１４ｄは、ユーザ、例えば、心臓を医学的に研究している研究者などによって、入力部１１から、Ｔ／Ｔユニット２０の状態と遷移を受け付ける画面を表示させる指示が入力されると、次の処理を実行する。すなわち、受付部１４ｄは、Ｔ／Ｔユニット２０の状態と遷移を受け付ける画面３０を表示部１２に表示させる。図１２の例は、画面３０を介して、「Ｃａ−ｏｆｆ」及び「Ｃａ−ｏｎ」の２つの状態が定義された場合を示す。また、「Ｉｎｉｔｉａｌ Ｓｔａｔｅ」の横のチェックボックスにチェックを入れると、チェックが入れられた状態が、生体シミュレーション処理開始時の状態となる。図１２の例は、生体シミュレーション処理開始時の状態が、「Ｃａ−ｏｆｆ」であることを示す。

また、図１３の例に示すように、受付部１４ｄは、ユーザによって、入力部１１から、ミオシンヘッド２２の状態と遷移を受け付ける画面を表示させる指示が入力されると、次の処理を実行する。すなわち、受付部１４ｄは、ミオシンヘッド２２の状態と遷移を受け付ける画面３１を表示部１２に表示させる。図１３の例は、画面３１を介して、「Ｎ＿ＸＢ」、「Ｐ＿ＸＢ」、「ＸＢ＿ＰｒｅＲ」及び「ＸＢ＿ＰｏｓｔＲ」の４つの状態が定義された場合を示す。また、「Ｉｎｉｔｉａｌ Ｓｔａｔｅ」の横のチェックボックスにチェックを入れると、チェックが入れられた状態が、生体シミュレーション処理開始時の状態となる。図１３の例は、生体シミュレーション処理開始時の状態が、「Ｎ＿ＸＢ」であることを示す。また、図１３の例に示す画面３１は、４つの状態のそれぞれが「Ｎｏｂｉｎｄｉｎｇ（非結合状態）」または「Ｂｉｎｄｉｎｇ（結合状態）」のいずれであるかを設定するためのラジオボタンが設けられている。図１３の例は、「Ｎ＿ＸＢ」及び「Ｐ＿ＸＢ」が非結合状態であり、「ＸＢ＿ＰｒｅＲ」及び「ＸＢ＿ＰｏｓｔＲ」が結合状態であることを示す。

また、図１４及び図１５の例に示すように、受付部１４ｄは、ユーザによって、入力部１１から、関数を定義するための画面を表示させる指示が入力されると、次の処理を実行する。すなわち、受付部１４ｄは、関数を定義するための画面３２や画面３３を表示部１２に表示させる。図１４の例は、画面３２を介して、状態が「Ｃａ−ｏｆｆ」である場合には、「Ｋ＿ＮＰ０」を返し、状態が「Ｃａ−ｏｎ」である場合には、「Ｋ＿ＮＰ１」を返すようなＴ／Ｔユニット２０の状態関数ｋｎｐが定義された場合を示す。また、図１５の例は、画面３３を介して、状態が「Ｎ＿ＸＢ」である場合には、「０」を返し、状態が「Ｐ＿ＸＢ」、「ＸＢ＿ＰｒｅＲ」及び「ＸＢ＿ＰｏｓｔＲ」である場合には、「１」を返すようなミオシンヘッド２２の状態関数ｎｇが定義された場合を示す。なお、定義された状態関数は、遷移率の定義において参照される。

また、図１６の例に示すように、受付部１４ｄは、ユーザが入力部１１を操作して、「Ｎ＿ＸＢ」から「Ｐ＿ＸＢ」へ向かう矢印が選択されると、「Ｎ＿ＸＢ」から「Ｐ＿ＸＢ」への遷移率を定義するための画面３４を表示部１２に表示させる。図１６の例は、「Ｎ＿ＸＢ」から「Ｐ＿ＸＢ」への遷移率を定義する場合に、ｇｅｔ（ＴＴ，ｋｎｐ）によりミオシンヘッド２２の直上のアクチン切片を配下とするＴ／Ｔユニット２０の状態に応じた値を返す関数ｋｎｐを用いる場合を示す。また、図１６の例は、ｇｅｔ（ＭＨ，ｎｇ，−１）及びｇｅｔ（ＭＨ，ｎｇ，１）によりミオシンヘッド２２の左および右隣のミオシンヘッド２２上の関数値ｎｇを参照し、ＧＡＭＭＡとして１より大きなパラメータが適用された場合を示す。これにより、近傍のミオシンヘッド２２間の協調性が表現される。さらに、ｘｉ＿ｏｖｅｒｌａｐ（）はサルコメア長からきまるアクチンフィラメント２１のオーバラップ状態を反映した関数であり、先の図６の例に示す関数χ_ＲＡと関数χ_ＬＡとを掛け合わせた関数である。

また、図１７の例に示すように、受付部１４ｄは、ユーザが入力部１１を操作して、「ＸＢ＿ＰｒｅＲ」から「ＸＢ＿ＰｏｓｔＲ」へ向かう矢印が選択されると、次の処理を行う。すなわち、受付部１４ｄは、ミオシンヘッド２２の回転を伴う、「ＸＢ＿ＰｒｅＲ」から「ＸＢ＿ＰｏｓｔＲ」への遷移率を定義するための画面３５を表示部１２に表示させる。図１７の例では、遷移がミオシンヘッド２２の回転を伴う遷移であれば、「Ｍｙｏｓｉｎ ｈｅａｄ ｓｗｉｎｇ」の横に存在するチェックボックスがユーザによりチェックされる。また、図１７の例では、回転にともなうミオシンアーム２４の伸びの増分が、Ｘ＿ＳＷＩＮＧとしてユーザにより入力される。遷移率の定義に際しては、ミオシンヘッド２２の元の状態「ＸＢ＿ＰｒｅＲ」のミオシンアーム２４の伸びＬが関数ａｒｍ＿ｌｅｎｇｔｈ（）で呼び出される。そして、遷移後の状態におけるミオシンアームの伸びａｒｍ＿ｌｅｎｇｔｈ（）＋Ｘ＿ＳＷＩＮＧでのミオシンアーム２４の弾性エネルギーがｓｐｒｉｎｇ＿ｅｎｅｒｇｙを呼び出されることにより求められる。そして、遷移後のミオシンヘッド２２の内部の化学エネルギーＥ＿ＸＢＰｏｓｔＲと足し合わせたトータルエネルギーから定まるボルツマン因子より遷移率ｒが定められる。

また、図１８の例に示すように、受付部１４ｄは、ユーザが入力部１１を操作して、「ＸＢ＿ＰｏｓｔＲ」から「Ｎ＿ＸＢ」へ向かう矢印が選択されると、次の処理を行う。すなわち、受付部１４ｄは、ミオシンヘッド２２の解離を伴う、「ＸＢ＿ＰｏｓｔＲ」から「Ｎ＿ＸＢ」への遷移率を定義するための画面３６を表示部１２に表示させる。図１８は、結合状態から非結合状態への遷移に対しては、ミオシンヘッド２２からＡＤＰが解離し、新たにＡＴＰが結合するものと仮定して、「ＡＴＰ ｂｉｎｄｉｎｇ」と、「ＡＤＰ ｒｅｌｅａｓｅ」のチェックボックスがユーザによりチェックされた場合を示す。本実施例に係る生体シミュレーション装置１０は、上述したようなユーザによる指定をもとして生体シミュレーション実行時に分子の消費量および生成量を計算することができる。

また、受付部１４ｄは、上述の入力で定義されたＴ／Ｔユニットとミオシンヘッドの状態と状態間の遷移情報をもとにモンテカルロシミュレーションを実行するコードを自動生成する。図１９は、ミオシンヘッドの状態と状態間遷移の定義をもとにモンテカルロステップを実行するコード（モンテカルロコード）を自動生成した結果を示す図である。図１９に示すように、かかるコードは、最初に[０，１]間の乱数ｒを生成し、現状態が何であるかに応じて処理を実行する。もし、現状態が結合状態であれば処理ｕｐｄａｔｅ＿ＸＢ（）により、先のステップＳ３０２で説明したように、ミオシンヘッド２２の変数（ＸＢ_Ｄ、Ｌ）を更新する。次に、コードは、Ｎ個ある状態のうち現状態に応じて、ｐｒｏｃｅｓｓ＿ｔｒａｎｓｉｔｉｏｎ＿ｓｔａｔｅ１、・・・、ｐｒｏｃｅｓｓ＿ｔｒａｎｓｉｔｉｏｎ＿ｓｔａｔｅＮのひとつを実行する。この処理では、乱数ｒの値に応じてｕｐｄａｔｅ＿ｓｔａｔｅ（）により状態を更新し、ｕｐｄａｔｅ＿ｖａｒｉａｂｌｅｓ（）により先のステップＳ３０５で説明したように、遷移の種類に応じた変数の更新を行う。

図１９に示すコードと、心筋モデル１３ａの有限要素コードとを、図８に示す処理のように組み合わせ、心筋モデル１３ａとの連成解析を通して心拍動のシミュレーションを実行するときの例を図２０Ａ及び図２０Ｂに示す。ここでは、左心室のモデルとして図２０Ｂに示す断面を回転して得られる回転対称な連続体を採用している。図２０Ａは、受付部１４ｄにより表示部１２に表示された画面を介して、「＃ ｏｆ ｅｌｅｍｅｎｔｓ ｉｎ Ｒ」および「＃ ｏｆ ｅｌｅｍｅｎｔｓ ｉｎ Ｌ」で貫壁方向および長手方向の要素数がユーザにより指定された場合を示す。また、図２０Ａは、「＃ ｏｆ ｈｅａｒｔ ｂｅａｔｓ」で、シミュレートする拍動回数がユーザにより指定された場合を示す。また、図２０Ａは、「α＿ｂａｓｅ」および「α＿ａｐｅｘ」で、図２０Ｂの例に示すように、線維方向の分布がユーザにより指定された場合を示す。

図２１Ａ〜図２１Ｃは、シミュレーション結果の一例を示す図である。図２１Ａの例に示すシミュレーション結果は、各拍動での拍出量(Ｅｊｅｃｔｉｏｎ)、拍出量とエネルギーの比（Ｅｊｅｃｔｉｏｎ Ｆｒａｃｔｉｏｎ）、心筋モデルの成した仕事量(Ｍｕｓｃｌｅ Ｗｏｒｋ)を出力したものである。さらに、図２１Ａの例に示すシミュレーション結果では、血液拍出で成した仕事量(Ｅｊｅｃｔｉｏｎ ｗｏｒｋ)、ＡＴＰ消費量(ＡＴＰ ｃｏｎｓｕｍｐｔｉｏｎ)、効率(Ｅｆｆｅｃｉｅｎｃｙ)が出力されている。また、図２１Ｂの例に示すシミュレーション結果は、心室内容積と心室内圧力の時間変化および血液拍出に要する仕事率とＡＴＰ消費率を加水分解エネルギーで換算したものの時間変化を出力したものである。また、図２１Ｃの例に示すシミュレーション結果は、１拍動において、各要素が成した仕事の積算値の時間変化を示したものである。図２１Ａ〜図２１Ｃに示すシミュレーション結果から、本実施例に係る生体シミュレーション装置１０によれば、サルコメアモデル１３ｂにおける分子レベルの現象と心臓の拍動性能の関係を解析することが可能となる。

上述してきたように、本実施例に係る生体シミュレーション装置１０は、生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、複数の状態間の遷移率が定義されたサルコメアモデル１３ｂを用いて、次の処理を行う。すなわち、生体シミュレーション装置１０によれば、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態を計算する。そして、生体シミュレーション装置１０は、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態に基づいて、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動を計算する。それゆえ、生体シミュレーション装置１０によれば、複数のアクチン及びミオシンのそれぞれの確率的な挙動を計算することができる。したがって、生体シミュレーション装置１０によれば、精緻なシミュレーション結果を得ることができる。

また、本実施例に係る生体シミュレーション装置１０は、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動に基づいて、サルコメアの挙動を計算し、サルコメアの挙動に基づいて、筋肉の挙動を計算する。それゆえ、本実施例に係る生体シミュレーション装置１０によれば、筋肉の挙動を計算するシミュレーションと、サルコメアの挙動を計算するシミュレーションとを連成させることができる。

また、本実施例に係る生体シミュレーション装置１０は、Δｔ間隔で、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動を計算する。そして、生体シミュレーション装置１０は、Δｔよりも長いΔＴ間隔で、ΔＴの間に計算された複数のサルコメアの挙動に基づいて、筋肉の挙動を計算する。このように、本実施例に係る生体シミュレーション装置１０は、ΔＴの間に複数回計算された複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動を用いて、筋肉の挙動を計算する。そのため、生体シミュレーション装置１０によれば、筋肉の挙動を計算する時間刻みΔＴを有限要素解析の都合に合わせて選ぶことができ、心臓の拍動のような秒オーダーの現象を適度な時間で実行することが可能となる。このような時間間隔のギャップがあるにも拘わらす、式（１０）に示したようにサルコメアモデルと心筋モデルで仕事量を一致させているので精度の良い連成解析が可能となる。さらに、式（７）で示しているように、時刻Ｔ＋ΔＴでのストレッチ速度から時間幅［Ｔ，Ｔ＋ΔＴ］でのミオシンアームの平均伸びＬ_ＡＶＲを算出しているので、サルコメアモデルと心筋モデルの仕事量が強く連成され適度な大きさの時間ステップΔＴでの安定した解析が可能となる。したがって、本実施例に係るシミュレーション装置１０によれば、精緻なシミュレーション結果を適切な計算時間内で得ることができる。

また、本実施例に係る生体シミュレーション装置１０が扱うサルコメアモデル１３ｂは、アクチンにミオシンが結合している場合には、当該ミオシンの所定範囲内に存在する他のミオシンが、アクチンに結合する状態となる確率が高くなるように定義されている。さらに、直上にあるアクチンフィラメントのオーバーラップ状態に応じて、個々のミオシンヘッドの状態遷移を定義している。すなわち、現実のサルコメアに類似する挙動を示すようにサルコメアモデル１３ｂが定義されている。そして、シミュレーション装置１０は、かかるサルコメアモデル１３ｂを用いて、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態を計算する。それゆえ、シミュレーション装置１０は、現実のサルコメアと類似するような状態を計算することができる。

また、生体シミュレーション装置１０は、ミオシンがアクチンに結合した際にミオシンアームが伸びた長さ、ミオシンヘッドの回転によってミオシンアームが伸びた長さ、ミオシンのすべり速度の積分値に基づいて、ミオシンが伸びた長さを計算する。

また、本実施例に係る生体シミュレーション装置１０は、サルコメアモデル１３ｂに定義される複数のアクチンの複数の状態及びアクチンの状態間の遷移、並びに、複数のミオシンの複数の状態及びミオシンの状態間の遷移を受け付ける。そして、生体シミュレーション装置１０は、受け付けた内容をもとに、モンテカルロステップを実行するモンテカルロコードを生成する。それゆえ、例えば、医学的に心臓の研究をしているがプログラミングが苦手な研究者などのユーザであっても、プログラミングを行うことなく、任意のサルコメアモデルのシミュレーションを行うことができる。

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

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

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

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

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

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

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

ＲＯＭ３２０には、ＯＳなどの基本プログラムが記憶されている。また、ＨＤＤ３３０には、上記の実施例１で示す心筋モデル処理部１４ａ、サルコメアモデル処理部１４ｂ、第１の計算部１５ａ、第２の計算部１５ｂと同様の機能を発揮する生体シミュレーションプログラム３３０ａが予め記憶される。なお、生体シミュレーションプログラム３３０ａについては、適宜分離しても良い。

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

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

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

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

１０ 生体シミュレーション装置
１３ａ 心筋モデル
１３ｂ サルコメアモデル
１４ａ 心筋モデル処理部
１４ｂ サルコメアモデル処理部
１５ａ 第１の計算部
１５ｂ 第２の計算部

Claims (6)

1. コンピュータに、
   生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、該複数の状態間の遷移率が定義されたモデルを用いて、それぞれの前記アクチンについて、前記アクチンに前記ミオシンが結合しているか否かに基づいて前記アクチンの状態を算出し、それぞれの前記ミオシンについて、前記アクチンとの結合状態に基づいて前記ミオシンの状態を算出し、
   算出したそれぞれの前記アクチンおよび前記ミオシンの各状態に基づいて前記ミオシンが伸びた長さを算出し、前記ミオシンが伸びた長さに基づいて、それぞれの前記アクチンおよび前記ミオシンの各状態間の状態遷移を算出し、
   算出したそれぞれの前記アクチンおよび前記ミオシンの各状態間の状態遷移に基づいて、前記サルコメアの仕事を算出し、
   算出した前記サルコメアの仕事に基づいて前記筋肉の収縮力および剛性を算出し、算出した前記筋肉の収縮力および剛性に基づいて、前記筋肉の変位ベクトル、速度ベクトルおよび加速度ベクトルを算出する、処理を実行させ、
   前記状態遷移を算出する処理は、前記ミオシンが前記アクチンに結合した際に該ミオシンに含まれるミオシンアームが伸びた長さ、該ミオシンが該アクチンに結合してからの該ミオシンに含まれるミオシンヘッドの回転によって該ミオシンアームが伸びた長さ、及び、前記ミオシンが前記アクチンに結合してから現在までの前記ミオシンのすべり速度の積分値に基づいて、該ミオシンが伸びた長さを算出する、
   ことを特徴とする生体シミュレーションプログラム。
2. 前記状態遷移を算出する処理は、第１の時間間隔で、前記それぞれの前記アクチンおよび前記ミオシンの各状態間の状態遷移を算出し、
   前記変位ベクトル、前記速度ベクトルおよび前記加速度ベクトルを算出する処理は、前記第１の時間間隔よりも長い第２の時間間隔で、該第２の時間の間に算出された複数の前記サルコメアの仕事に基づいて、前記筋肉の収縮力および剛性を算出し、算出した前記筋肉の収縮力および剛性に基づいて、前記筋肉の変位ベクトル、前記速度ベクトルおよび前記加速度ベクトルを算出する、
   ことを特徴とする請求項１に記載の生体シミュレーションプログラム。
3. 前記アクチンの状態及び前記ミオシンの状態を算出する処理は、アクチンにミオシンが結合している場合には、該ミオシンの所定範囲内に存在する他のミオシンが、アクチンに結合している状態となる確率が高くなるように定義された前記モデルを用いて、前記アクチンの状態及び前記ミオシンの状態を算出する、
   ことを特徴とする請求項１または２に記載の生体シミュレーションプログラム。
4. 前記コンピュータに、
   前記モデルに定義される前記複数のアクチンの複数の状態及び該アクチンの状態間の遷移、並びに、前記複数のミオシンの複数の状態及び該ミオシンの状態間の遷移を受け付ける、
   処理をさらに実行させることを特徴とする請求項１〜３のいずれか１項に記載の生体シミュレーションプログラム。
5. コンピュータが、
   生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、該複数の状態間の遷移率が定義されたモデルを用いて、それぞれの前記アクチンについて、前記アクチンに前記ミオシンが結合しているか否かに基づいて前記アクチンの状態を算出し、それぞれの前記ミオシンについて、前記アクチンとの結合状態に基づいて前記ミオシンの状態を算出し、
   算出したそれぞれの前記アクチンおよび前記ミオシンの各状態に基づいて前記ミオシンが伸びた長さを算出し、前記ミオシンが伸びた長さに基づいて、それぞれの前記アクチンおよび前記ミオシンの各状態間の状態遷移を算出し、
   算出したそれぞれの前記アクチンおよび前記ミオシンの各状態間の状態遷移に基づいて、前記サルコメアの仕事を算出し、
   算出した前記サルコメアの仕事に基づいて前記筋肉の収縮力および剛性を算出し、算出した前記筋肉の収縮力および剛性に基づいて、前記筋肉の変位ベクトル、速度ベクトルおよび加速度ベクトルを算出する、処理を実行し、
   前記状態遷移を算出する処理は、前記ミオシンが前記アクチンに結合した際に該ミオシンに含まれるミオシンアームが伸びた長さ、該ミオシンが該アクチンに結合してからの該ミオシンに含まれるミオシンヘッドの回転によって該ミオシンアームが伸びた長さ、及び、前記ミオシンが前記アクチンに結合してから現在までの前記ミオシンのすべり速度の積分値に基づいて、該ミオシンが伸びた長さを算出する、
   ことを特徴とする生体シミュレーション方法。
6. 生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、該複数の状態間の遷移率が定義されたモデルを用いて、それぞれの前記アクチンについて、前記アクチンに前記ミオシンが結合しているか否かに基づいて前記アクチンの状態を算出し、それぞれの前記ミオシンについて、前記アクチンとの結合状態に基づいて前記ミオシンの状態を算出する第１の算出部と、
   算出したそれぞれの前記アクチンおよび前記ミオシンの各状態に基づいて前記ミオシンが伸びた長さを算出し、前記ミオシンが伸びた長さに基づいて、それぞれの前記アクチンおよび前記ミオシンの各状態間の状態遷移を算出する第２の算出部と、
   算出したそれぞれの前記アクチンおよび前記ミオシンの各状態間の状態遷移に基づいて、前記サルコメアの仕事を算出し、算出した前記サルコメアの仕事に基づいて前記筋肉の収縮力および剛性を算出し、算出した前記筋肉の収縮力および剛性に基づいて、前記筋肉の変位ベクトル、速度ベクトルおよび加速度ベクトルを算出する第３の算出部と、を有し、
   前記第２の算出部は、前記ミオシンが前記アクチンに結合した際に該ミオシンに含まれるミオシンアームが伸びた長さ、該ミオシンが該アクチンに結合してからの該ミオシンに含まれるミオシンヘッドの回転によって該ミオシンアームが伸びた長さ、及び、前記ミオシンが前記アクチンに結合してから現在までの前記ミオシンのすべり速度の積分値に基づいて、該ミオシンが伸びた長さを算出する、
   ことを特徴とする生体シミュレーション装置。

Images