(57)【要約】      （修正有）【課題】精緻なシミュレーションの結果を適切な計算時間内に得る。【解決手段】生体シミュレーション装置１０は、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態に基づいて、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動を計算する。前述の計算を、筋肉の連続体モデルの各要素に埋め込まれた複数のサルコメアモデルに対して各要素の変形を考慮したサルコメア長のもとで実行する。連続体の変形を求める有限要素解析の際の収縮力の計算においては、サルコメアモデルでの計算結果を参照して収縮力を求める。このサルコメアモデルと筋肉の連続体モデルの連成計算を、サルコメアモデルにおけるモンテカルロステップの細かい時間刻みと筋肉の連続体モデルにおける有限要素解析の粗い時間刻みの異なる時間刻みに対して、両者の間で仕事量が一致するように実行して、計算精度を確保する。【選択図】図１
2014149649.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータに、
生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、該複数の状態間の遷移率が定義されたモデルを用いて、前記複数のアクチンのそれぞれの状態及び前記複数のミオシンのそれぞれの状態を計算し、
前記複数のアクチンのそれぞれの状態及び前記複数のミオシンのそれぞれの状態に基づいて、前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動を計算し、
前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動に基づいて、前記サルコメアの挙動を計算し、
前記サルコメアの挙動に基づいて、前記筋肉の挙動を計算する
処理を実行させることを特徴とする生体シミュレーションプログラム。
【請求項２】
前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動を計算する処理は、第１の時間間隔で、前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動を計算し、
前記筋肉の挙動を計算する処理は、前記第１の時間よりも長い第２の時間間隔で、該第２の時間の間に計算された複数の前記サルコメアの挙動に基づいて、前記筋肉の挙動を計算する
ことを特徴とする請求項１記載の生体シミュレーションプログラム。
【請求項３】
前記複数のアクチンのそれぞれの状態及び前記複数のミオシンのそれぞれの状態を計算する処理は、アクチンにミオシンが結合している場合には、該ミオシンの所定範囲内に存在する他のミオシンが、アクチンに結合している状態となる確率が高くなるように定義された前記モデルを用いて、前記複数のアクチンのそれぞれの状態及び前記複数のミオシンのそれぞれの状態を計算する
ことを特徴とする請求項１または２記載の生体シミュレーションプログラム。
【請求項４】
前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動を計算する処理は、前記ミオシンが前記アクチンに結合した際に該ミオシンに含まれるミオシンアームが伸びた長さ、該ミオシンが該アクチンに結合してからの該ミオシンに含まれるミオシンヘッドの回転によって該ミオシンアームが伸びた長さ、及び、前記ミオシンが前記アクチンに結合してから現在までの前記ミオシンのすべり速度の積分値に基づいて、該ミオシンが伸びた長さを計算する
ことを特徴とする請求項１〜３のいずれか１項記載の生体シミュレーションプログラム。
【請求項５】
前記コンピュータに、
前記モデルに定義される前記複数のアクチンの複数の状態及び該アクチンの状態間の遷移、並びに、前記複数のミオシンの複数の状態及び該ミオシンの状態間の遷移を受け付ける
処理をさらに実行させることを特徴とする請求項１〜４のいずれか１項記載の生体シミュレーションプログラム。
【請求項６】
コンピュータが、
生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、該複数の状態間の遷移率が定義されたモデルを用いて、前記複数のアクチンのそれぞれの状態及び前記複数のミオシンのそれぞれの状態を計算し、
前記複数のアクチンのそれぞれの状態及び前記複数のミオシンのそれぞれの状態に基づいて、前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動を計算し、
前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動に基づいて、前記サルコメアの挙動を計算し、
前記サルコメアの挙動に基づいて、前記筋肉の挙動を計算する
処理を実行することを特徴とする生体シミュレーション方法。
【請求項７】
生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、該複数の状態間の遷移率が定義されたモデルを用いて、前記複数のアクチンのそれぞれの状態及び前記複数のミオシンのそれぞれの状態を計算する第１の計算部と、
前記複数のアクチンのそれぞれの状態及び前記複数のミオシンのそれぞれの状態に基づいて、前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動を計算する第２の計算部と、
前記複数のアクチンのそれぞれの挙動及び前記複数のミオシンのそれぞれの挙動に基づいて、前記サルコメアの挙動を計算し、前記サルコメアの挙動に基づいて、前記筋肉の挙動を計算する第３の計算部と
を有することを特徴とする生体シミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、生体シミュレーションプログラム、生体シミュレーション方法及び生体シミュレーション装置に関する。
【背景技術】
【０００２】
分子生物学の分野では、実験事実に基づきサルコメア内のミオシン・アクチン間の架橋運動を記述した多様なサルコメアモデルが提案されている。例えば、サルコメア内のミオシン・アクチン間の結合に寄与する分子に対して複数の代表的な状態が定義され、近傍分子間の相互作用および分子のエネルギーなどが考慮されて状態間の遷移率が定義されたサルコメアモデルがある。かかるサルコメア内の分子の運動と筋肉の連続体モデルを連成させた解析を行う際には、分子モデルに対してはその平均的な挙動を想定し状態濃度を変数とする常微分方程式を導く方法が一般的である。また、サルコメアモデルの計算結果から連続体モデルでの収縮力を算出し、かかる結果をもとに連続体としての筋肉の運動を計算する方法がとられることが一般的である。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】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ＲおよびＸＢＰｏｓｔＲにおいて、アクチンからのミオシンヘッド２２の外れにくさを示す１以下の係数である。
【００２３】
ここで、サルコメア長（ＳＬ）に応じてアクチンフィラメント２１のオーバラップ状態が変化することについて説明する。まず、ＳＬは、次のような方法によって計算できる。すなわち、下記の式（１）を用いて、時刻Ｔにおける心筋モデル１３ａの各要素１３ａ＿１の線維方向の伸び（線維方向ストレッチ）λ（Ｔ）を計算する。
【数１】
2014149649.tif 000003
ただし、Ｆ（Ｔ）は、時刻Ｔでの心筋モデル１３ａの変形勾配テンソルである。
【００２４】
続いて、下記の式（２）を用いて、各要素１３ａ＿１の線維方向の速度（ストレッチ速度）
【数２】
2014149649.tif 000004
を計算する。
【数３】
2014149649.tif 000005
【００２５】
そして、下記の式（３）を用いて、時刻Ｔにおけるサルコメアモデル１３ｂのサルコメア長ＳＬ（Ｔ）を計算する。
【数４】
2014149649.tif 000006
ただし、ＳＬ０は、時刻Ｔ＝０のときの無負荷状態におけるサルコメア長である。
【００２６】
そして、下記の式（４）を用いて、サルコメアモデル１３ｂのすべり速度
【数５】
2014149649.tif 000007
を計算する。
【数６】
2014149649.tif 000008
【００２７】
図６を参照して、サルコメア長（ＳＬ）に応じてアクチンフィラメント２１のオーバラップ状態が変化することについて説明する。図６は、ＳＬに応じてアクチンフィラメント２１のオーバラップ状態が変化することを説明するための図である。図６の例では、時刻Ｔ１、Ｔ２、Ｔ３のそれぞれの場合におけるサルコメア長ＳＬ（Ｔ１）、ＳＬ（Ｔ２）、ＳＬ（Ｔ３）が示されている。図６の上の例に示すように、サルコメア長（ＳＬ）が長すぎると中央部のミオシンヘッドに対して結合できるアクチンがなくなる。また、図６下の例に示すように、サルコメア長（ＳＬ）が短すぎた場合も、アクチンフィラメント２１がオーバラップする部分が中央付近にでき、この部分に位置するミオシンヘッドがアクチンに結合しにくくなる。そのため、本実施例では、アクチンフィラメントのオーバラップの状態をミオシンヘッド２２の状態遷移に反映するための関数χＬＡ及びχＲＢを構成し、ミオシンヘッド２２の結合状態への状態遷移において、関数χＬＡ及びχＲＢを参照して遷移率を決定することもできる。
【００２８】
本実施例では、従来の常微分方程式が用いられたモデルのように平均的な振る舞いを一つの代表ユニットで記述するのではなく、個々のＴ／Ｔユニット２０および個々のミオシンヘッド２２の状態遷移が近傍の状態も参照しながら制御される。それゆえ、本実施例によれば、ミクロのモデルの状態遷移をより現実に近いかたちで平均化などにともなう誤差を抑制してシミュレーションを行い、以下に述べるように心筋の連続体モデルとの正確な連成解析を行うことが可能となる。
【００２９】
記憶部１３は、例えば、フラッシュメモリなどの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。なお、記憶部１３は、上記の種類の記憶装置に限定されるものではなく、ＲＡＭ（Random  Access  Memory)、ＲＯＭ（Read  Only  Memory)等であってもよい。
【００３０】
制御部１４は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。図１に示すように、複数の心筋モデル処理部１４ａと、複数のサルコメアモデル処理部１４ｂと、表示制御部１４ｃと、受付部１４ｄとを有する。
【００３１】
複数の心筋モデル処理部１４ａのそれぞれは、心筋モデル１３ａ内の複数の要素のそれぞれに対応し、各心筋モデル処理部１４ａは、対応する心筋モデル１３ａの挙動を計算する。また、複数のサルコメアモデル処理部１４ｂのそれぞれは、複数のサルコメアモデル１３ｂのそれぞれに対応し、各サルコメアモデル処理部１４ｂは、対応するサルコメアモデル１３ｂの挙動を計算する。制御部１４は、演算処理部であるコアを複数有するマルチコアのＣＰＵ（Central  Processing  Unit）を１つ又は複数有する。若しくは、制御部１４は、演算処理部であるコアを１つ有するシングルコアのＣＰＵを複数有してもよい。心筋モデル処理部１４ａ、サルコメアモデル処理部１４ｂ、表示制御部１４ｃ及び受付部１４ｄの各部は、後述の生体シミュレーションプログラムを各コアが実行することにより、各コア上で実現される。
【００３２】
サルコメアモデル処理部１４ｂは、第１の計算部１５ａ、第２の計算部１５ｂを有する。
【００３３】
ここで、制御部１４で実行される処理の一例について図７を参照して説明する。図７は、実施例に係る生体シミュレーション処理の手順を示すフローチャートである。生体シミュレーション処理は、例えば、入力部１１から制御部１４に、生体シミュレーション処理を実行する指示が入力された場合に実行される。なお、生体シミュレーション処理では、心筋モデル１３ａの有限要素解析が時間ΔＴ間隔で実行される。また、サルコメアモデルのモンテカルロシミュレーション処理では、サルコメアモデル１３ｂのモンテカルロ法を用いたシミュレーションが、時間区分［Ｔ，Ｔ＋ΔＴ］がｎ個のステップの細かい時間区分に分割された時間Δｔ（ΔＴ／ｎ）間隔で実行される。ΔＴの値としては、例えば、２．５ミリ秒を採用することができる。また、Δｔの値としては、例えば、５マイクロ秒を採用することができる。モンテカルロ法では、遷移率と時間区分Δｔの積が１を超えないことが要請されるので、上記のような細かい時間刻みとなる。一方で、有限要素解析も細かい時間刻みで計算した方が精度の上から好ましいが、陰解法における線形求解に時間がかかるため上記のような時間区分ΔＴが適当である。
【００３４】
図７に示すように、心筋モデル処理部１４ａ及びサルコメアモデル処理部１４ｂは、シミュレーション処理を実行する（Ｓ１０１）。図８は、実施例に係るシミュレーション処理の手順を示すフローチャートである。図８に示すように、第１の計算部１５ａは、時間幅［Ｔ，Ｔ＋ΔＴ］における初期化を行う（Ｓ２０１）。例えば、Ｓ２０１で、第１の計算部１５ａは、変数ＳｕｍＬＩＲに「０」を設定する。また、Ｓ２０１で、第１の計算部１５ａは、変数ＸＢＤに「０」を設定する。また、Ｓ２０１で、第１の計算部１５ａは、変数ＸＢＤ０に「０」を設定する。また、Ｓ２０１で、第１の計算部１５ａは、フラグｆｌａｇＤ０に「１」を設定する。また、Ｓ２０１で、第１の計算部１５ａは、変数ＬＳ０に「Ｌｓ（Ｔ）」を設定する。ここで、Ｌｓ（Ｔ）は、時刻Ｔにおけるミオシンアーム２４の伸びのうちフィラメント間のすべり速度による寄与分のみを抜き出したものである。この変数は、処理開始時に「０」に設定され、後述するように毎回の有限要素ステップの終了後に正しい値に更新される。
【００３５】
図９を参照して、ミオシンアーム２４の伸びについての計算方法の一例について説明する。図９は、ミオシンアームの伸びについての計算方法の一例について説明するための図である。図９の例は、時刻ｔｂでミオシンヘッド２２がアクチンフィラメント２１に結合し、時刻Ｔに至るまでの間、結合が維持されている場合を示す。また、図９の例は、ミオシンヘッド２２が時刻ｔｂから時刻Ｔまでの間で首振り運動をした場合を示す。図９の例に示す場合において、第１の計算部１５ａは、下記の式（５）を用いて、Ｌｓ（Ｔ）を計算する。
【数７】
2014149649.tif 000009
ただし、ＬＩＮＩＴは、結合時のミオシンアーム２４の初期伸びである。また、ＬＲＯＴは、結合後のミオシンヘッド２２の首振り運動による伸びである。
【００３６】
図８に戻り、第１の計算部１５ａは、変数ｋの値を「１」に設定する（Ｓ２０２）。続いて、第１の計算部１５ａ及び第２の計算部１５ｂは、モンテカルロシミュレーション処理を実行する（Ｓ２０３）。なお、Ｓ２０３で実行されるモンテカルロシミュレーション処理は、時刻（Ｔ＋ｋΔｔ）における処理である。図１０は、実施例に係るモンテカルロシミュレーション処理の手順を示すフローチャートである。図１０に示すように、第１の計算部１５ａは、ミオシンヘッド２２の状態を判定し、結合状態であるか否かを判定する（Ｓ３０１）。例えば、Ｓ３０１では、第１の計算部１５ａは、ＰＸＢ、ＸＢＰｒｅＲ、ＸＢＰｏｓｔＲの何れかである場合には、結合状態であると判定し、ＮＸＢである場合には、結合状態でないと判定する。
【００３７】
結合状態でない場合（Ｓ３０１；Ｎｏ）には、第１の計算部１５ａは、Ｓ３０３へ進む。一方、結合状態である場合（Ｓ３０１；Ｙｅｓ）には、第１の計算部１５ａは、各種の前処理を行う（Ｓ３０２）。例えば、Ｓ３０２では、第１の計算部１５ａは、変数ＸＢＤの値を１つインクリメントする。また、Ｓ３０２では、第１の計算部１５ａは、フラグｆｌａｇＤ０の値が「１」である場合には、変数ＸＢＤ０の値を１つインクリメントする。そして、Ｓ３０２では、第１の計算部１５ａは、下記の式（６）を用いて、ミオシンアームの伸びＬを計算する。
【数８】
2014149649.tif 000010
ただし、ＬＩＲは、上述したＬＩＮＩＴとＬＲＯＴとの和を示す。ここで、ＬＩＮＩＴは結合状態への遷移直後のミオシンアームの伸びであり、例えば熱揺らぎを考慮してミオシンアームの弾性エネルギーから定まるボルツマン分布をもとに計算される。
【００３８】
そして、第１の計算部１５ａは、乱数を生成し、生成した乱数を用いて、ミオシンヘッド２２の状態を計算する（Ｓ３０３）。これにより、第１の計算部１５ａは、ミオシンヘッド２２の確率的な挙動を計算することができる。
【００３９】
続いて、第２の計算部１５ｂは、ミオシンヘッド２２の状態が遷移したか否かを判定する（Ｓ３０４）。ミオシンヘッド２２の状態が遷移しない場合（Ｓ３０４：Ｎｏ）には、第２の計算部１５ｂは、Ｓ３０６へ進む。
【００４０】
一方、ミオシンヘッド２２の状態が遷移した場合（Ｓ３０４：Ｙｅｓ）には、第２の計算部１５ｂは、状態遷移処理を行う（Ｓ３０５）。例えば、Ｓ３０５では、第２の計算部１５ｂは、ミオシンヘッド２２の状態が、非結合状態ＮＸＢから結合状態ＰＸＢへ遷移した場合には、変数ＬＩＮＩＴの値を変数ＬＩＲに設定し、フラグｆｌａｇＤ０の値を「０」に設定する。また、Ｓ３０５では、第２の計算部１５ｂは、ミオシンヘッド２２の状態が、結合状態ＰＸＢ、ＸＢＰｒｅＲ、ＸＢＰｏｓｔＲのいずれかから、非結合状態ＮＸＢへ遷移した場合には、変数ＬＩＲ、変数ＬＳ０及び変数ＸＢｄのそれぞれの値を「０」に設定する。また、Ｓ３０５では、第２の計算部１５ｂは、ミオシンヘッド２２の状態が、結合状態ＰＸＢ、ＸＢＰｒｅＲ、ＸＢＰｏｓｔＲのいずれかの状態のままで、ミオシンヘッド２２が首振り運動を行い、回転した場合には、次の処理を行う。すなわち、第２の計算部１５ｂは、変数ＬＩＲの値と、ミオシンヘッド２２の回転により伸びたミオシンアーム２４の長さＸＳＷＩＮＧとの和（ＬＩＲ＋ＸＳＷＩＮＧ）を、変数ＬＩＲに設定する。
【００４１】
そして、第２の計算部１５ｂは、後処理を実行し（Ｓ３０６）、処理結果を内部メモリに格納し、リターンする。例えば、Ｓ３０６では、第２の計算部１５ｂは、変数ＳｕｍＬＩＲの値と、変数ＬＩＲの値との和（ＳｕｍＬＩＲ＋ＬＩＲ）を変数ＳｕｍＬＩＲに設定する。また、Ｓ３０６では、第２の計算部１５ｂは、変数ＳｕｍＸＢＤの値と、変数ＸＢＤの値との和（ＳｕｍＸＢＤ＋ＸＢＤ）を変数ＳｕｍＸＢＤに設定する。上述した処理が有限要素解析における１ステップの間にｎステップ回行われることで、変数ＳｕｍＸＢＤには、有限要素解析における１ステップの間における、ミオシンヘッド２２の結合状態の連続回数が設定される。また、変数ＳｕｍＬＩＲには、有限要素解析における１ステップの間におけるＬＩＲの和が設定され、後述するようにこれらの和は有限要素解析におけるアクティブ応力を計算する際に用いられる。
【００４２】
図８の説明に戻り、第２の計算部１５ｂは、変数ｋの値が、ｎを超えたか否かを判定する（Ｓ２０４）。変数ｋの値がｎを超えていない場合（Ｓ２０４；Ｎｏ）には、第２の計算部１５ｂは、変数ｋの値を１つインクリメントし（Ｓ２０５）、Ｓ２０３に戻る。一方、変数ｋの値がｎを超えた場合（Ｓ２０４；Ｙｅｓ）には、心筋モデル処理部１４ａは、有限要素解析を実行する（Ｓ２０６）。図１１は、実施例に係る有限要素解析の手順を示すフローチャートである。図１１に示すように、心筋モデル処理部１４ａは、１ミオシンアームの平均諸量を計算する（Ｓ４０１）。例えば、Ｓ４０１では、心筋モデル処理部１４ａは、下記の式（７）を用いて、有限要素解析の１ステップの間における各ミオシンアーム２４の伸びの平均値ＬＡＶＲを計算する。
【数９】
2014149649.tif 000011
【００４３】
また、Ｓ４０１では、心筋モデル処理部１４ａは、下記の式（８）を用いて、有限要素解析の１ステップの間におけるミオシンアーム２４の力の平均値ＦＭＨ、及び、剛性の平均値ＫＭＨを計算する。
【数１０】
2014149649.tif 000012
ここで、Ｗａｒｍは、ミオシンアーム２４の伸びＬに対する関数として与えられたバネの弾性エネルギーである。
【００４４】
また、下記の式（９）は、与えられた線維方向のストレッチの変分δλに対するミオシンアームの仮想的な仕事δＷＭＨを表したものである。
【数１１】
2014149649.tif 000013
【００４５】
そして、心筋モデル処理部１４ａは、下記の式（１０）に基づき、心筋の連続体モデル上での任意の歪みの変分δEに対して、有限要素解析の１ステップにおける単位体積あたりの心筋モデル１３ａがした仮想的な仕事とミオシンアームの集合がした仮想的な仕事が等しくなるようにアクティブ応力Ｓａｃｔｉｖｅを計算する。以下、その方法の一例について説明する。
【数１２】
2014149649.tif 000014
式（１０）において、ｉｍは、ミオシンフィラメント２３上のミオシンヘッド２２のインデックスであり、ｉｓは、要素１３ａ＿１内のサルコメアモデル１３ｂのインデックスである。また、ＲＳＡは、筋肉に対するサルコメアが占める割合であり、ＳＡ０は、１サルコメアモデル１３ｂあたりの断面積である。また、Ｓａｃｔｉｖｅは、連続体モデルのアクティブ応力テンソル、δＥはＧｒｅｅｎ−Ｌａｇｒａｎｇｅ歪みテンソルの変分であり、この二つのテンソルの積が連続体の仮想仕事を表す。
【００４６】
心筋モデル処理部１４ａは、式（１０）に基づき下記の式（１１）を用いて、単位断面積あたりの心筋モデル１３ａの収縮力Ｆを計算する（Ｓ４０２）。
【数１３】
2014149649.tif 000015
【００４７】
また、Ｓ４０２では、心筋モデル処理部１４ａは、下記の式（１２）を用いて、単位断面積あたりの心筋モデル１３の剛性Ｋを計算する。
【数１４】
2014149649.tif 000016
【００４８】
ここで、下記の式（１３）と（１４）はＧｒｅｅｎ−Ｌａｇｒａｎｇｅ歪みテンソルの変分と線維方向のストレッチλの変分の関係を示すものである。式（１３）より、式（１１）で計算したＦを用いてアクティブ応力テンソルＳａｃｔｉｖｅを下記の式（１５）に示すように計算すれば式（１０）が成立することがわかる。
【数１５】
2014149649.tif 000017
【数１６】
2014149649.tif 000018
【数１７】
2014149649.tif 000019
【００４９】
それゆえ、心筋モデル処理部１４ａは、式（１５）を用いて、Ｓａｃｔｉｖｅを計算する（Ｓ４０３）。そして、有限要素解析をＮｅｗｍａｒｋ−β法のような陰解法で実行する場合には、心筋モデル処理部１４ａは、陰解法で用いられる剛性行列を、下記の式（１６）及び式（１７）を用いて計算する（Ｓ４０４）。
【数１８】
2014149649.tif 000020
【数１９】
2014149649.tif 000021
【００５０】
これらの式は、線維方向ストレッチλおよびその時間微分とＧｒｅｅｎ−Ｌａｇｒａｎｇｅ歪みテンソルＥおよびその時間微分の関係式（１３）および式（１４）から導かれたものである。Ｓ４０４では、心筋モデル処理部１４ａは、式（１０）により算出したアクティブ応力Ｓａｃｔｉｖｅから等価節点力を、式（１６）および式（１７）をもとに剛性行列を計算し、これらをすべての要素にわたって重ね合わせることによりＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復における線形方程式の右辺ベクトルと係数行列を作成する。残差ベクトルに対応する右辺ベクトルが十分小さければ（Ｓ４０５肯定）、Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復を終了して処理結果を内部メモリに格納しリターンする。そうでなければ（Ｓ４０５否定）、心筋モデル処理部１４ａは、線形方程式を解き、得られた解をもとに連続体モデルの変位ベクトル、速度ベクトル、加速度ベクトルを更新し（Ｓ４０６）、Ｓ４０１へ戻り、次のＮｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復に移る。ただし、これらベクトルの時間ステップＴ＋ΔＴにおける初期値は、これらベクトルの時刻Ｔにおける値とする。式（７）にみるように各ミオシンアームの時間幅［Ｔ，Ｔ＋ΔＴ］における平均伸びＬＡＶＲを時刻Ｔ＋ΔＴにおける線維方向ストレッチλの時間微分から計算しているので、収縮力の計算において連続体モデルの運動のフィードバックが強連的にかかり、剛性行列に対しても時間幅［Ｔ，Ｔ＋ΔＴ］での架橋運動の影響が適切に反映される結果となり、安定した計算が可能となる。
【００５１】
Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ反復が収束したら有限要素解析を終了し、図８に戻り、心筋モデル処理部１４ａは、下記の式（１８）を用いて、変数Ｌｓを次の時間ステップでの計算のために更新する（Ｓ２０７）。この更新によりＬＳには時刻Ｔ＋ΔＴまでにバネの伸びＬの中ですべり速度が寄与した分が正しく設定し直される。
【数２０】
2014149649.tif 000022
【００５２】
受付部１４ｄは、サルコメアモデル１３ｂに定義されるＴ／Ｔユニット２０の状態と遷移、及び、ミオシンヘッド２２の状態と遷移の定義を受け付ける。そして、受付部１４ｄは、受け付けた内容を各種のモデルに反映させる。図１２〜図１８は、実施例に係る受付部が受け付ける情報を説明するための図である。図１２の例に示すように、受付部１４ｄは、ユーザ、例えば、心臓を医学的に研究している研究者などによって、入力部１１から、Ｔ／Ｔユニット２０の状態と遷移を受け付ける画面を表示させる指示が入力されると、次の処理を実行する。すなわち、受付部１４ｄは、Ｔ／Ｔユニット２０の状態と遷移を受け付ける画面３０を表示部１２に表示させる。図１２の例は、画面３０を介して、「Ｃａ−ｏｆｆ」及び「Ｃａ−ｏｎ」の２つの状態が定義された場合を示す。また、「Ｉｎｉｔｉａｌ  Ｓｔａｔｅ」の横のチェックボックスにチェックを入れると、チェックが入れられた状態が、生体シミュレーション処理開始時の状態となる。図１２の例は、生体シミュレーション処理開始時の状態が、「Ｃａ−ｏｆｆ」であることを示す。
【００５３】
また、図１３の例に示すように、受付部１４ｄは、ユーザによって、入力部１１から、ミオシンヘッド２２の状態と遷移を受け付ける画面を表示させる指示が入力されると、次の処理を実行する。すなわち、受付部１４ｄは、ミオシンヘッド２２の状態と遷移を受け付ける画面３１を表示部１２に表示させる。図１３の例は、画面３１を介して、「Ｎ＿ＸＢ」、「Ｐ＿ＸＢ」、「ＸＢ＿ＰｒｅＲ」及び「ＸＢ＿ＰｏｓｔＲ」の４つの状態が定義された場合を示す。また、「Ｉｎｉｔｉａｌ  Ｓｔａｔｅ」の横のチェックボックスにチェックを入れると、チェックが入れられた状態が、生体シミュレーション処理開始時の状態となる。図１３の例は、生体シミュレーション処理開始時の状態が、「Ｎ＿ＸＢ」であることを示す。また、図１３の例に示す画面３１は、４つの状態のそれぞれが「Ｎｏｂｉｎｄｉｎｇ（非結合状態）」または「Ｂｉｎｄｉｎｇ（結合状態）」のいずれであるかを設定するためのラジオボタンが設けられている。図１３の例は、「Ｎ＿ＸＢ」及び「Ｐ＿ＸＢ」が非結合状態であり、「ＸＢ＿ＰｒｅＲ」及び「ＸＢ＿ＰｏｓｔＲ」が結合状態であることを示す。
【００５４】
また、図１４及び図１５の例に示すように、受付部１４ｄは、ユーザによって、入力部１１から、関数を定義するための画面を表示させる指示が入力されると、次の処理を実行する。すなわち、受付部１４ｄは、関数を定義するための画面３２や画面３３を表示部１２に表示させる。図１４の例は、画面３２を介して、状態が「Ｃａ−ｏｆｆ」である場合には、「Ｋ＿ＮＰ０」を返し、状態が「Ｃａ−ｏｎ」である場合には、「Ｋ＿ＮＰ１」を返すようなＴ／Ｔユニット２０の状態関数ｋｎｐが定義された場合を示す。また、図１５の例は、画面３３を介して、状態が「Ｎ＿ＸＢ」である場合には、「０」を返し、状態が「Ｐ＿ＸＢ」、「ＸＢ＿ＰｒｅＲ」及び「ＸＢ＿ＰｏｓｔＲ」である場合には、「１」を返すようなミオシンヘッド２２の状態関数ｎｇが定義された場合を示す。なお、定義された状態関数は、遷移率の定義において参照される。
【００５５】
また、図１６の例に示すように、受付部１４ｄは、ユーザが入力部１１を操作して、「Ｎ＿ＸＢ」から「Ｐ＿ＸＢ」へ向かう矢印が選択されると、「Ｎ＿ＸＢ」から「Ｐ＿ＸＢ」への遷移率を定義するための画面３４を表示部１２に表示させる。図１６の例は、「Ｎ＿ＸＢ」から「Ｐ＿ＸＢ」への遷移率を定義する場合に、ｇｅｔ（ＴＴ，ｋｎｐ）によりミオシンヘッド２２の直上のアクチン切片を配下とするＴ／Ｔユニット２０の状態に応じた値を返す関数ｋｎｐを用いる場合を示す。また、図１６の例は、ｇｅｔ（ＭＨ，ｎｇ，−１）及びｇｅｔ（ＭＨ，ｎｇ，１）によりミオシンヘッド２２の左および右隣のミオシンヘッド２２上の関数値ｎｇを参照し、ＧＡＭＭＡとして１より大きなパラメータが適用された場合を示す。これにより、近傍のミオシンヘッド２２間の協調性が表現される。さらに、ｘｉ＿ｏｖｅｒｌａｐ（）はサルコメア長からきまるアクチンフィラメント２１のオーバラップ状態を反映した関数であり、先の図６の例に示す関数χＲＡと関数χＬＡとを掛け合わせた関数である。
【００５６】
また、図１７の例に示すように、受付部１４ｄは、ユーザが入力部１１を操作して、「ＸＢ＿ＰｒｅＲ」から「ＸＢ＿ＰｏｓｔＲ」へ向かう矢印が選択されると、次の処理を行う。すなわち、受付部１４ｄは、ミオシンヘッド２２の回転を伴う、「ＸＢ＿ＰｒｅＲ」から「ＸＢ＿ＰｏｓｔＲ」への遷移率を定義するための画面３５を表示部１２に表示させる。図１７の例では、遷移がミオシンヘッド２２の回転を伴う遷移であれば、「Ｍｙｏｓｉｎ  ｈｅａｄ  ｓｗｉｎｇ」の横に存在するチェックボックスがユーザによりチェックされる。また、図１７の例では、回転にともなうミオシンアーム２４の伸びの増分が、Ｘ＿ＳＷＩＮＧとしてユーザにより入力される。遷移率の定義に際しては、ミオシンヘッド２２の元の状態「ＸＢ＿ＰｒｅＲ」のミオシンアーム２４の伸びＬが関数ａｒｍ＿ｌｅｎｇｔｈ（）で呼び出される。そして、遷移後の状態におけるミオシンアームの伸びａｒｍ＿ｌｅｎｇｔｈ（）＋Ｘ＿ＳＷＩＮＧでのミオシンアーム２４の弾性エネルギーがｓｐｒｉｎｇ＿ｅｎｅｒｇｙを呼び出されることにより求められる。そして、遷移後のミオシンヘッド２２の内部の化学エネルギーＥ＿ＸＢＰｏｓｔＲと足し合わせたトータルエネルギーから定まるボルツマン因子より遷移率ｒが定められる。
【００５７】
また、図１８の例に示すように、受付部１４ｄは、ユーザが入力部１１を操作して、「ＸＢ＿ＰｏｓｔＲ」から「Ｎ＿ＸＢ」へ向かう矢印が選択されると、次の処理を行う。すなわち、受付部１４ｄは、ミオシンヘッド２２の解離を伴う、「ＸＢ＿ＰｏｓｔＲ」から「Ｎ＿ＸＢ」への遷移率を定義するための画面３６を表示部１２に表示させる。図１８は、結合状態から非結合状態への遷移に対しては、ミオシンヘッド２２からＡＤＰが解離し、新たにＡＴＰが結合するものと仮定して、「ＡＴＰ  ｂｉｎｄｉｎｇ」と、「ＡＤＰ  ｒｅｌｅａｓｅ」のチェックボックスがユーザによりチェックされた場合を示す。本実施例に係る生体シミュレーション装置１０は、上述したようなユーザによる指定をもとして生体シミュレーション実行時に分子の消費量および生成量を計算することができる。
【００５８】
また、受付部１４ｄは、上述の入力で定義されたＴ／Ｔユニットとミオシンヘッドの状態と状態間の遷移情報をもとにモンテカルロシミュレーションを実行するコードを自動生成する。図１９は、ミオシンヘッドの状態と状態間遷移の定義をもとにモンテカルロステップを実行するコード（モンテカルロコード）を自動生成した結果を示す図である。図１９に示すように、かかるコードは、最初に[０，１]間の乱数ｒを生成し、現状態が何であるかに応じて処理を実行する。もし、現状態が結合状態であれば処理ｕｐｄａｔｅ＿ＸＢ（）により、先のステップＳ３０２で説明したように、ミオシンヘッド２２の変数（ＸＢＤ、Ｌ）を更新する。次に、コードは、Ｎ個ある状態のうち現状態に応じて、ｐｒｏｃｅｓｓ＿ｔｒａｎｓｉｔｉｏｎ＿ｓｔａｔｅ１、・・・、ｐｒｏｃｅｓｓ＿ｔｒａｎｓｉｔｉｏｎ＿ｓｔａｔｅＮのひとつを実行する。この処理では、乱数ｒの値に応じてｕｐｄａｔｅ＿ｓｔａｔｅ（）により状態を更新し、ｕｐｄａｔｅ＿ｖａｒｉａｂｌｅｓ（）により先のステップＳ３０５で説明したように、遷移の種類に応じた変数の更新を行う。
【００５９】
図１９に示すコードと、心筋モデル１３ａの有限要素コードとを、図８に示す処理のように組み合わせ、心筋モデル１３ａとの連成解析を通して心拍動のシミュレーションを実行するときの例を図２０Ａ及び図２０Ｂに示す。ここでは、左心室のモデルとして図２０Ｂに示す断面を回転して得られる回転対称な連続体を採用している。図２０Ａは、受付部１４ｄにより表示部１２に表示された画面を介して、「＃  ｏｆ  ｅｌｅｍｅｎｔｓ  ｉｎ  Ｒ」および「＃  ｏｆ  ｅｌｅｍｅｎｔｓ  ｉｎ  Ｌ」で貫壁方向および長手方向の要素数がユーザにより指定された場合を示す。また、図２０Ａは、「＃  ｏｆ  ｈｅａｒｔ  ｂｅａｔｓ」で、シミュレートする拍動回数がユーザにより指定された場合を示す。また、図２０Ａは、「α＿ｂａｓｅ」および「α＿ａｐｅｘ」で、図２０Ｂの例に示すように、線維方向の分布がユーザにより指定された場合を示す。
【００６０】
図２１Ａ〜図２１Ｃは、シミュレーション結果の一例を示す図である。図２１Ａの例に示すシミュレーション結果は、各拍動での拍出量(Ｅｊｅｃｔｉｏｎ)、拍出量とエネルギーの比（Ｅｊｅｃｔｉｏｎ  Ｆｒａｃｔｉｏｎ）、心筋モデルの成した仕事量(Ｍｕｓｃｌｅ  Ｗｏｒｋ)を出力したものである。さらに、図２１Ａの例に示すシミュレーション結果では、血液拍出で成した仕事量(Ｅｊｅｃｔｉｏｎ  ｗｏｒｋ)、ＡＴＰ消費量(ＡＴＰ ｃｏｎｓｕｍｐｔｉｏｎ)、効率(Ｅｆｆｅｃｉｅｎｃｙ)が出力されている。また、図２１Ｂの例に示すシミュレーション結果は、心室内容積と心室内圧力の時間変化および血液拍出に要する仕事率とＡＴＰ消費率を加水分解エネルギーで換算したものの時間変化を出力したものである。また、図２１Ｃの例に示すシミュレーション結果は、１拍動において、各要素が成した仕事の積算値の時間変化を示したものである。図２１Ａ〜図２１Ｃに示すシミュレーション結果から、本実施例に係る生体シミュレーション装置１０によれば、サルコメアモデル１３ｂにおける分子レベルの現象と心臓の拍動性能の関係を解析することが可能となる。
【００６１】
上述してきたように、本実施例に係る生体シミュレーション装置１０は、生体の筋肉に含まれるサルコメア内の複数のアクチン及び複数のミオシンの複数の状態、及び、複数の状態間の遷移率が定義されたサルコメアモデル１３ｂを用いて、次の処理を行う。すなわち、生体シミュレーション装置１０によれば、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態を計算する。そして、生体シミュレーション装置１０は、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態に基づいて、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動を計算する。それゆえ、生体シミュレーション装置１０によれば、複数のアクチン及びミオシンのそれぞれの確率的な挙動を計算することができる。したがって、生体シミュレーション装置１０によれば、精緻なシミュレーション結果を得ることができる。
【００６２】
また、本実施例に係る生体シミュレーション装置１０は、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動に基づいて、サルコメアの挙動を計算し、サルコメアの挙動に基づいて、筋肉の挙動を計算する。それゆえ、本実施例に係る生体シミュレーション装置１０によれば、筋肉の挙動を計算するシミュレーションと、サルコメアの挙動を計算するシミュレーションとを連成させることができる。
【００６３】
また、本実施例に係る生体シミュレーション装置１０は、Δｔ間隔で、複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動を計算する。そして、生体シミュレーション装置１０は、Δｔよりも長いΔＴ間隔で、ΔＴの間に計算された複数のサルコメアの挙動に基づいて、筋肉の挙動を計算する。このように、本実施例に係る生体シミュレーション装置１０は、ΔＴの間に複数回計算された複数のアクチンのそれぞれの挙動及び複数のミオシンのそれぞれの挙動を用いて、筋肉の挙動を計算する。そのため、生体シミュレーション装置１０によれば、筋肉の挙動を計算する時間刻みΔＴを有限要素解析の都合に合わせて選ぶことができ、心臓の拍動のような秒オーダーの現象を適度な時間で実行することが可能となる。このような時間間隔のギャップがあるにも拘わらす、式（１０）に示したようにサルコメアモデルと心筋モデルで仕事量を一致させているので精度の良い連成解析が可能となる。さらに、式（７）で示しているように、時刻Ｔ＋ΔＴでのストレッチ速度から時間幅［Ｔ，Ｔ＋ΔＴ］でのミオシンアームの平均伸びＬＡＶＲを算出しているので、サルコメアモデルと心筋モデルの仕事量が強く連成され適度な大きさの時間ステップΔＴでの安定した解析が可能となる。したがって、本実施例に係るシミュレーション装置１０によれば、精緻なシミュレーション結果を適切な計算時間内で得ることができる。
【００６４】
また、本実施例に係る生体シミュレーション装置１０が扱うサルコメアモデル１３ｂは、アクチンにミオシンが結合している場合には、当該ミオシンの所定範囲内に存在する他のミオシンが、アクチンに結合する状態となる確率が高くなるように定義されている。さらに、直上にあるアクチンフィラメントのオーバーラップ状態に応じて、個々のミオシンヘッドの状態遷移を定義している。すなわち、現実のサルコメアに類似する挙動を示すようにサルコメアモデル１３ｂが定義されている。そして、シミュレーション装置１０は、かかるサルコメアモデル１３ｂを用いて、複数のアクチンのそれぞれの状態及び複数のミオシンのそれぞれの状態を計算する。それゆえ、シミュレーション装置１０は、現実のサルコメアと類似するような状態を計算することができる。
【００６５】
また、生体シミュレーション装置１０は、ミオシンがアクチンに結合した際にミオシンアームが伸びた長さ、ミオシンヘッドの回転によってミオシンアームが伸びた長さ、ミオシンのすべり速度の積分値に基づいて、ミオシンが伸びた長さを計算する。
【００６６】
また、本実施例に係る生体シミュレーション装置１０は、サルコメアモデル１３ｂに定義される複数のアクチンの複数の状態及びアクチンの状態間の遷移、並びに、複数のミオシンの複数の状態及びミオシンの状態間の遷移を受け付ける。そして、生体シミュレーション装置１０は、受け付けた内容をもとに、モンテカルロステップを実行するモンテカルロコードを生成する。それゆえ、例えば、医学的に心臓の研究をしているがプログラミングが苦手な研究者などのユーザであっても、プログラミングを行うことなく、任意のサルコメアモデルのシミュレーションを行うことができる。
【００６７】
さて、これまで開示の装置に関する実施例について説明したが、本発明は上述した実施例以外にも、種々の異なる形態にて実施されてよいものである。
【００６８】
また、実施例において説明した各処理のうち、自動的に行われるものとして説明した処理の全部または一部を手動的に行うこともできる。また、本実施例において説明した各処理のうち、手動的に行われるものとして説明した処理の全部または一部を公知の方法で自動的に行うこともできる。
【００６９】
また、各種の負荷や使用状況などに応じて、実施例において説明した各処理の各ステップでの処理を任意に細かくわけたり、あるいはまとめたりすることができる。また、ステップを省略することもできる。
【００７０】
また、各種の負荷や使用状況などに応じて、実施例において説明した各処理の各ステップでの処理の順番を変更できる。
【００７１】
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的状態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。
【００７２】
［生体シミュレーションプログラム］
また、上記の実施例で説明した生体シミュレーション装置１０の各種の処理は、あらかじめ用意されたプログラムをパーソナルコンピュータやワークステーションなどのコンピュータシステムで実行することによって実現することもできる。そこで、以下では、図２２を用いて、上記の実施例で説明した生体シミュレーション装置と同様の機能を有するプログラムを実行するコンピュータの一例を説明する。図２２は、生体シミュレーションプログラムを実行するコンピュータを示す図である。
【００７３】
図２２に示すように、コンピュータ３００は、マルチコアの複数のＣＰＵ３１０、ＲＯＭ３２０、Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ（ＨＤＤ）３３０、ＲＡＭ３４０を有する。これら３１０〜３４０は、バス３５０を介して接続される。
【００７４】
ＲＯＭ３２０には、ＯＳなどの基本プログラムが記憶されている。また、ＨＤＤ３３０には、上記の実施例１で示す心筋モデル処理部１４ａ、サルコメアモデル処理部１４ｂ、第１の計算部１５ａ、第２の計算部１５ｂと同様の機能を発揮する生体シミュレーションプログラム３３０ａが予め記憶される。なお、生体シミュレーションプログラム３３０ａについては、適宜分離しても良い。
【００７５】
そして、ＣＰＵ３１０が、生体シミュレーションプログラム３３０ａを、ＨＤＤ３３０から読み出して実行する。
【００７６】
なお、上記した生体シミュレーションプログラム３３０ａについては、必ずしも最初からＨＤＤ３３０に記憶させておく必要はない。
【００７７】
例えば、コンピュータ３００に挿入されるフレキシブルディスク（ＦＤ）、ＣＤ−ＲＯＭ、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの「可搬用の物理媒体」に生体シミュレーションプログラム３３０ａを記憶させておく。そして、コンピュータ３００がこれらから生体シミュレーションプログラム３３０ａを読み出して実行するようにしてもよい。  
【００７８】
さらには、公衆回線、インターネット、ＬＡＮ、ＷＡＮなどを介してコンピュータ３００に接続される「他のコンピュータ（またはサーバ）」などに生体シミュレーションプログラム３３０ａを記憶させておく。そして、コンピュータ３００がこれらから生体シミュレーションプログラム３３０ａを読み出して実行するようにしてもよい。
【符号の説明】
【００７９】
１０        生体シミュレーション装置
１３ａ      心筋モデル
１３ｂ      サルコメアモデル
１４ａ      心筋モデル処理部
１４ｂ      サルコメアモデル処理部
１５ａ      第１の計算部
１５ｂ      第２の計算部
【図１】
2014149649.tif 000023
【図２】
2014149649.tif 000024
【図３】
2014149649.tif 000025
【図４】
2014149649.tif 000026
【図５】
2014149649.tif 000027
【図６】
2014149649.tif 000028
【図７】
2014149649.tif 000029
【図８】
2014149649.tif 000030
【図９】
2014149649.tif 000031
【図１０】
2014149649.tif 000032
【図１１】
2014149649.tif 000033
【図１２】
2014149649.tif 000034
【図１３】
2014149649.tif 000035
【図１４】
2014149649.tif 000036
【図１５】
2014149649.tif 000037
【図１６】
2014149649.tif 000038
【図１７】
2014149649.tif 000039
【図１８】
2014149649.tif 000040
【図１９】
2014149649.tif 000041
【図２０Ａ】
2014149649.tif 000042
【図２０Ｂ】
2014149649.tif 000043
【図２１Ａ】
2014149649.tif 000044
【図２１Ｂ】
2014149649.tif 000045
【図２１Ｃ】
2014149649.tif 000046
【図２２】
2014149649.tif 000047