(57)【要約】【課題】シミュレーションにおけるモデルのパラメータを効率的、かつ有効に決定する。【解決手段】モデルパラメータ決定プログラムによって、解析対象となる系のモデルにおける変数のデータの入力を受け取るステップと、受け取ったデータと系のモデルとに基づいて、該モデルのパラメータに関する制約式を生成するステップと、該制約式に対して限定記号消去アルゴリズムを適用して、解の存在の有無、および／または実行可能なパラメータの存在範囲を求めるステップとを計算機に実行させる。【選択図】図１
2010033615.tif 000002
【特許請求の範囲】
【請求項１】
系に対するシミュレーションを行う計算機に使用されるプログラムであって、
解析対象となる系のモデルに基づき行なわれたシミュレーション結果に基づく変数のデータを入力装置を介して受け取り、第一のファイルに格納し、または実験の結果に基づく変数のデータを入力装置を介して受け取り、第二のファイルに格納するステップと、
前記モデルの、実行可能範囲を決定したいパラメータの指定を受け取るステップと、
前記モデルにおけるモデル式或いは変数に対応する誤差変数と、前記第一のファイルに格納されたシミュレーション結果に基づく変数のデータまたは前記第二のファイルに格納された実験結果に基づく変数のデータと、前記系のモデルとに基づいて、前記指定された実行可能範囲を決定したいパラメータと前記誤差変数とこれらの誤差変数の範囲を表す変数とに関する制約式を生成するステップと、
該制約式に対して限定記号消去アルゴリズムを適用して、前記誤差変数の範囲を表す変数の可能領域を獲得するステップと、
前記獲得した誤差変数の範囲を表す変数の可能領域のデータを誤差変数の範囲を表す変数のデータとして前記制約式に入力し、限定記号消去アルゴリズムを適用して、解の存在の有無、および／または実行可能なパラメータの存在範囲を求めるステップと、
該求められた解の存在の有無、および／または実行可能なパラメータの存在範囲を出力装置に出力するステップと、
を計算機に実行させるモデルパラメータ決定プログラム。
【請求項２】
系に対するシミュレーションにおけるモデルパラメータを決定する装置であって、
解析対象となる系のモデルに基づき行なわれたシミュレーションまたは実験の結果に基づく変数のデータの入力を受け取るデータ入力手段と、
前記データ入力手段で受け取ったシミュレーション結果に基づく変数のデータを格納する第一のファイルと、
前記データ入力手段で受け取った実験の結果に基づく変数のデータを格納する第二のファイルと、
前記モデルの、実行可能範囲を決定したいパラメータの指定を受け取り、前記モデルにおけるモデル式或いは変数に対応する誤差変数と、前記第一のファイルに格納されたシミュレーション結果に基づく変数のデータまたは前記第二のファイルに格納された実験結果に基づく変数のデータと、前記系のモデルとに基づいて、前記指定された実行可能範囲を決定したいパラメータと前記誤差変数とこれらの誤差変数の範囲を表す変数とに関する制約式を生成する定式化部と、
前記定式化部によって生成された制約式に対して限定記号消去アルゴリズムを適用して、前記誤差変数の範囲を表す変数の可能領域を獲得し、前記獲得した誤差変数の範囲を表す変数の可能領域のデータを誤差変数の範囲を表す変数のデータとして前記制約式に入力し、限定記号消去アルゴリズムを適用して、解の存在の有無、および／または実行可能なパラメータの存在範囲を求めるパラメータ決定手段と、
該パラメータ決定手段で求められた解の存在の有無、および／または実行可能なパラメータの存在範囲を出力する出力装置とを備えることを特徴とするモデルパラメータ決定装置。
【請求項３】
系に対するシミュレーションにおけるモデルパラメータを決定する方法において、
計算機が、
解析対象となる系のモデルに基づき行なわれたシミュレーション結果に基づく変数のデータを入力装置を介して受け取り、第一のファイルに格納し、または実験の結果に基づく変数のデータを入力装置を介して受け取り、第二のファイルに格納し、
前記モデルの、実行可能範囲を決定したいパラメータの指定を受け付け、
前記モデルにおけるモデル式或いは変数に対応する誤差変数と、前記第一のファイルに格納されたシミュレーション結果に基づく変数のデータまたは前記第二のファイルに格納された実験結果に基づく変数のデータと、前記系のモデルとに基づいて、前記指定された実行可能範囲を決定したいパラメータと前記誤差変数とこれらの誤差変数の範囲を表す変数とに関する制約式を生成し、
該制約式に対して限定記号消去アルゴリズムを適用して、前記誤差変数の範囲を表す変数の可能領域を獲得し、
前記獲得した誤差変数の範囲を表す変数の可能領域のデータを誤差変数の範囲を表す変数のデータとして前記制約式に入力し、限定記号消去アルゴリズムを適用して、解の存在の有無、および／または実行可能なパラメータの存在範囲を求め、
該求められた解の存在の有無、および／または実行可能なパラメータの存在範囲を出力装置に出力する処理を行うことを特徴とするモデルパラメータ決定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、理工学分野における様々なモデルに対する実行可能パラメータの決定方式に係り、さらに詳しくは、例えば生体系シミュレーションのモデルと、変数データとに基づいて作成される制約式に対して、限定記号消去アルゴリズムを適用して実行可能なパラメータの存在範囲などを求めるためのモデルパラメータ決定プログラム、および決定装置に関する。
【背景技術】
【０００２】
本発明は、理工学の各種の分野において一階述語論理式、すなわち限定記号付の式によって表現される広範な範囲の応用問題を取り扱うものであるが、ここでは生体系シミュレーションにおける応用を主要な例として、本発明の内容を説明する。
【０００３】
生体や細胞に関する各種の系、すなわち生体系のモデル、例えば解糖系反応などのモデルに対してシミュレーションを行い、生体系の解析や設計を行う分野において、モデルのパラメータを効率的に決定することが望まれている。
【０００４】
すなわち本発明で扱う１つの問題は、例えばハイブリッドペトリネット（ＨＰＮ）によってモデル化された、対象となる生体系のモデルに対応して数値シミュレーション、または実測により生体系の変数の時系列データの値が与えられた時に、モデルに基づいて生体系の性質、あるいは本質を把握することである。具体的には、例えばモデルのパラメータとしての反応係数の実行可能範囲を求めること、モデルのパラメータのいくつかの実測値がわかっている時に残りのパラメータの実行可能領域を求めること、モデルの解の存在を確認すること、生体系に対する新しい知見を発見することである。
【０００５】
近年バイオ研究が盛んになるにつれて、このような要求を効率的に満たす有効な解決手法の開発が望まれている。そのような手法は、産業上の生体系設計過程、例えば疾病発症メカニズムの解明、予防医療、テーラーメイド医療などで、極めて広範な範囲において必要とされるものである。しかしながら従来においては、このような要求を効率的に満たす有効な手法は存在せず、生体系のモデルを構築して数値シミュレーションを行うためのツール、例えばゲノミックオブジェクトネット（ＧＯＮ）、ビジュアルオブジェクトネット（ＶＯＮ）、Ｅ−ＣＥＬＬなどを用いて、パラメータの値を様々に設定し、数値シミュレーションを繰り返し、トライアルアンドエラーの方式で望ましいパラメータの値を決定したり、系の性質を推測したりする方法が行われており、実際の作業としては大変に困難で、コストのかかる作業となっている。
【０００６】
このようなモデル化に用いられるハイブリッドペトリネット、シミュレーションのためのツールとしてのゲノミックオブジェクトネットや、ビジュアルオブジェクトネットについて、これらを紹介する文献や論文として、後掲の非特許文献１〜３がある。しかしながらこれらの文献の内容については、本発明と直接の関連がないため、その説明を省略する。
【０００７】
図２０は従来における生体系シミュレーション方式の説明図である。同図においては、例えば解糖系に対してハイブリッドペトリネットなどを用いてモデルが作成され、そのモデルに対してシミュレーション用ツールを用いて様々なパラメータの値について数値シミュレーションが繰り返し行われ、トライアルアンドエラーの方式で生体系の設計や、生体系の本質の理解が行われている。
【０００８】
しかしながらこのような方法では、パラメータの値の決定が設計者の知識や経験に強く依存したり、とりあえずの解で妥協しがちになったり、存在しない解を追って無意味な繰り返しシミュレーションを行ったりしてしまうという問題点があり、また複数の要求を同時に満たすパラメータの決定が非常に困難であり、生体系の本質の理解や、新しい知見の抽出が非常に困難であるという問題点があった。具体的には、同じパラメータで健康な状態と病気の状態を表すモデルが必要とされているが、未だにそのようなモデルが実現されていないという問題点があった。
【０００９】
次に本発明では、後述するように、例えば生体系を一階述語論理式で記述される制約問題として扱い、限定記号消去法によって解の存在の有無や、実行可能パラメータの存在範囲などを求めるが、この限定記号消去法を制御工学などの分野に応用した従来技術として後掲の特許文献１がある。
【００１０】
この文献では、半正定値計画（ＳＤＰ）問題、または拡張半正定値計画（ＥＳＤＰ）問題として与えられた問題を定式化し、その問題を一階述語論理式に変換し、解析的に解を求める技術が開示されている。しかしながら、本発明が対象とする生体系シミュレーションなどでは、解析的に解を求めるよりも、実行可能なパラメータの範囲などを求めることが重要であり、この文献の技術をそのまま応用することは困難であるという問題点があった。
【先行技術文献】
【特許文献】
【００１１】
【特許文献１】特開平１１−３２８２３９号公報「制御系解析・設計装置および制御系解析・設計処理方法」
【非特許文献】
【００１２】
【非特許文献１】ｈｔｔｐ：／／ｗｗｗ．ｇｅｎｏｍｉｃｏｂｊｅｃｔ．ｎｅｔ／ｍｅｍｂｅｒ３／ｉｎｄｅｘ．ｈｔｍｌ，「Ｇｅｎｏｍｉｃ  Ｏｂｊｅｃｔ  Ｎｅｔ  Ｐｒｏｊｅｃｔｓ」
【非特許文献２】ｈｔｔｐ：／／ｗｗｗ．ｇｅｎｏｍｅ．ｉｂ．ｓｃｉ．ｙａｍａｇｕｃｈｉ−ｕ．ａｃ．ｊｐ／〜ａｔｓｕｓｈｉ／ｐｈａｇｅ／ｋａｒｕｉ００．ｐｄｆ，「ハイブリッドペトリネットによるλファージの誘発と逆調節機構の表現」
【非特許文献３】ｈｔｔｐ：／／ｗｗｗ．Ｓｙｓｔｅｍｔｅｃｈｎｉｋ．ｔｕ−ｉｌｍｅｎａｕ．ｄｕ／〜ｄｒａｔｈ／ｖｉｓｕａｌ＿Ｅ．ｈｔｍ，「Ｖｉｓｕａｌ  Ｏｂｊｅｃｔ  Ｎｅｔ  ＋＋」
【非特許文献４】穴井宏和「Ｑｕａｎｔｉｆｉｅｒ  Ｅｌｉｍｉｎａｔｉｏｎ−アルゴリズム・実装・応用−」数式処理  Ｖｏｌ．１０，Ｎｏ．１，ｐｐ．３−１２（２００３）
【発明の概要】
【発明が解決しようとする課題】
【００１３】
本発明の課題は、前述の問題点に鑑み、対象となる生体系や制御系などのモデルや、数値シミュレーション、または実測によって得られる系における変数のデータの値に対応して、モデルのパラメータ、例えば反応係数の実行可能範囲を効率的に求め、その結果を用いて系に対する新しい知見を発見するための、シミュレーションにおけるモデルパラメータ決定プログラム、および決定装置を提供することである。
【課題を解決するための手段】
【００１４】
図１は本発明のモデルパラメータ決定プログラムの原理的な機能ブロック図である。このプログラムは生体系や制御系などのシミュレーションを行う計算機によって使用されるプログラムであり、図１の１において系のモデルにおける変数のデータの入力が受け取られ、２で受け取ったデータと系のモデルとに基づいて、そのモデルのパラメータに関する制約式が作成され、３でその制約式に対して限定記号消去（ＱＥ）アルゴリズムが適用されて、解の存在の有無や、実行可能なパラメータの存在範囲が求められる。
【００１５】
発明の実施の形態においては、図１の４で実行可能なパラメータの存在範囲をディスプレイ上に表示することも、また５で系に関する新しい知見として、限定記号消去アルゴリズムの適用結果に基づいたサンプル解、モデルの変数間の条件式、またはパラメータの間の関係式を含む知見を出力することもできる。
【００１６】
さらに実施の形態においては、前記変数のデータとして異なる２つの時刻における変数のデータを用いることもでき、また前記制約式が、前記モデルのパラメータに加えて、生体系などのシミュレーションモデルの変数に対応する誤差変数、または制約式内の方程式に対応する誤差変数を含むこともできる。
【００１７】
次に実施の形態においては、図１に原理的な機能ブロック図を示したプログラムの内容に対応するモデルパラメータ決定方法を用いることができる。
さらに実施の形態においては、モデルパラメータ決定装置は、生体系などのモデルにおける変数のデータの入力を受け取るデータ入力手段と、受け取ったデータと系のモデルとに基づいて、該モデルのパラメータに関する制約式を生成する制約式生成手段と、該制約式に対して限定記号消去アルゴリズムを適用して、解の存在の有無および／または実行可能なパラメータの存在範囲を求めるパラメータ決定手段とを備える。またこの装置は、求められた実行可能なパラメータの存在範囲を表示する表示手段をさらに備えることもできる。
【００１８】
すなわち本発明においては、まず与えられたモデルや系の変数のデータを用いて、モデルのパラメータが満たすべき制約式として、方程式や不等式の系が作成される。このモデルのパラメータについての制約問題は、一般的には非線形で非凸な制約問題であり、数値的な手法を用いて解くことは非常に困難である。そこで制約問題を解く手法としての限定記号消去（ＱＥ）アルゴリズムを適用することによって、非線形で非凸な制約問題を扱うことが可能となる。さらにパラメータをそのまま扱うことができるため、モデルのパラメータの実行可能領域を、例えば等式や不等式の論理和や論理積の形式で求めることができる。また解がない場合にも、それを正確に判定することができる。
【発明の効果】
【００１９】
本発明によれば、モデルの全てのパラメータについてその実行可能領域が得られ、また解の存在も確認することができる。系に対して複数の要求、例えば生体系における濃度の範囲、グルコースの流入速度、後述する解糖系のモデルにおける乳酸とピルビン酸の流出速度などについての要求があった場合にも、それぞれの要求に対応するパラメータの実行可能領域を重ね合わせることによって、複数の要求を同時に満たすパラメータを求めることもできる。さらにパラメータが領域として求められるため、パラメータの系への寄与など、系の本質の理解が深まり、系に関する新しい知見の発見にもつながる。
【００２０】
また本発明は系のパラメータ設計手法や装置として用いることもでき、単なる数値シミュレーションと比べて設計手順が系統的となり、よりよい解を選択することができ、正確な解を得ることができ、また複数の仕様の設計が容易であるという特徴を持ち、系に対するロバストな設計手法や装置を提供することができる。
【図面の簡単な説明】
【００２１】
【図１】本発明のモデルパラメータ決定プログラムの原理的な機能ブロック図である。
【図２】限定記号消去（ＱＥ）法の概要の説明図である。
【図３】制約問題に対するＱＥ法の適用の説明図である。
【図４】制約問題に対するＱＥ法の適用例を示す図である。
【図５】本発明における生体系シミュレーション方式の全体説明図である。
【図６】図５において本発明に特有の処理を実行するモデルパラメータ決定装置の機能構成ブロック図である。
【図７】図６における定式化部における制約式生成方法の説明図である。
【図８】解糖系に対するモデルを示す図である。
【図９】図８における変数、パラメータ、および反応速度を示す図である。
【図１０】図８に対応する無酸素解糖系を示す図である。
【図１１】時刻ｔにおける解糖系のモデルを示す図である。
【図１２】時刻ｔ’における解糖系のモデルを示す図である。
【図１３】制約式の定式化処理フローチャートである。
【図１４】実施例２におけるパラメータの実行可能領域を示す図である。
【図１５】実施例４におけるパラメータの実行可能領域を示す図である。
【図１６】実施例５におけるパラメータの実行可能領域を示す図である。
【図１７】機械系モデルの例としてのばね−質量−ダンパ系を示す図である。
【図１８】図１７の系に対する実験とシミュレーションの結果を示す図である。
【図１９】実施例７におけるパラメータの実行可能領域を示す図である。
【図２０】従来における生体系シミュレーション方式の説明図である。
【発明を実施するための形態】
【００２２】
本発明は、系に対するシミュレーションの実行時に、系のモデルのパラメータについて、例えば制約条件式を構成し、限定記号消去（ＱＥ）法を適用することによって、対象となる系のパラメータの実行可能範囲を求めたり、解そのものの存否を決定するものであり、まずこの限定記号消去法についてその概要を説明する。
【００２３】
多くの産業上の問題や数学の問題は方程式、不等式、限定記号、ブール演算記号などからなる式として記述される。このような式は一階述語論理式と呼ばれ、限定記号消去（ＱＥ）のアルゴリズムは、与えられた一階述語論理式に対して、限定記号のない等価な式を構成するためのアルゴリズムである。
【００２４】
この限定記号消去（Ｑｕａｎｔｉｆｉｅｒ  Ｅｌｉｍｉｎａｔｉｏｎ）法についてはその概要を紹介した前掲の非特許文献４がある。
図２はこのＱＥ法の概要の説明図である。同図において入力は多項式や不等式などを用いた一階述語論理式であり、出力は限定記号のないパラメータの実行可能領域であり、全ての変数が限定されている場合にはそのｔｒｕｅ、またはｆａｌｓｅ、すなわち解が存在するかしないか、また存在する場合にはサンプルの解も出力として得ることができる。このような問題は決定問題といわれる。
【００２５】
図２において変数ｘについてｘ2＋ｂｘ＋ｃ＝０という限定のついた問題に対して、限定記号のない等価な式としてｂ2−４ｃ≧０という式が得られている。
いくつかの変数に対して限定記号がない場合には、ＱＥアルゴリズムによって与えられた一階述語論理式と等価な限定記号のない式が得られるが、得られた式は限定記号がついていない残りの変数の可能な範囲を記述するものである。そのような範囲が存在しない場合、ｆａｌｓｅが出力される。このような問題は一般限定記号消去問題と呼ばれる。
【００２６】
次にＱＥアルゴリズムを用いた制約問題の解法とその簡単な具体例について図３、図４を用いて説明する。図３は制約問題の一般的解法の説明図である。同図において制約条件式として、ｆi（ｘ1，．．，ｘn）に関する不等式や等式などが与えられる。そのような条件式を満足する変数ｘ1，．．，ｘnを見つける問題がＱＥアルゴリズムを用いて解かれ、出力として実行可能性、すなわちｔｒｕｅ、またはｆａｌｓｅが出力され、また実行可能範囲が得られ、利点として任意のブーリアン結合を用いることができることなどが示されている。
【００２７】
図４においては具体例が示されており、例１ではｘとｙとの両方について限定記号がついているため、ＱＥアルゴリズムによってｔｒｕｅであることと、サンプル解とが出力されている。例２においてはｘだけに限定記号がついているため、ＱＥアルゴリズムによって、もう１つの変数としてのｙの実行可能範囲が出力されている。
【００２８】
限定記号消去法、すなわちＱＥアルゴリズムは単なる制約問題だけではなく、いわゆる最適化問題、すなわちある制約条件のもとで目的関数を最適化する目的関数付き制約問題の解法にも適用可能であるが、本実施形態では後述するように単なる制約問題に対してＱＥアルゴリズムを適用するために、最適化問題に関する説明は省略する。
【００２９】
続いて本発明が適用される第１の例として生体系のシミュレーションについて説明する。図５は本発明における生体系シミュレーション方式の全体説明図である。本方式においては、まず生物学の知識を用いて生体系のモデル化が行われる。このモデル化では前述のように、例えばハイブリッドペトリネットなどが用いられ、シミュレーションに適したモデル化が行われる。続いて生物学にＩＴ技術を応用した従来の技術としてのバイオシミュレーションが行われる。このバイオシミュレーションにおいては、例えば必要に応じてファイルに格納された実験値が用いられ、シミュレーションの結果はファイルに格納される。
【００３０】
バイオシミュレーションより先の部分が本発明の特有の処理であり、まずシミュレーション結果や実験値などを用いて、例えば実行可能範囲を決定したいパラメータとしての反応係数などを用いて制約問題、または最適化問題における制約式が定式化部によって作成される。作成された制約式は数式処理エンジンに与えられる。
【００３１】
数式処理エンジンによる処理としては、まず限定記号消去を行うＱＥ部によって、解の存在の有無と、反応係数などのパラメータの実行可能範囲などが求められ、その結果に対応して知見取得コンポーネント部によって、反応係数の絞込みや、反応係数の間の関係解析などが行われ、生体系のパラメータやその動作に関する知見が取得され、その結果はディスプレイ上に、例えばパラメータの実行可能範囲として表示される。そしてその結果は再びバイオシミュレーションによって利用され、シミュレーションの精度向上に寄与することになる。
【００３２】
図６は図５の定式化部と、数式処理エンジンとに対応する処理を実行するモデルパラメータ決定装置の構成ブロック図である。同図において計算機システム１０に対して、入力装置１１を介してシミュレーション結果が、入力装置１３を介して実験値が入力され、それぞれファイル１２および１４に格納される。
【００３３】
これらのファイル内のデータは、実行可能範囲を決定したいパラメータとしての反応係数などに対する定式化を行う定式化部１５に入力される。定式化部１５は、後述するように実験値などに対する誤差を変数として扱うために誤差変数を導入する誤差変数導入部、ある時刻ｔにおける基質に対する酵素の反応の速度を求める時刻ｔ反応速度計算部、制約問題あるいは最適化問題としての制約式を作成する制約式生成部、後述するように変数、例えば基質の濃度の変数や誤差変数に対して制約を付け加えるために、決定したい濃度変数と誤差変数を含む式をＥとして抽出するＥ抽出部、およびＥに対する制約付加部を備えている。
【００３４】
ＱＥ部１６、知見取得コンポーネント部１７による処理は、数式処理エンジン部１８を用いて行われる。ＱＥ部１６は解の有無を判定する解の存在判定部、および反応係数などのパラメータの可能範囲計算部を備え、また知見取得コンポーネント部１７は反応係数の絞込み部と、反応係数間の関係解析部とを備え、それぞれの処理結果は出力装置１９、例えばプリンタやディスプレイによって出力される。
【００３５】
図７は定式化部１５による制約式生成方法の説明図である。まず定式化部１５に対して生体系モデルと、少なくとも２つの時刻ｔおよびｔ’＝ｔ＋Δｔにおける各基質の濃度変数の数値データが与えられる。生体系モデルとしては、後述するようにモデル内の各基質の濃度の変数ｄ１からｄｎまでと、モデルのパラメータｐ１からｐｎまでとが、また各濃度変数の反応速度を与える式ｖ１〜ｖｎが与えられ、濃度変数の数値データとしては時刻ｔにおける各基質の濃度ｓ１からｓｎ、時刻ｔ’における数値データｓ１’からｓｎ’までが与えられる。
【００３６】
定式化部１５に対しては、パラメータの集合Ｐのうちで実行可能範囲を決定したい１つ以上のパラメータＰｆの指定が行われるとともに、その他のパラメータが後述するように実際の数値として与えられ、定式化部１５は以上の入力データを用いて制約式を作成し、ＱＥ部１６に出力する。
【００３７】
続いて、本実施形態における生体系問題へのＱＥアルゴリズムの適用について、具体例を用いてさらに説明する。図８はほとんどの生物に存在する糖分解の基本的代謝経路としての解糖系のモデルの一部であり、図９は図８のモデルにおける変数、パラメータ、および反応速度式を示す。解糖系は一般にＤ−グルコース、すなわちブドウ糖からエネルギーを作り出す経路であり、その化学的詳細は図１０に示される。
【００３８】
図８においては、図１０の解糖系に対応するモデルが概略的に示されている。図８の一番左上の二重丸は、図１０のＤ−フルクトース  ６−リン酸に相当し、Ｄ−フルクトース
１，６−ビスリン酸（ＦＢＰ）からグリセルアルデヒド  ３−リン酸（ＧＡＰ）、１，３−ビスホスホグリセリン酸（ＢＰＧ）の順序で反応が進行し、ピルビン酸（Ｐｙｒ）、および乳酸（Ｌａｃ）が形成されるまでの経路と、途中までのフィードバック経路とが示されている。
【００３９】
図８において二重丸は図１０との対応で説明したようにある物質、すなわち酵素の作用対象となる基質を表し、二重丸内の数字はその基質の濃度を表す。濃度の単位は、例えばｍｏｌ／ｌであるが、ここではある値を１として正規化した数値を示している。
【００４０】
例えば図８の一番左上のＦ６Ｐ、すなわちＤ−フルクトース  ６−リン酸の現在の濃度は１６であり、それに対して矢印の先のＰＦＫ、すなわち６−ホスホフルクトキナーゼという酵素を作用させることによってＦＢＰ、すなわちＤ−フルクトース  １，６−ビスリン酸が得られ、その反応速度がＦ６Ｐの濃度（変数）をａとして、ａ／（０．１＋ａ）であることが示されている。またＦＢＰの濃度がｒ＝２である状態に対して、反応速度がＴ２＝ｋｓ×ｒをもつある酵素を作用させると、グリセルアルデヒド  ３−リン酸（ＧＡＰ）が得られることが示されている。
【００４１】
さらに図８の右上にある二重丸はグルコースを表し、これらのグルコースの一部が必要に応じて解糖系において用いられる。ここでは解糖系への基本的な入力として０．２の濃度をもつグルコースが用いられ、そのグルコースに対して反応速度Ｔ０＝ｓ１をもつ酵素の作用によってＦ６Ｐが得られることが示されている。そのグルコースの右隣の二重丸は前述の反応速度Ｔ２の式内のｋｓの値を与える。なお図８はある計算開始時刻ｔ＝０におけるモデルを示し、モデル内の多くの基質の濃度はまだ０である。
【００４２】
図８において（  ）内の変数は各基質の濃度であり、Ｆ６Ｐに対するａ、ＦＢＰに対するｒ、ＧＡＰに対するｇａｐ、ＢＰＧに対するｂｐｇ、Ｐｙｒに対するｐｙｒ、Ｌａｃに対するｌａｃ、ＮＡＤＨ２に対するｎａｄｈ２、およびＮＡＤに対するｎａｄの６つが変数として用いられる。これらの変数はここでは濃度を表す変数であるが、一般には例えば制御系などにおいて一階常微分方程式に対応する状態変数に相当するものである。
【００４３】
各基質の間の反応速度を示すＴ０からＴ９は図９に示されているが、これらの式において用いられているｖｍ１からｖｍ６、およびｋｍ１からｋｍ６までがシミュレーションのパラメータである。これらのパラメータのうちｖｍはそれぞれの反応に対応する最大速度を示し、ｋｍは速度の式としてのミカエリス−メンテンの式におけるミカエリス定数である。なお、これらの反応速度の式の内、Ｔ０＝ｓ１は０．２であり、グルコースの流入速度を示し、Ｔ６は疲労物質としての乳酸の流出速度を示す。その他の式は酵素反応速度としてのミカエリス−メンテンの式であり、ひとつの基質の濃度Ｘで決定される式と、２つの基質の濃度Ｘ、Ｙによって決定される式との２種類がある。
【００４４】
図８においては合計１２個のパラメータ、すなわちｖｍ１からｖｍ６およびｋｍ１からｋｍ６までのパラメータの値がすべて図の中に記入されている。例えばＴ１の式におけるｖｍ１の値は１であり、またｋｍ１の値は０．１であることが、ＰＦＫの長方形の横に書かれたＴ１の式から読み取ることができる。ＱＥアルゴリズムでは、これらの１２個のパラメータのうち、少なくとも一部のパラメータの実行可能範囲を求めるために、少なくとも一部の変数に対する制約式を作成し、限定記号消去を行うことになる。
【００４５】
図１１および図１２は、図８では時刻ｔが０であった解糖系に対応する時刻ｔ＝５０秒後のモデルと、５０．０１秒後のモデルとを示す。当然のこととして図８、図１１および図１２で各基質の濃度は変化している。
【００４６】
本実施形態では図１１と図１２で示したように２つの時刻の各基質の濃度、すなわち変数の値を用いて制約式を作成することにする。すなわち、ｔ＝５０（図１１）における各濃度変数の値から、０．０１秒後の濃度変数の値を計算する式を各変数ごとに作成し、それらがｔ＝５０．０１（図１２）の時の濃度の値にそれぞれ一致するはずであるとして制約式を作成する。このように２つの時刻における変数の値を使うことによって、ＱＥアルゴリズムの計算が二次式の計算の程度となり、計算量を大幅に削減できることが判明した。
【００４７】
１つの時刻の濃度の式と、反応速度が正という条件を用いて制約式を作成し、パラメータの実行可能範囲を求めることもできるが、本実施形態では２つの時刻における変数の値を用いて制約式を作成する。
【００４８】
図８で説明した１２個のパラメータのうち、ｖｍ１からｖｍ５およびｋｍ１からｋｍ５までの１０個のパラメータについては、図８内の反応速度の式の中に示されている固定値を用いることにし、残りの２つのパラメータｖｍ６とｋｍ６との実行可能範囲を求めることにする。前述のように２つの時刻から求められた制約式のうちで、ｖｍ６とｋｍ６との２つのパラメータを含む式を取り出してＱＥアルゴリズムを適用し、ｖｍ６とｋｍ６との２つのパラメータの実行可能領域を求める。
【００４９】
しかしながら、数値シミュレーションや実験によるデータを変数の値として用いる場合には、数値計算の誤差や観測誤差が含まれるために、実際には解があっても、数式処理によって正確な計算を行うＱＥアルゴリズムにおいては解無しという結果が得られてしまうことがある。そこで本実施形態では、観測データにそれぞれ微小な誤差が含まれているものとし、誤差の変数を用いて制約式を作成し、ＱＥを適用することにする。実際にはそれぞれの濃度の変数に対する誤差変数として、各変数の先頭にｅをつけたｅａ，ｅｒ，ｅｇａｐ，ｅｂｐｇ，ｅｐｙｒ，ｅｌａｃ，ｅｎａｄｈ２，およびｅｎａｄの８つを用いて、各変数の観測データに誤差変数を加えて制約式を作成する。例えば変数ｇａｐに対しては、図１１における濃度の値の０．５２５０５にｅｇａｐを加えた結果を用いる。これらの制約式に、実行可能な範囲を求めたい２つのパラメータｖｍ６とｋｍ６を含む式と、各誤差変数に対して具体的な範囲を与える式、例えばｅｇａｐが−０．０１と０．０１との間の値であるという式を合わせて全体的な制約式を構成する。
【００５０】
図１３は以上のような処理を実行するための制約式定式化処理のフローチャートである。同図においてまずステップＳ１でそれぞれの変数、ここでは濃度変数に対して誤差変数ｅｓ１からｅｓｎが導入され、図７で説明した時刻ｔ’におけるそれぞれのデータの値に誤差変数が加えられ、ステップＳ２で各濃度変数の反応速度式に時刻ｔにおける各濃度変数データ、すなわちｓｔが代入され、時刻ｔにおける各変数の反応速度が計算され、その速度ｖ１ｔからｖｎｔまでが求められる。
【００５１】
続いてステップＳ３で濃度変数についての制約式の作成が行われる、すなわち時刻ｔ’における濃度変数の数値データ、例えば観測データに誤差変数が加えられたものと、時刻ｔにおける数値データに２つの時刻の差Δｔの間の反応の進行量に相当する濃度差を加算したものとの差が０であるという制約式が各濃度変数について作成される。そしてこれらのうち、図７で説明した実行可能範囲を決定したい１つ以上のパラメータＰｆ中のパラメータを含む式だけを取り出したものがＥとされ、ステップＳ４でＥの制約式にさらに濃度変数、および誤差変数についての制約が付け加えられてＱＥ部１６に出力される。濃度変数についての制約としては、物理的な当然の制約として各濃度変数の値が正であるという制約が加えられ、それぞれの濃度変数に対応する誤差変数に対しては絶対値としてその変数の範囲をきめるｅｍａｘが指定される。このｅｍａｘとしては後述するように具体的な数値を与えても良く、また変数として与えても良い。
【００５２】
図１３のフローチャートによって、前述のように２つのパラメータｖｍ６とｋｍ６の実行可能範囲を求めるための、各誤差変数の値を含む制約式が求められた。その制約式φを以下に示す。なおｖｍ６とｋｍ６との２つのパラメータは物理的に正の値をもつものとする。
【００５３】
φ(vm6,km6,egap,ebpg,enadh2,enad)＝
752390214379046110000*vm6+(642102601000000000000000*egap-4092911580889000000000)*km6
+75239021437904611000000*egap-479591052424276853779 = 0   and 
-2340591908156415067446250*vm6+(1997500928885875000000000000*ebpg
+5100547808120217000000000)*km6+234059190815641506744625000*ebpg
+597661846070320436514387 = 0   and 
-1814225129161183878996213129675*vm6+(1548290570465983757542500000000000*enadh2
+2265570353992712908187000000000)*km6+181422512916118387899621312967500*enadh2
+265470496720724021649561902057 = 0   and 
18142251291611838789962131296750*vm6+(15482905704659837575425000000000000*enad
-22036387311740735578853000000000)*km6+1814225129161183878996213129675000*enad
-2582135962040792861335770495383 = 0 and
vm6 >0 and km6 >0 and
-1/100 <egap <1/100 and 
-1/100 <ebpg <1/100   and 
-1/100 <enadh2 <1/100 and 
-1/100 <enad <1/100 
【００５４】
この制約問題にＱＥを適用する。これにより，数値計算あるいは観測の誤差を考慮に入れた形で，解の可能領域を求めたり，解の存否を決定したりすることが可能となる。
【００５５】
具体的には，実施例１としてまず以下のＱＥ問題を考える。
ψ１＝∃ｖｍ６  ∃ｋｍ６  ∃ｅｇａｐ  ∃ｅｂｐｇ  ∃ｅｎａｄｈ２  ∃ｅｎａｄ（φ）
全ての変数が限定記号付きなので，ＱＥアルゴリズムは，与えられた一階述語論理式ψ１がｔｒｕｅ（実数の解がある）か、ｆａｌｓｅであるかを決定する。この場合，答えは，ｔｒｕｅとなる。ｔｒｕｅの場合には，ＱＥは，サンプルの解も与えるが，この場合は，変数の具体的な数値ではなく変数間の条件式が得られる。以下がその解である。ここで，ｅｐｓｉｌｏｎ１は正の微小量（ｐｏｓｉｔｉｖｅ  ｉｎｆｉｎｉｔｅｓｉｍａｌ）である。ｉｎｆｉｎｉｔｙ１は，ある値を超えた大きな数である。
【００５６】
[true;
km6 = infinity1,
100000000000*epsilon1*infinity1 + 11717601100*epsilon1
vm6 = --------------------------------------------------------
117176011
1997500928885875000*epsilon1 - 5100547808120217
ebpg = -------------------------------------------------,
1997500928885875000
- 642102601000000*epsilon1 + 4092911580889
egap = ---------------------------------------------,
642102601000000
- 15482905704659837575425000*epsilon1 + 22036387311740735578853
enad = ------------------------------------------------------------------,
15482905704659837575425000
1548290570465983757542500*epsilon1 - 2265570353992712908187
enadh2 = -------------------------------------------------------------,]
1548290570465983757542500
【００５７】
すなわち、実施例１に対しては解が存在することがわかり、正の微小量ｅｐｓｉｌｏｎ１が決まれば、他の誤差変数の値も決まり、２つのパラメータｖｍ６とｋｍ６との間の関係も判明する。またｅｐｓｉｌｏｎ１と、ある値を超えた大きな数としてのｉｎｆｉｎｉｔｙ１については具体的な範囲を決定することも可能である。このような関係を表示することによって、生体系のシミュレーションにおける新しい知見を得ることも可能となる。
【００５８】
次に実施例２として、次のψ２で与えられるＱＥ問題を考える。
ψ２＝∃ｅｇａｐ  ∃ｅｂｐｇ  ∃ｅｎａｄｈ２  ∃ｅｎａｄ（φ）
実施例２においては、実施例１と比較して２つのパラメータｖｍ６とｋｍ６に対しては限定記号がないため、ＱＥアルゴリズムによって一階述語論理式ψ２と等価な式として、ｖｍ６およびｋｍ６についての限定記号のない式が得られるが、その式はｖｍ６とｋｍ６についての実行可能範囲を与えることになる。ＱＥアルゴリズムによって得られた結果を次に示す。
【００５９】
176865444358339111333103000000000000000000*km6^2
- 18142251291611838789962131296750000000000*km6*vm6
+ 41448774507305263302595803584266000000000*km6
- 2125836636910673029782829386411422264250*vm6
+ 2428401028802260556551431104671896121463 >0 and
132792669734857640175397000000000000000000*km6^2
+ 18142251291611838789962131296750000000000*km6*vm6
+ 31120230659142091857252721602734000000000*km6
+ 2125836636910673029782829386411422264250*vm6
+ 1823272245019085503014227668150948407037 >0 and 
17748476058652550483612000000000000000000*km6^2
- 1814225129161183878996213129675000000000*km6*vm6
+ 4159391251763815801291550063464000000000*km6
- 212583663691067302978282938641142226425*vm6
+ 243690437534990324867056242221754181052 >0 and
13217335350667124667238000000000000000000*km6^2
+ 1814225129161183878996213129675000000000*km6*vm6
+ 3097509264880919714693302455236000000000*km6
+ 212583663691067302978282938641142226425*vm6
+ 181476889847144281089509635060530271798 >0 and 
25075557096978967000000000000000000*km6^2
- 2340591908156415067446250000000000*km6*vm6
+ 5876507508453471007921274000000000*km6
- 274261223176647081663647531908750*vm6
+ 344292854226063255906182144679007 >0 and
14874461480738533000000000000000000*km6^2
+ 2340591908156415067446250000000000*km6*vm6
+ 3485860124172189261863726000000000*km6
+ 274261223176647081663647531908750*vm6
+ 204229592127230907421112919138493 >0 and 
10513937590889000000000000000000*km6^2
- 752390214379046110000000000000*km6*vm6
+ 2463962533606645927558000000000*km6
- 88162084036371465154867210000*vm6
+ 144358650470740106440320705569 >0 and
2328114429111000000000000000000*km6^2
+ 752390214379046110000000000000*km6*vm6
+ 545598323909538512442000000000*km6
+ 88162084036371465154867210000*vm6
+ 31965517602002823869413714431 >0 and
km6 >0 and   vm6 >0
【００６０】
これらの結果のうちで、上の８つの不等式は２つのパラメータｖｍ６とｋｍ６とを含むものであり、ｋｍ６についての２次不等式となっている。これらの不等式の左辺の２次式を因数分解することによって２つの１次式の積が得られるが、２つの１次式のうち一方は必ず正であることが判明し、したがって他方の１次式によって２つのパラメータｖｍ６とｋｍ６との実行可能範囲が決定される。これらの１次式に対応する１次不等式は次の８つの不等式となる。
【００６１】
G1=176865444358339111333103000000000*km6-18142251291611838789962131296750*vm6+20724387253652631651297901792133＞0
G2=132792669734857640175397000000000*km6+18142251291611838789962131296750*vm6+15560115329571045928626360801367＞0
G3=17748476058652550483612000000000*km6-1814225129161183878996213129675*vm6+2079695625881907900645775031732＞0
G4=13217335350667124667238000000000*km6+1814225129161183878996213129675*vm6+1548754632440459857346651227618＞0
G5=25075557096978967000000000*km6-2340591908156415067446250*vm6+2938253754226735503960637＞0
G6=14874461480738533000000000*km6+2340591908156415067446250*vm6+1742930062086094630931863＞0
G7=10513937590889000000000*km6-752390214379046110000*vm6+1231981266803322963779
＞0
G8=2328114429111000000000*km6+752390214379046110000*vm6+272799161954769256221＞0
【００６２】
図１４は実施例２の結果としてのｖｍ６とｋｍ６との実行可能範囲の表示結果である。同図においては、２つのパラメータｖｍ６とｋｍ６とがともに正であるという物理的に当然な事実と比較すると、８つの不等式のうちの一部の不等式によって、第１象限の一部が切り取られているだけに過ぎないという、実用的にはあまり有効でない結果が得られている。
【００６３】
シミュレーションにおいては２つのパラメータｖｍ６およびｋｍ６をさらに絞り込むことが必要となる。この２つのパラメータの絞込みを実施例３として説明する。この実施例３においては、２つのパラメータｖｍ６とｋｍ６の実行可能範囲をさらに絞り込むために、まず各濃度変数に対する誤差の範囲をどれだけ小さくできるかをＱＥアルゴリズムによって求めることにする。誤差の範囲を表す変数の絶対値をｅｍａｘとし、各濃度変数に対する誤差が＋ｅｍａｘと−ｅｍａｘとの間に入るという制約を追加して、次のψ３で表されるＱＥ問題を考える。すなわちこの実施例３におけるＱＥ問題は誤差の範囲がｅｍａｘ
で指定される以外は実施例１におけるＱＥ問題と同じである。
【００６４】
ψ3＝∃vm6 ∃km6 ∃egap ∃ebpg ∃enadh2 ∃enad（τ）
ここで、
τ(vm6,km6,egap,ebpg,enadh2,enad,emax)＝
(752390214379046110000*vm6+(642102601000000000000000*egap-4092911580889000000000)*km6+75239021437904611000000*egap-479591052424276853779= 0 and
-2340591908156415067446250*vm6+(1997500928885875000000000000*ebpg+5100547808120217000000000)*km6+234059190815641506744625000*ebpg+597661846070320436514387= 0 and
-1814225129161183878996213129675*vm6+(1548290570465983757542500000000000*enadh2+2265570353992712908187000000000)*km6+181422512916118387899621312967500*enadh2+265470496720724021649561902057= 0 and
18142251291611838789962131296750*vm6+(15482905704659837575425000000000000*enad-22036387311740735578853000000000)*km6+1814225129161183878996213129675000*enad-2582135962040792861335770495383= 0 and
vm6 >0 and   km6 >0 and emax >0 and
-emax <egap <emax and 
-emax <ebpg <emax and 
-emax <enadh2 <emax and 
-emax <enad <emax
【００６５】
ψ３にＱＥを適用すると、限定記号のついていない変数  ｅｍａｘ  の実行可能領域が以下のように得られる。
55941746971386000000*emax - 138482701358119397 ＞ 0
すなわち、ｅｍａｘ  が  138482701358119397/55941746971386000000 ≒  0.00247548  より大きければｖｍ６，ｋｍ６についての制約式が解を持つことがわかる。
【００６６】
そこで誤差の範囲としてｅｍａｘ＝０．００５として、ＱＥ問題ψ２、すなわち実施例２におけるＱＥ問題を解く場合を実施例４として説明する。実施例２と同様に２つのパラメータｖｍ６とｋｍ６とについての１次式を求めた結果を以下に示す。
【００６７】
R1=49725457917519961727989000000000*km6-9071125645805919394981065648375*vm6+5826630803923356128158418071879> 0
R2=27689070605779226149136000000000*km6+9071125645805919394981065648375*vm6+3244494841882563266822647576496> 0
R3=20014046412645263391799000000000*km6-3628450258322367757992426259350*vm6+2345166122602631922295336933789> 0
R4=10951764996674411759051000000000*km6+3628450258322367757992426259350*vm6+1283284135719735835697089325561> 0
R5=7544026226274796000000000*km6-1170295954078207533723125*vm6+883978900074263985118756> 0
R6=2443478418154579000000000*km6+1170295954078207533723125*vm6+286317054003943548604369> 0
R7=7303424585889000000000*km6-752390214379046110000*vm6+855786159613799908779 
>0
R8=882398575889000000000*km6-752390214379046110000*vm6+103395945234753798779< 
0
【００６８】
図１５は実施例４としての２つのパラメータｖｍ６とｋｍ６との実行可能範囲の表示結果である。同図において、ｖｍ６とｋｍ６との実行可能範囲が第１象限内において図１４よりかなり狭まった領域であることがハッチされた領域として示されている。
【００６９】
同様にｅｍａｘ＝０．００２５とした場合の、２つのパラメータｖｍ６とｋｍ６との実行可能範囲の絞込みを実施例５として説明する。実施例５においても、実施例４と同様にＱＥ問題ψ２に対してＱＥアルゴリズムが適用され、実施例４におけると同様に以下の１次不等式が得られた。
【００７０】
L1=121487303146780659034831000000000*km6-36284502583223677579924262593500*vm6+14235397569887505117652606639141> 0
L2=33341753899817716719419000000000*km6+36284502583223677579924262593500*vm6+3906853721724333672309524657609> 0
L3=24545187120630689208173000000000*km6-7256900516644735515984852518700*vm6+2876107116044079965594460737903> 0 
L4=6420624288688985942677000000000*km6+7256900516644735515984852518700*vm6+752343142278287792397965521447> 0 
L5=20188600260669809000000000*km6-4681183816312830134892500*vm6+2365619646218848406751899< 0
L6=213590971811059000000000*km6-4681183816312830134892500*vm6+25027738062433339305649> 0 
L7=2487655078389000000000*km6-752390214379046110000*vm6+291493498829515326279< 0
【００７１】
図１６は実施例５における２つのパラメータｖｍ６とｋｍ６との実行可能範囲の表示結果である。同図において、実行可能範囲は非常に狭まり、ほとんど一本の太い直線によって示されている。
【００７２】
以上で生体系についての実施例の説明を終了し、他の分野についての実施例を図１７に示す機械系を対象として説明する。同図において質量Ｍによる慣性力とダンバー（粘性定数Ｄ）の制御力とばね（係数ｋ）の力の和が外力ｆに等しく、次式が成立する。
【００７３】
【数１】
2010033615.tif 000003
速度ｖを位置ｘ（被制御変数）で置き換えて次式を得る。
【００７４】
【数２】
2010033615.tif 000004
簡単のため外力（強制項、操作変数）を０とし、また差分式とするために
【００７５】
【数３】
2010033615.tif 000005
として次の連立常微分方程式が得られる。
【００７６】
【数４】
2010033615.tif 000006
【００７７】
この系を対象とするＱＥ問題を考える。なおこの機械系は制御系、電子回路、粘弾性皮膜などの生態系において相似な問題としてよく利用される基本的なモデルを示すものである。
【００７８】
この機械系に対する実験データとシミュレーション結果を図１８に示す。これらの結果は前述のようにｆ＝０とし、またＭ＝１として、ｔ＝０〜２０の範囲でΔｔ＝０．１刻みとして、実験結果（右側）はｘとｘ’＝ｖとを、シミュレーション結果（左側）はｘを示す。
【００７９】
生体系に対すると同様に２つの時刻として図１８にイタリックで示したｔ＝４．９と５．０をとり、ｋとｄ（Ｄ）、誤差変数ｅ１、ｅ２、および誤差範囲の最大絶対値ｅｍａｘを含む制約式として次式が得られる。
【００８０】
μ（ｋ，ｄ，ｅ１，ｅ２，ｅｍａｘ）＝
(1-d/10)*(467990039/1000000000)-(k/10)*(-6865249/500000000)-(113172839/250000000)+ e1 = 0 and
(1/10)*(467990039/1000000000)+(-6865249/500000000)- (32334129/1000000000) +e2 = 0 and
-emax <e1 <emax and
-emax <e2 <emax and
emax >=0 
【００８１】
なおここで誤差変数ｅ１、ｅ２は上から２つの方程式のそれぞれに対する誤差を表わすものであり、前述の生体系の場合のように各変数に対する誤差を表わすものではない。このように誤差変数は各変数に対して設定することも、また制約式内の各方程式に対して設定することも可能である。
【００８２】
この制約式にＱＥを適用する。具体的には実施例６として次のＱＥ問題を考える。
ψ４＝∃ｋ∃ｄ ∃ｅ１ ∃ｅ２（μ）
ψ４ではｅｍａｘだけに限定記号がついていないため、ＱＥアルゴリズムはｅｍａｘの範囲を決定し、その他の変数については変数間の関係を与える。その結果を以下に示す。
【００８３】
10000000000*emax - 7343769 >0,
10000000000*emax - 467990039*epsilon1 + 13730498*k + 152986830
d =----------------------------------------------------------------,
467990039
10000000000*emax - 467990039*epsilon1
e1 =---------------------------------------,
10000000000
- 7343769
e2 = -------------
10000000000      ここでｅｐｓｉｌｏｎ１は正の微小量
【００８４】
実施例６の結果としてｅｍａｘの範囲が得られたので、その限界値よりわずかに大きいｅｍａｘの値を用いて制約式を作成し、ｅ１とｅ２に限定記号をつけてＱＥ問題とする場合を実施例７として最後に説明する。
その制約式は以下のようになる。
【００８５】
λ（ｋ，ｄ，ｅ１，ｅ２）＝
(1-d/10)*(467990039/1000000000)-(k/10)*(-6865249/500000000)-(113172839/250000000)+ e1 = 0 and
(1/10)*(467990039/1000000000)+(-6865249/500000000)- (32334129/1000000000) +e2 = 0 and
-7343770/10000000000< e1 <7343770/10000000000 and
-7343770/10000000000< e2 <7343770/10000000000 and 
k>0and d>0 
【００８６】
実施例７におけるＱＥ問題は
ψ５＝∃ｅ１  ∃ｅ２（λ）
であり、ＱＥアルゴリズムは次の結果を与える。
【００８７】
467990039*d - 13730498*k - 145643060 >0 and 
467990039*d - 13730498*k - 160330600 <0 and 
d >0  and k >0,
467990039*d - 13730498*k - 152986830              - 7343769
e1 = --------------------------------------,      e2 =-------------
10000000000                          10000000000
【００８８】
上の４つの不等式によってｋとｄ（Ｄ）の実行可能領域が与えられる。その範囲を図１９に示す。同図において、第１象限の２本の直線の間の狭い範囲としてｋとＤの実行可能領域が示されている。
【００８９】
（付記１）  系に対するシミュレーションを行う計算機によって使用されるプログラムにおいて、
解析対象となる系のモデルに基づき行なわれたシミュレーションまたは実験の結果に基づく変数データの入力を入力装置を介して受け取るステップと、
該受け取ったデータと前記系のモデルとに基づいて、定式化部が該モデルのパラメータに関する制約式を生成するステップと、
該制約式に対して限定記号消去アルゴリズムを適用して、限定記号消去部が解の存在の有無、および／または実行可能なパラメータの存在範囲を求めるステップとを計算機に実行させるためのモデルパラメータ決定プログラム。
【００９０】
（付記２）  前記制約式が、前記モデルのパラメータに加えて、前記モデルにおける変数に対応する誤差変数を含むことを特徴とする付記１記載のモデルパラメータ決定プログラム。
【００９１】
（付記３）  前記変数のデータとして、複数の時刻における変数のデータを用いることを特徴とする付記１記載のモデルパラメータ決定プログラム。
【００９２】
（付記４）  前記実行可能なパラメータの存在範囲を出力装置に表示するステップをさらに備えることを特徴とする付記１記載のモデルパラメータ決定プログラム。
【００９３】
（付記５）  前記制約式に対する限定記号消去アルゴリズムの適用結果に基づいて、サンプル解、モデルの変数間の条件式、またはパラメータ間の関係式を含む系に関する新しい知見を出力装置に出力するステップをさらに備えることを特徴とする付記１記載のモデルパラメータ決定プログラム。
【００９４】
（付記６）  前記制約式が、前記モデルのパラメータに加えて、該制約式内の各方程式に対応する誤差変数を含むことを特徴とする付記１記載のモデルパラメータ決定プログラム。
【００９５】
（付記７）  前記解析対象となる系が生体系であることを特徴とする付記１記載のモデルパラメータ決定プログラム。
（付記８）  前記解析対象となる系が、機械系を含み、常微分方程式で表わされる物理系であることを特徴とする付記１記載のモデルパラメータ決定プログラム。
【００９６】
（付記９）  系に対するシミュレーションを行う計算機によって使用されるプログラムにおいて、
系のモデルにおける変数のデータの入力を受け取るステップと、
該受け取ったデータと系のモデルとに基づいて、該モデルのパラメータに関する制約式を生成するステップと、
該制約式に対して限定記号消去アルゴリズムを適用して、解の存在の有無、および／または実行可能なパラメータの存在範囲を求めるステップとを計算機に実行させるためのモデルパラメータ決定プログラム。
【００９７】
（付記１０）  系に対するシミュレーションにおけるモデルパラメータを決定する装置において、
解析対象となる系のモデルに基づき行なわれたシミュレーションまたは実験の結果に基づく変数データの入力を受け取るデータ入力手段と、
該受け取ったデータと前記系のモデルとに基づいて、該モデルのパラメータに関する制約式を生成する制約式生成手段と、
該制約式に対して限定記号消去アルゴリズムを適用して、解の存在の有無、および／または実行可能なパラメータの存在範囲を求めるパラメータ決定手段とを備えることを特徴とするモデルパラメータ決定装置。
【００９８】
（付記１１）  前記求められた実行可能なパラメータの存在範囲を表示する表示手段をさらに備えることを特徴とする付記１０記載のモデルパラメータ決定装置。
【００９９】
（付記１２）  系に対するシミュレーションにおけるモデルパラメータを決定する方法において、
解析対象となる系のモデルに基づき行なわれたシミュレーションまたは実験の結果に基づく変数データの入力を入力装置を介して受け取り、
該受け取ったデータと前記系のモデルとに基づいて、定式化部が該モデルのパラメータに関する制約式を生成し、
該制約式に対して限定記号消去アルゴリズムを適用して、限定記号消去部が解の存在の有無、および／または実行可能なパラメータの存在範囲を求めることを特徴とするモデルパラメータ決定方法。
【産業上の利用可能性】
【０１００】
本発明は生体系シミュレーションやバイオ研究を必要とする産業は当然のこととして、一階述語論理式で記述されるシステムなどを取り扱うすべての産業において利用可能である。
【符号の説明】
【０１０１】
１０  計算機システム
１１，１３  入力装置
１２，１４  ファイル
１５  定式化部
１６  限定記号消去（ＱＥ）部
１７  知見取得コンポーネント部
１８  数式処理エンジン部
１９  出力装置
【図１】
2010033615.tif 000007
【図５】
2010033615.tif 000008
【図６】
2010033615.tif 000009
【図７】
2010033615.tif 000010
【図８】
2010033615.tif 000011
【図９】
2010033615.tif 000012
【図１０】
2010033615.tif 000013
【図１１】
2010033615.tif 000014
【図１２】
2010033615.tif 000015
【図１３】
2010033615.tif 000016
【図１４】
2010033615.tif 000017
【図１５】
2010033615.tif 000018
【図１６】
2010033615.tif 000019
【図１８】
2010033615.tif 000020
【図２０】
2010033615.tif 000021
【図２】
2010033615.tif 000022
【図３】
2010033615.tif 000023
【図４】
2010033615.tif 000024
【図１７】
2010033615.tif 000025
【図１９】
2010033615.tif 000026