(57)【要約】【課題】  個々の細胞中のｍＲＮＡ量が時間的に、また、細胞相互間においても、大きな揺らぎを持って変化する系に適用可能な、各細胞内に存在している、特定のｍＲＮＡの量（Ｘi）の経時的な変化を確率論的に記述するモデル化方法の提供。【解決手段】  確率論的に推定可能な形態で表記するモデルとして、  一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生し；一つの細胞内における該特定の遺伝子の転写事象発生一回に付随して、シャープなピーク状のｍＲＮＡ量の時間的な変化が生じ；一つの細胞内において、一定の時間内における該特定の遺伝子の転写事象発生回数は、ポワソン分布に類する確率に従い、また、前記確率と反しない範囲で、前後する二回の転写事象発生間における時間間隔は、実質的にランダムであるという条件を満たすモデルを用いる。【選択図】    なし【特許請求の範囲】
【請求項１】
同種の細胞複数個を含んでなる生細胞群において、前記細胞が有する特定の遺伝子の転写に起因して生成されるｍＲＮＡに関して、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を確率論的に推定可能な形態で表記するモデル化を行う方法であって、
該生細胞群に含まれる同種の細胞複数個のうち、前記特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞数は、時間的に大きな揺らぎを示し、かつ、
前記特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞複数においても、個々の細胞内における前記ｍＲＮＡ量を比較すると、細胞間で大きな揺らぎを示している生細胞群に対して、
その確率論的に推定可能な形態で表記するモデルとして、
前記生細胞群における、一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生し、
一つの細胞内における該特定の遺伝子の転写事象発生一回に付随して、前記ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化が生じ、
一つの細胞内において、一定の時間内における該特定の遺伝子の転写事象発生回数は、ポワソン分布に類する確率分布に従い、
また、前記確率分布と反しない範囲で、前後する二回の転写事象発生間における時間間隔は、実質的にランダムであるという条件を満たすモデルを用いる
ことを特徴とする、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法。
【請求項２】
同種の細胞複数個を含んでなる生細胞群において、前記細胞が有する特定の遺伝子の転写に起因して生成されるｍＲＮＡに関して、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を確率論的に推定可能な形態で表記するモデル化を行う方法であって、
前記生細胞群における、一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生し、
個々の細胞内における前記ｍＲＮＡ量の時間的な変動は、
前記一定の時間を短い時間間隔（Δｔ）に分割し、該時間間隔（Δｔ）内における前記ｍＲＮＡ量（Ｘi）の変化量：ΔＸi、
転写に起因するｍＲＮＡの生成過程に相当する速度項：ｋgenerate＊Ｘi
（但し、ｋgenerateは、ある範囲内で時間的に変動するが、その時間的な平均値＜ｋgenerate＞は、特定の値に集束するような、時間的に変動する速度定数である）、
酵素的分解に起因するｍＲＮＡの消失過程に相当する速度項：−ｋdegrade＊Ｘi
（但し、ｋdegradeは、ある範囲内で時間的に変動するが、その時間的な平均値＜ｋdegrade＞は、特定の値に集束するような、時間的に変動する速度定数である）、
と定義する際、
前記ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、
前記ｍＲＮＡ量（Ｘi）に依存する項として、前記二つの速度項の和：ｋgenerate＊Ｘi−ｋdegrade＊Ｘiを含み、
前記時間的に変動する速度定数ｋgenerateとｋdegradeは、
前記時間的な平均値＜ｋgenerate＞と＜ｋdegrade＞とは、＜ｋdegrade＞≧＜ｋgenerate＞の条件を満足し、かつ互いに独立な確率函数に従って、時間的に変動する函数とし、
個々の細胞に対して、独立に選択される該速度定数ｋgenerateとｋdegradeを割り振ることにより、確率論的に推定可能な形態の表記がなされる
ことを特徴とする、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法。
【請求項３】
個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、下記式（Ｉ）：
式（Ｉ）  ：  ΔＸi＝ｋgenerate＊Ｘi−ｋdegrade＊Ｘi
で示され、前記ｍＲＮＡ量（Ｘi）に依存する項である前記二つの速度項の和：ｋgenerate＊Ｘi−ｋdegrade＊Ｘiのみ含むものとする
ことを特徴とする、請求項２に記載の確率論的モデル化の方法。
【請求項４】
個々の細胞内における前記ｍＲＮＡ量（Ｘi）を該ｍＲＮＡ分子数と定義する際、該ｍＲＮＡ分子数の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、下記式（II）：
式（II）  ：  ΔＸi＝ｋgenerate＊Ｘi−ｋdegrade＊Ｘi＋Ｃ
（但し、Ｃは、ある範囲内で時間的に変動するが、その時間的な平均値＜Ｃ＞は、特定の値に集束するような、時間的に変動する定数項である）
かつ、
前記時間的な平均値＜ｋgenerate＞と＜ｋdegrade＞、ならびに＜Ｃ＞とは、
＜ｋdegrade＞≧＜ｋgenerate＞の条件を満足し、
前記時間的に変動する定数項Ｃは、
０≧＜ｋgenerate＞−＜ｋdegrade＞＋＜Ｃ＞の条件を満足するように、独立な確率函数に従って、時間的に変動する定数とされている、
ことを特徴とする、請求項２に記載の確率論的モデル化の方法。
【請求項５】
個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、前記式（Ｉ）において、
前記時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝Ａ＋α＋Ｒ1
ｋdegrade＝Ａ＋Ｒ2
（但し、
Ａは、Ａ＞０  の条件を満たす定数、
αは、（Ａ＋α）＞０≧α  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝０であり、標準偏差σ（Ｒ1）は、（Ａ＋α）＞３σ（Ｒ1）の条件を満たす正規分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝０であり、標準偏差σ（Ｒ2）は、（Ａ＋α）＞３σ（Ｒ2）の条件を満たす正規分布に従う確率変数である）
を選択してなる、下記式（１）：
式（１）：  ΔＸi＝（Ａ＋α＋Ｒ1）＊Ｘi−（Ａ＋Ｒ2）＊Ｘi
で示される
ことを特徴とする、請求項３に記載の確率論的モデル化の方法。
【請求項６】
個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、前記式（II）において、
前記時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝α＋Ｒ1
ｋdegrade＝β＋Ｒ2
前記時間的に変動する定数項Ｃとして、
Ｃ＝Ｒ3
（但し、
αは、α≧０  の条件を満たす定数、
βは、β≧０  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝λ1であり、ポワソン分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝λ2であり、ポワソン分布に従う確率変数、
Ｒ3は、平均値＜Ｒ3＞は、＜Ｒ3＞＝λ3であり、ポワソン分布に従う確率変数であって、
（β＋λ2）≧（α＋λ1）  ならびに  ０≧（α＋λ1）−（β＋λ2）＋λ3  の条件を満たす）
を選択してなる、下記式（２）：
式（２）：  ΔＸi＝（α＋Ｒ1）＊Ｘi−（β＋Ｒ2）＊Ｘi＋Ｒ3
で示される
ことを特徴とする、請求項５に記載の確率論的モデル化の方法。
【請求項７】
個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、前記式（II）において、
前記時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝α1＋α2＊Ｒ1
ｋdegrade＝β1＋β2＊Ｒ2
前記時間的に変動する定数項Ｃとして、
Ｃ＝Ｒ3
（但し、
α1は、α1＞０  の条件を満たす定数、α2は、α2＞０  の条件を満たす定数、
β1は、β1≧０  の条件を満たす定数、β2は、β2＞０  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝λ1であり、ポワソン分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝λ2であり、ポワソン分布に従う確率変数、
Ｒ3は、平均値＜Ｒ3＞は、＜Ｒ3＞＝λ3であり、ポワソン分布に従う確率変数であって、
（β1＋β2＊λ2）≧（α1＋α2＊λ1）  ならびに  ０≧（α1＋α2＊λ1）−（β1＋β2＊λ2）＋λ3  の条件を満たす）
を選択してなる、下記式（３）：
式（３）：  ΔＸi＝（α1＋α2＊Ｒ1）＊Ｘi−（β1＋β2＊Ｒ2）＊Ｘi＋Ｒ3
で示される
ことを特徴とする、請求項５に記載の確率論的モデル化の方法。
【請求項８】
同種の細胞複数個を含んでなる生細胞群において、前記細胞が有する特定の遺伝子の転写に起因して生成されるｍＲＮＡに関して、細胞あたりのｍＲＮＡ量平均を、該生細胞群から採取されるサンプル群を用いて、定量する方法であって、
前記生細胞群を構成する個々の細胞中における、前記特定の遺伝子に由来するｍＲＮＡ量は、時間的に揺らぎを示し、また、個々の細胞間においても、前記特定の遺伝子の発現量に揺らぎを示す際、
個々の細胞中における、前記特定の遺伝子に由来するｍＲＮＡ量の時間的な揺らぎの分布、あるいは、個々の細胞間における、前記特定の遺伝子の発現量の揺らぎの分布の大きさに基づき決定される、基準の細胞数に対して、該生細胞群から採取されるサンプル群を構成する細胞数を前記基準の細胞数より多くするサンプリング方法を用い、
前記細胞数を有するサンプル群において測定される、細胞あたりのｍＲＮＡ量平均を、該生細胞群全体における細胞あたりのｍＲＮＡ量平均とする
ことを特徴とする、生細胞群における、特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の定量方法。
【請求項９】
個々の細胞中における、前記特定の遺伝子に由来するｍＲＮＡ量の時間的な揺らぎの分布、あるいは、個々の細胞間における、前記特定の遺伝子の発現量の揺らぎの分布の大きさを、請求項１に記載の遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法を利用して、推定する
ことを特徴とする、請求項８に記載のｍＲＮＡ量の定量方法。
【請求項１０】
所定の時間内において、前記生細胞群における、前記特定の遺伝子の転写がなされている細胞数の比率の多寡に従って、前記基準の細胞数を決定する
ことを特徴とする、請求項８または９に記載のｍＲＮＡ量の定量方法。
【請求項１１】
所定の時間内において、前記生細胞群における、前記特定の遺伝子の転写がなされている細胞数の比率に反比例するように、前記基準の細胞数を決定する
ことを特徴とする、請求項１０に記載のｍＲＮＡ量の定量方法。
【請求項１２】
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１０，０００個以上である
ことを特徴とする、請求項１０または１１に記載のｍＲＮＡ量の定量方法。
【請求項１３】
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１００，０００個以上である
ことを特徴とする、請求項１０または１１に記載のｍＲＮＡ量の定量方法。
【請求項１４】
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１，０００，０００個以上である
ことを特徴とする、請求項１０または１１に記載のｍＲＮＡ量の定量方法。
【請求項１５】
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１０，０００，０００個以上である
ことを特徴とする、請求項１０または１１に記載のｍＲＮＡ量の定量方法。
【請求項１６】
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１００，０００，０００個以上である
ことを特徴とする、請求項１０または１１に記載のｍＲＮＡ量の定量方法。
【請求項１７】
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１，０００，０００，０００個以上である
ことを特徴とする、請求項１０または１１に記載のｍＲＮＡ量の定量方法。
【請求項１８】
請求項１〜７のいずれか一項に記載の遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法に基づき、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象を反映する、前記ｍＲＮＡ量の時間的変動の確率論的なシミュレーションを、コンピュータ上において実現するための演算処理用制御プログラム。
【請求項１９】
請求項１〜７のいずれか一項に記載の遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法に基づき、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象を反映する、前記ｍＲＮＡ量の時間的変動の確率論的なシミュレーションを、コンピュータ上において実現するための演算処理用制御プログラムを、前記コンピュータにより読み取り可能な形態で格納した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、所謂、システム・バイオロジーの分野において、生細胞群中における遺伝子産物の産生過程を取り扱う際に利用可能な解析手法に関する。特には、生細胞群中に含まれる、各細胞内において、各種のタンパク質をコードする遺伝子（ＤＮＡ）に基づき、対応するＲＮＡ分子の転写に起因する、ｍＲＮＡの生成過程と、その後のｍＲＮＡの分解過程とにより支配される、ｍＲＮＡ量の変動現象を確率論的に取り扱う際、利用可能な近似的なモデル化手法と、かかるモデル化に基づき、シミュレーションされるｍＲＮＡ量の変動の確率論的な予測と、実際に、全体の生細胞群中からサンプリングされた部分細胞群において観測されるｍＲＮＡ量の変動とを利用して、全体の生細胞群におけるｍＲＮＡ量の変動を高い確度で推定する方法に関する。
【背景技術】
【０００２】
一般に、生細胞内において、そのゲノム上にコードされているタンパク質の発現過程では、先ず、当該タンパク質をコードするＤＮＡに基づき、ＲＮＡ分子の転写がなされ、イントロン配列を含む場合には、スプライシング過程を経て、ｍＲＮＡの生成がなされる。次いで、ｍＲＮＡに必要に応じて、様々な化学的修飾がなされた後、リボソームにおいて、ｍＲＮＡの塩基配列に従って、対応するアミノ酸配列を有するペプチド鎖へと翻訳される。翻訳されたペプチド鎖は、適正にホールディングされ、同時に、必要に応じて、翻訳後修飾が施され、最終的に、所定の三次元構造を有するタンパク質分子となる。この一連の過程は、ウイルスを除く、生物全般の生命活動を規律するものであり、セントラル・ドクマと称されている。
【０００３】
この過程においては、タンパク質の生合成には、対応するｍＲＮＡの生成が必須であり、また、ｍＲＮＡの生成量は、その後、翻訳されるペプチド鎖の量を決定する主要な因子となっている。すなわち、対象となるタンパク質をコードするＤＮＡに基づく、ＲＮＡ分子転写過程の制御を介して、タンパク質産生の制御がなされている。ＲＮＡ分子転写過程における、ＲＮＡ分子の合成自体は、ＲＮＡポリメラーゼの酵素反応によるが、転写制御には、転写因子タンパク質などの制御因子が関与している。
【０００４】
転写促進がなされ、対応するｍＲＮＡの生成が進むと、それに付随して、タンパク質の産生も進む。その後、転写が抑制、停止されると、新たなｍＲＮＡの生成がなされず、また、タンパク質の産生も抑制、停止される。最終的には、一旦生成されたｍＲＮＡは、細胞内に存在する分解酵素の作用により、分解を受ける。一方、産生されたタンパク質も、その機能を終了すると、細胞内に存在する分解酵素の作用により、分解を受けるものの、ｍＲＮＡの分解と、対応するタンパク質の分解とは、必ずしも同期していない。
【０００５】
タンパク質は、対応する遺伝子の発現に由来する産物であり、新たな産生に伴う、該タンパク質量の増加は、遺伝子の発現量、例えば、転写量を反映する間接的な指標となる。また、新たに生成されるｍＲＮＡの増加量は、遺伝子に基づく、転写過程を直接的に反映する指標となる。この特徴を考慮して、遺伝子の発現量変化を追跡する際、対応するタンパク質量の増加に代えて、転写に起因して、生成されるｍＲＮＡ量を定量し、指標として利用している。
【０００６】
実際は、全体が均質な条件下に置かれている、同一種類の生細胞群を対象とする際、多数の細胞中から、その一部のサンプリングした上で、このサンプル群に含まれる複数個の細胞中に存在する目的のｍＲＮＡ量の総和を定量して、そのサンプル群に含まれる細胞数で除し、一細胞当たりの平均値を算定している。また、全体の生細胞群が均一な分布を有していると見なせる場合、該サンプル群の「一細胞当たりの平均値」を、全体の生細胞群における代表値として利用している。
【０００７】
一方、外見上は同一の環境下にある生細胞群について、その全体群から採取されたサンプル群について、全ｍＲＮＡ試料を調製して、各ｍＲＮＡについて、その含有量を定量し、「一細胞当たりの平均値」の対比を行うと、相当のバラツキを示すことが判明してきた。具体的には、ＤＮＡマイクロアレイなどを利用することで、一つのサンプル中に存在する、複数種のｍＲＮＡについて、そのｍＲＮＡ量を一度に測定することが可能な方法が開発され、現在では、非常に多種類のｍＲＮＡ量を同時に定量すること可能となってきている。その結果、遺伝子の種類に依って、そのｍＲＮＡ量の相違が非常に大きいことが判ってきている。この現象は、多くの生物種に共通して観測されており、その際、多数のサンプル群において、実測された各ｍＲＮＡ量の「一細胞当たりの平均値（ｋ）」と、その平均値（ｋ）が観測される頻度：確率関数Ｐ（ｋ）との間には、概ね、Ｐ（ｋ）∝ｋ-r（なお、ｒは、２程度）の関係が見出されることが報告されている（非特許文献１）。この種の変動を生じさせる要因としては、個々の細胞内における遺伝子の発現制御は、各細胞毎に独立してなされており、また、発現制御の機構自体、細胞内における種々の生体反応が関与しており、その両者が、複雑に寄与している結果と推定される。
【０００８】
例えば、細胞の発育条件などを変化させると、異なる発育条件下にある一群の細胞から採取される、各ｍＲＮＡ量の「一細胞当たりの平均値（ｋ）」は、発育条件に応じて、異なっている。しかし、実測された各ｍＲＮＡ量の「一細胞当たりの平均値（ｋ）」と、その平均値（ｋ）が観測される頻度：確率関数Ｐ（ｋ）との間には、概ね、Ｐ（ｋ）∝ｋ-r（なお、ｒは、２程度）の関係が保存されていることが報告されている。その際、外的な要因に付随して、ｍＲＮＡ量の減少（増加）が引き起こされているものに関して、ｍＲＮＡ量の減少（増加）比率の比較を行うと、対応する遺伝子がより高い発現状態にある際に測定されるｍＲＮＡ量の「一細胞当たりの平均値（ｋ）」が多いものほど、そのｍＲＮＡ量の減少（増加）比率が大きいという傾向が報告されている（非特許文献１）。この種の外的要因に付随して、複数種の遺伝子において、そのｍＲＮＡ量の変動が起こる機構は、個々の細胞内におけるこれらお遺伝子の発現制御は、互いに密接に絡み合っており、また、発現制御の機構自体、細胞内における種々の生体反応が関与しており、その両者が、複雑に寄与している結果と推定されてきた。換言すると、従来は、個々の遺伝子の発現、つまり、個々の遺伝子発現に由来する各ｍＲＮＡ量の変化は、その生物中で起こっている多種多様な生化学反応の結果、必然的に生じているという理解が主流である。この考え方に基づき、上記の外的環境の変化に伴い、個々の遺伝子発現量には影響が現れるにもかかわらず、実測された各ｍＲＮＡ量の「一細胞当たりの平均値（ｋ）」と、その平均値（ｋ）が観測される頻度：確率関数Ｐ（ｋ）との間には、依然として“べき乗関係”が存在する現象を説明する上では、個々の遺伝子の発現制御を司っている「ネットワークメカニズム」を解明することの必要性が指摘されている（非特許文献１）。
【非特許文献１】Ｐｒｏｃｅｅｄｉｎｇｓ  ｏｆ  ｔｈｅ  Ｎａｔｉｏｎａｌ  Ａｃａｄｅｍｙ  ｏｆ  Ｓｃｉｅｎｃｅｓ  Ｍａｒｃｈ  １６，  ２００４  Ｖｏｌ．  １０１  Ｎｏ．  １１  ３７６５−３７６９
【発明の開示】
【発明が解決しようとする課題】
【０００９】
例えば、細胞のセルサイクルに関連している、ハウスキーピング・タンパク質の発現を考慮すると、細胞分裂サイクルに対応して、各細胞内においては、遺伝子発現促進に付随して、転写の促進が起こると、対応するｍＲＮＡの量が急速に上昇し、その後、転写抑制がなされると、ｍＲＮＡの生成は停止し、一方、細胞内に存在するＲＮＡ分子の酵素的分解反応に伴って、ｍＲＮＡの残量は指数関数的に減衰していく。結果的に、一回の細胞分裂サイクルの期間内、限られた時間内にシャープなスパイク状のピークを示すような過渡的なｍＲＮＡ量変化が起こる。一方、細胞分裂自体は、必ずしも、周期的に起こる現象ではなく、従って、このシャープなスパイク状のｍＲＮＡ量変化が生じる時間間隔には大きなバラツキが見出される。さらには、最終的に産生されるハウスキーピング・タンパク質量は、ある許容範囲内で変動可能であり、すなわち、各スパイク状のｍＲＮＡ量変化のピーク高さも、ある範囲内での変動は許容されている。
【００１０】
このような、各細胞内における遺伝子発現が起こるまでの、時間間隔（ｔinterval）、ならびに、各回の遺伝子発現において、その転写量に主に依存するｍＲＮＡ量の極大値（Ｘmaximum）の双方ともに、ポワソン分布に類する確率関数に従って、変動する系と見なせる場合、各細胞内に存在している、特定のｍＲＮＡの量（Ｘi）の経時的な変化を確率論的に記述するモデルは、未だ提案されていない。加えて、このような、スパイク状のｍＲＮＡ量変化各回における、ｍＲＮＡ量の極大値（Ｘmaximum）、次回のスパイク状のｍＲＮＡ量変化が起こるまでの時間間隔（ｔinterval）は、散発的に発生する事象の頻度を記述する、ポワソン分布に類する確率分布に従って変動するものの、長期間にわたり、スパイク状のｍＲＮＡ量変化を積分すると、その積分値（＜Ｘi＞）となる頻度：確率関数Ｐ（＜Ｘi＞）との間には、概ね、Ｐ（＜Ｘi＞）∝＜Ｘi＞-r（ｒは、２程度）の関係が成立するような、特定のｍＲＮＡの量（Ｘi）の経時的な変化を確率論的に記述するモデルは、未だ提案されていない。
【００１１】
仮に、前述する要件を満足する、特定のｍＲＮＡの量（Ｘi）の経時的な変化を確率論的に記述するモデル化と、かかるモデルを近似的に表記でき、数値演算に適用可能なモデル函数が見出されると、全体が均質な条件下に置かれている、同一種類の生細胞群を対象とする際、個々の細胞における、特定のｍＲＮＡの量（Ｘi）の経時的な変化を個別に、該モデル函数を適用して、シミュレーションした上で、それらのシミュレーション計算結果を統合して、その全体群から採取されたサンプル群において、特定のｍＲＮＡの量を定量した際に、期待される「一細胞当たりの平均値」を推定することが可能となる。さらには、推定されるサンプル群において観測される「一細胞当たりの平均値」と、実測された「一細胞当たりの平均値」とを対比することで、全体の生細胞群における「代表値」を高い確度で推定することも可能となる。
【００１２】
本発明は、前記の課題を解決するもので、本発明の目的は、各細胞内における遺伝子発現が起こるまでの、時間間隔（ｔinterval）、ならびに、各回の遺伝子発現において、その転写量に主に依存するｍＲＮＡ量の極大値（Ｘmaximum）の双方ともに、離散的な系における散発的な事象の出現頻度を記述する、ポワソン分布に類する確率分布に従って、変動する系と見なせる場合、各細胞内に存在している、特定のｍＲＮＡの量（Ｘi）の経時的な変化を確率論的に記述するモデル化と、かかるモデルを近似的に表記でき、数値演算に適用可能なモデル函数を提供することにある。特には、本発明は、スパイク状のｍＲＮＡ量変化各回における、ｍＲＮＡ量の極大値（Ｘmaximum）、次回のスパイク状のｍＲＮＡ量変化が起こるまでの時間間隔（ｔinterval）は、離散的な系における散発的な事象の出現頻度の記述に適している、ポワソン分布に類する確率分布に従って変動するものの、長期間にわたり、スパイク状のｍＲＮＡ量変化を積分すると、その該積分値（＜Ｘi＞）となる頻度：確率関数Ｐ（＜Ｘi＞）との間には、概ね、Ｐ（＜Ｘi＞）∝＜Ｘi＞-r（ｒは、２程度）の関係が成立するような、特定のｍＲＮＡの量（Ｘi）の経時的な変化を確率論的に記述するモデルを近似的に表記でき、数値演算に適用可能なモデル函数を提供することにある。
【課題を解決するための手段】
【００１３】
本発明者は、上記の課題を解決すべく、鋭意研究を進めた結果、下記のような知見を得た。
【００１４】
全体が均質な条件下に置かれている、同一種類の生細胞群を対象とする際、各細胞内における、特定の遺伝子の転写に起因して生成されるｍＲＮＡ量が時間的に大きな揺らぎを示し、同時に、各細胞内で転写に起因して生成されるｍＲＮＡ量は、細胞相互においても、大きな揺らぎを示している系は、その生細胞群に含まれる同種の細胞複数個のうち、前記特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞数は、時間的に大きな揺らぎを示し、かつ、この特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞複数においても、個々の細胞内における前記ｍＲＮＡ量を比較すると、細胞間で大きな揺らぎを示している系として、取り扱えることを見出した。さらには、前記の系において、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を確率論的に推定可能な形態で表記するモデルとして、
・一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生する。
・一つの細胞内における該特定の遺伝子の転写事象発生一回に付随して、ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化が生じている。
・一つの細胞内において、一定の時間内における該特定の遺伝子の転写事象発生回数は、ポワソン分布に類する確率分布に従い、また、この確率分布と反しない範囲で、前後する二回の転写事象発生間における時間間隔は、実質的にランダムとなっている。
この３条件を満足するモデルを利用できることを見出した。
【００１５】
本発明者は、これらの一連の検討結果に基づき、本発明を完成するに至った。
【００１６】
すなわち、本発明にかかる遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法は、
同種の細胞複数個を含んでなる生細胞群において、前記細胞が有する特定の遺伝子の転写に起因して生成されるｍＲＮＡに関して、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を確率論的に推定可能な形態で表記するモデル化を行う方法であって、
該生細胞群に含まれる同種の細胞複数個のうち、前記特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞数は、時間的に大きな揺らぎを示し、かつ、
前記特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞複数においても、個々の細胞内における前記ｍＲＮＡ量を比較すると、細胞間で大きな揺らぎを示している生細胞群に対して、
その確率論的に推定可能な形態で表記するモデルとして、
前記生細胞群における、一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生し、
一つの細胞内における該特定の遺伝子の転写事象発生一回に付随して、前記ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化が生じ、
一つの細胞内において、一定の時間内における該特定の遺伝子の転写事象発生回数は、ポワソン分布に類する確率分布に従い、また、前記確率分布と反しない範囲で、前後する二回の転写事象発生間における時間間隔は、実質的にランダムであるという条件を満たすモデルを用いる
ことを特徴とする、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法である。
【発明の効果】
【００１７】
本発明にかかる確率論的モデル化の方法は、全体が均質な条件下に置かれている、同一種類の生細胞群を対象とする際、各細胞内における、特定の遺伝子の転写に起因して生成されるｍＲＮＡ量が時間的に大きな揺らぎを示し、同時に、各細胞内で転写に起因して生成されるｍＲＮＡ量は、細胞相互においても、大きな揺らぎを示している系において、かかる生細胞群全体における、特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の時間的変化、また、細胞相互における揺らぎに由来する分布を無理なく説明できるという効果がある。
【発明を実施するための最良の形態】
【００１８】
上記の問題点を解決するために、本発明では、図１により表記される転写産物反応モデルに基づき、個々の細胞において、個々のｍＲＮＡ量が時間的に、また、細胞相互間においても、大きな揺らぎを持って変化する事象を確率論的にモデル化している。
【００１９】
上述するように、本発明にかかる遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法は、
同種の細胞複数個を含んでなる生細胞群において、前記細胞が有する特定の遺伝子の転写に起因して生成されるｍＲＮＡに関して、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を確率論的に推定可能な形態で表記するモデル化を行う方法であって、
該生細胞群に含まれる同種の細胞複数個のうち、前記特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞数は、時間的に大きな揺らぎを示し、かつ、
前記特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞複数においても、個々の細胞内における前記ｍＲＮＡ量を比較すると、細胞間で大きな揺らぎを示している生細胞群に対して、
その確率論的に推定可能な形態で表記するモデルとして、
前記生細胞群における、一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生し、
一つの細胞内における該特定の遺伝子の転写事象発生一回に付随して、前記ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化が生じ、
一つの細胞内において、一定の時間内における該特定の遺伝子の転写事象発生回数は、ポワソン分布に類する確率分布に従い、また、前記確率分布と反しない範囲で、前後する二回の転写事象発生間における時間間隔は、実質的にランダムであるという条件を満たすモデルを用いる
ことを特徴とする、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法であるが、
この本発明にかかるモデル化の方法は、下記する形態とすることが可能である。
【００２０】
まず、本発明のかかるモデル化の方法は、
一つの細胞内における該特定の遺伝子の転写事象発生一回に付随して、前記ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化が生じ、
一つの細胞内において、一定の時間内における該特定の遺伝子の転写事象発生回数は、ポワソン分布に類する確率分布に従い、また、前記確率分布と反しない範囲で、前後する二回の転写事象発生間における時間間隔は、実質的にランダムである際に適用されるモデル化の方法であり、特には、
同種の細胞複数個を含んでなる生細胞群において、前記細胞が有する特定の遺伝子の転写に起因して生成されるｍＲＮＡに関して、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を確率論的に推定可能な形態で表記するモデル化を行う方法であって、
前記生細胞群における、一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生し、
個々の細胞内における前記ｍＲＮＡ量の時間的な変動は、
前記一定の時間を短い時間間隔（Δｔ）に分割し、該時間間隔（Δｔ）内における前記ｍＲＮＡ量（Ｘi）の変化量：ΔＸi、
転写に起因するｍＲＮＡの生成過程に相当する速度項：ｋgenerate＊Ｘi
（但し、ｋgenerateは、ある範囲内で時間的に変動するが、その時間的な平均値＜ｋgenerate＞は、特定の値に集束するような、時間的に変動する速度定数である）、
酵素的分解に起因するｍＲＮＡの消失過程に相当する速度項：−ｋdegrade＊Ｘi
（但し、ｋdegradeは、ある範囲内で時間的に変動するが、その時間的な平均値＜ｋdegrade＞は、特定の値に集束するような、時間的に変動する速度定数である）、
と定義する際、
前記ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、
前記ｍＲＮＡ量（Ｘi）に依存する項として、前記二つの速度項の和：ｋgenerate＊Ｘi−ｋdegrade＊Ｘiを含み、
前記時間的に変動する速度定数ｋgenerateとｋdegradeは、
前記時間的な平均値＜ｋgenerate＞と＜ｋdegrade＞とは、＜ｋdegrade＞≧＜ｋgenerate＞の条件を満足し、かつ互いに独立な確率函数に従って、時間的に変動する函数とし、
個々の細胞に対して、独立に選択される該速度定数ｋgenerateとｋdegradeを割り振ることにより、確率論的に推定可能な形態の表記がなされる
ことを特徴とする、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法として、具体化することが可能である。
【００２１】
その際、本発明の第一の形態にかかる確率論的モデル化の方法は、
個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、下記式（Ｉ）：
式（Ｉ）  ：  ΔＸi＝ｋgenerate＊Ｘi−ｋdegrade＊Ｘi
で示され、前記ｍＲＮＡ量（Ｘi）に依存する項である前記二つの速度項の和：ｋgenerate＊Ｘi−ｋdegrade＊Ｘiのみ含むものとする
ことを特徴とする確率論的モデル化の方法である。
【００２２】
また、本発明の第二の形態にかかる確率論的モデル化の方法は、
個々の細胞内における前記ｍＲＮＡ量（Ｘi）を該ｍＲＮＡ分子数と定義する際、該ｍＲＮＡ分子数の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、下記式（II）：
式（II）  ：  ΔＸi＝ｋgenerate＊Ｘi−ｋdegrade＊Ｘi＋Ｃ
（但し、Ｃは、ある範囲内で時間的に変動するが、その時間的な平均値＜Ｃ＞は、特定の値に集束するような、時間的に変動する定数項である）
かつ、
前記時間的な平均値＜ｋgenerate＞と＜ｋdegrade＞、ならびに＜Ｃ＞とは、
＜ｋdegrade＞≧＜ｋgenerate＞の条件を満足し、
前記時間的に変動する定数項Ｃは、
０≧＜ｋgenerate＞−＜ｋdegrade＞＋＜Ｃ＞の条件を満足するように、独立な確率函数に従って、時間的に変動する定数項とされている、
ことを特徴とする確率論的モデル化の方法である。
【００２３】
上記の本発明の第一の形態にかかる確率論的モデル化の方法においては、
個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、前記式（Ｉ）において、
前記時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝Ａ＋α＋Ｒ1
ｋdegrade＝Ａ＋Ｒ2
（但し、
Ａは、Ａ＞０  の条件を満たす定数、
αは、（Ａ＋α）＞０≧α  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝０であり、標準偏差σ（Ｒ1）は、（Ａ＋α）＞３σ（Ｒ1）の条件を満たす正規分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝０であり、標準偏差σ（Ｒ2）は、（Ａ＋α）＞３σ（Ｒ2）の条件を満たす正規分布に従う確率変数である）
を選択してなる、下記式（１）：
式（１）：  ΔＸi＝（Ａ＋α＋Ｒ1）＊Ｘi−（Ａ＋Ｒ2）＊Ｘi
で示される
ことを特徴とする確率論的モデル化の方法とすることが好ましい。
【００２４】
一方、上記の本発明の第二の形態にかかる確率論的モデル化の方法においては、
個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、前記式（II）において、
前記時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝α＋Ｒ1
ｋdegrade＝β＋Ｒ2
前記時間的に変動する定数Ｃとして、
Ｃ＝Ｒ3
（但し、
αは、α≧０  の条件を満たす定数、
βは、β≧０  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝λ1であり、ポワソン分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝λ2であり、ポワソン分布に従う確率変数、
Ｒ3は、平均値＜Ｒ3＞は、＜Ｒ3＞＝λ3であり、ポワソン分布に従う確率変数であって、
（β＋λ2）≧（α＋λ1）  ならびに  ０≧（α＋λ1）−（β＋λ2）＋λ3  の条件を満たす）
を選択してなる、下記式（２）：
式（２）：    ΔＸi＝（α＋Ｒ1）＊Ｘi−（β＋Ｒ2）＊Ｘi＋Ｒ3
で示される
ことを特徴とする確率論的モデル化の方法という態様を選択すると好ましい。
【００２５】
また、前記の態様をより汎用性を有するものとする上では、
上記の本発明の第二の形態にかかる確率論的モデル化の方法において、
個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する、ｍＲＮＡ量（Ｘi）の変化量：ΔＸiは、前記式（II）において、
前記時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝α1＋α2＊Ｒ1
ｋdegrade＝β1＋β2＊Ｒ2
前記時間的に変動する定数Ｃとして、
Ｃ＝Ｒ3
（但し、
α1は、α1≧０  の条件を満たす定数、α2は、α2＞０  の条件を満たす定数、
β1は、β1≧０  の条件を満たす定数、β2は、β2＞０  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝λ1であり、ポワソン分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝λ2であり、ポワソン分布に従う確率変数、
Ｒ3は、平均値＜Ｒ3＞は、＜Ｒ3＞＝λ3であり、ポワソン分布に従う確率変数であって、
（β1＋β2＊λ2）≧（α1＋α2＊λ1）  ならびに  ０≧（α1＋α2＊λ1）−（β1＋β2＊λ2）＋λ3  の条件を満たす）
を選択してなる、下記式（３）：
式（３）：    ΔＸi＝（α1＋α2＊Ｒ1）＊Ｘi−（β1＋β2＊Ｒ2）＊Ｘi＋Ｒ3
で示される
ことを特徴とする確率論的モデル化の方法という態様を選択することも可能である。
【００２６】
具体的には、上述する本発明の第一の形態ならびに第二の形態にかかる確率論的モデル化の方法に基づき、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象を反映する、前記ｍＲＮＡ量の時間的変動の確率論的なシミュレーションを実施する際には、上述するモデル函数系に従って、数値演算を行う。従って、本発明は、かかるモデル函数系に従う数値演算の実施に利用される演算処理用制御プログラムの発明をも提供する。すなわち、本発明にかかる演算処理用制御プログラムは、
上記の種々の構成を有する、本発明の第一の形態あるいは第二の形態にかかる遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法に基づき、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象を反映する、前記ｍＲＮＡ量の時間的変動の確率論的なシミュレーションを、コンピュータ上において実現するための演算処理用制御プログラムである。かかる演算処理用制御プログラムは、演算処理を実施するコンピュータにより読み取り可能な形態で、該演算処理用制御プログラムを格納した記録媒体の態様とすることも可能である。すなわち、本発明にかかる演算処理用制御プログラムを格納した記録媒体は、
上記の種々の構成を有する、本発明の第一の形態あるいは第二の形態にかかる遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法に基づき、遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象を反映する、前記ｍＲＮＡ量の時間的変動の確率論的なシミュレーションを、コンピュータ上において実現するための演算処理用制御プログラムを、前記コンピュータにより読み取り可能な形態で格納した記録媒体である。
【００２７】
加えて、本発明は、生細胞群における、特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の定量方法の発明をも提供しており、すなわち、本発明にかかるｍＲＮＡ量の定量方法は、
同種の細胞複数個を含んでなる生細胞群において、前記細胞が有する特定の遺伝子の転写に起因して生成されるｍＲＮＡに関して、細胞あたりのｍＲＮＡ量平均を、該生細胞群から採取されるサンプル群を用いて、定量する方法であって、
前記生細胞群を構成する個々の細胞中における、前記特定の遺伝子に由来するｍＲＮＡ量は、時間的に揺らぎを示し、また、個々の細胞間においても、前記特定の遺伝子の発現量に揺らぎを示す際、
個々の細胞中における、前記特定の遺伝子に由来するｍＲＮＡ量の時間的な揺らぎの分布、あるいは、個々の細胞間における、前記特定の遺伝子の発現量の揺らぎの分布の大きさに基づき決定される、基準の細胞数に対して、該生細胞群から採取されるサンプル群を構成する細胞数を前記基準の細胞数より多くするサンプリング方法を用い、
前記細胞数を有するサンプル群において測定される、細胞あたりのｍＲＮＡ量平均を、該生細胞群全体における細胞あたりのｍＲＮＡ量平均とする
ことを特徴とする、生細胞群における、特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の定量方法である。
【００２８】
その際、本発明にかかるｍＲＮＡ量の定量方法においては、
個々の細胞中における、前記特定の遺伝子に由来するｍＲＮＡ量の時間的な揺らぎの分布、あるいは、個々の細胞間における、前記特定の遺伝子の発現量の揺らぎの分布の大きさを、請求項１に記載の遺伝子の転写に起因する細胞内ｍＲＮＡ量変動現象の確率論的モデル化の方法を利用して、推定する
ことを特徴とするｍＲＮＡ量の定量方法とすることが好ましい。
【００２９】
特には、
所定の時間内において、前記生細胞群における、前記特定の遺伝子の転写がなされている細胞数の比率の多寡に従って、前記基準の細胞数を決定することが好ましい。具体的には、所定の時間内において、前記生細胞群における、前記特定の遺伝子の転写がなされている細胞数の比率に反比例するように、前記基準の細胞数を決定することがより好ましい。
【００３０】
例えば、前記の基準の細胞数を決定する方法を選択する場合、
生細胞群における、前記特定の遺伝子の転写がなされている細胞数の比率に応じて、
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１０，０００個以上である；
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１００，０００個以上である；
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１，０００，０００個以上である；
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１０，０００，０００個以上である；
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１００，０００，０００個以上である；
前記生細胞群から採取されるサンプル群を構成する細胞数の選択に用いる、前記基準の細胞数は、少なくとも１，０００，０００，０００個以上である；
様々なサンプリング条件を選択することが望ましい。
以下に、本発明をより詳しく、具体的な実施形態を参照しつつ説明する。
【００３１】
細胞内における特定の遺伝子の発現過程では、先ず、該遺伝子（ＤＮＡ）から転写によって、対応するＲＮＡ鎖が生成される。次いで、遺伝子（ＤＮＡ）の塩基配列が、アミノ酸配列のコードを行うエキソン領域と、その間に存在するイントロン配列を含む場合には、スプライシング過程を経て、ｍＲＮＡの生成がなされる。さらに、ｍＲＮＡに必要に応じて、様々な化学的修飾がなされた後、リボソームにおいて、ｍＲＮＡの塩基配列に従って、対応するアミノ酸配列を有するペプチド鎖へと翻訳される。翻訳されたペプチド鎖は、適正にホールディングされ、同時に、必要に応じて、翻訳後修飾が施され、最終的に、所定の三次元構造を有するタンパク質分子となる。所謂、セントラルドグマは、このＤＮＡ上の遺伝情報に基づき、対応するアミノ酸配列を有するタンパク質の生合成がなされる過程を表現している。すなわち、ｍＲＮＡの生成量は、遺伝子（ＤＮＡ）から転写によって、生成されるＲＮＡ鎖の量に比例している。
【００３２】
例えば、アルターネイティブ・スプライシングの過程は、対応する完全長のＲＮＡ鎖から、その間に存在するイントロン配列を除去するスプライシング過程において、異なるスプライシングが起こり、一部のエキソン領域も併せて除去される。その結果、生成されるｍＲＮＡでは、完全長のＲＮＡ鎖上には存在していた、一部のエキソン領域が欠落し、対応する部分アミノ酸配列を欠く、スプライシング変異体の翻訳に利用可能なものとなる。このように、対応する完全長のＲＮＡ鎖から、複数種のｍＲＮＡ分子が生成される場合であっても、ｍＲＮＡの生成量は、原理的には、遺伝子（ＤＮＡ）から転写によって、生成されるＲＮＡ鎖の量に比例している。
【００３３】
その後、生成されたｍＲＮＡの塩基配列に基づき、リボソームにおいて、対応するアミノ酸配列を有するペプチド鎖へと翻訳される。この翻訳過程に関与するリボソームは、ｍＲＮＡの種類を区別せず、対応するペプチド鎖への翻訳を行うため、生成されたｍＲＮＡ量と、合成されるペプチド鎖の量とは、原則的に比例している。リボソームにおける翻訳に利用されたｍＲＮＡ分子は、その後、離脱して、細胞内に一時期残留することになる。
【００３４】
一方、種々の研究結果から、遺伝子（ＤＮＡ）から転写過程に伴って、生成されたｍＲＮＡは、一時期細胞内に蓄積されるが、その後、消失し、最終的には、転写過程が開始する前と同様の極めて低い水準となることも知られている。このｍＲＮＡの消失過程は、ＲＮＡ分解酵素（ＲＮａｓｅ）が関与する酵素的分解過程に起因している。この酵素的分解過程は、ＲＮＡ分解酵素（ＲＮａｓｅ）と、その基質であるｍＲＮＡとの一次反応であり、その反応速度は、基質であるｍＲＮＡ量の一次に比例している。
【００３５】
一方、遺伝子（ＤＮＡ）から転写によって、対応するＲＮＡ鎖が生成される過程は、ＲＮＡ鎖の生成反応自体は、ＲＮＡポリメラーゼによる酵素反応であり、その際、鋳型となる遺伝子（ＤＮＡ）は、原理的に、一細胞中の存在量は、同種の細胞では等しい。この転写過程は、例えば、転写因子タンパク質による転写促進など、他の促進因子によって誘導される過程であり、転写促進因子の濃度に依存している。従って、細胞内において、転写促進因子の消費が進み、その濃度の枯渇が生じると、転写の起こる速度が急速に低下する。従って、細胞内において、転写促進因子の濃度の上昇とともに、急速に生成されたｍＲＮＡの量は増加するが、ある時間を経過すると、増加量（速度）が急速に低下する。
【００３６】
細胞内に存在するｍＲＮＡ量の時間的な変化は、この転写過程に起因する、ｍＲＮＡの生成過程と、上記のｍＲＮＡの酵素的分解過程に因るｍＲＮＡの分解過程とで決定されている。通常、その過渡的な変化は、転写過程が生じる前、ｍＲＮＡ量は極めて低い水準であるが、転写過程が誘起されると、急速に上昇し、ある時点でピークを向かえ、その後、また急速に低下し、当初の極めて低い水準に復する。すなわち、ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化を示すことになる。
【００３７】
本発明が対象とする系は、全体が均質な条件下に置かれている、同一種類の生細胞群であるが、その生細胞群中に含まれる各細胞内における、特定の遺伝子の転写に起因して生成されるｍＲＮＡ量に着目すると、各細胞内におけるｍＲＮＡ量は時間的に大きな揺らぎを示すものである。すなわち、各細胞内におけるｍＲＮＡ量（Ｘi）の時間的な変化を追跡すると、多くの時間は、実質的にｍＲＮＡ量は「０」となっているが、時折、前記遺伝子（ＤＮＡ）から転写が起こり、シャープなピーク状のｍＲＮＡ量の時間的な変化を示す系である。同時に、各細胞内で転写に起因して生成されるｍＲＮＡ量は、細胞相互においても、大きな揺らぎを示している系である。このような挙動を示す系は、生細胞群中に含まれる同種の細胞複数個のうち、ある時点では、前記特定の遺伝子の転写に起因して生成されるｍＲＮＡ量の変動を示している細胞は幾つか存在しているが、その細胞数は、時間的に大きな揺らぎを示し、また、この幾つかの細胞では、前記のシャープなピーク状のｍＲＮＡ量の時間的な変化が起こっているが、そのピークにおけるｍＲＮＡ量を比較すると、細胞間で大きな揺らぎを示している系と見做すことが可能である。
【００３８】
従って、本発明においては、前記の系において、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を確率論的に推定可能な形態で表記するモデルとして、
（ａ）  一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生する。
（ｂ）  一つの細胞内における該特定の遺伝子の転写事象発生一回に付随して、ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化が生じている。
（ｃ）  一つの細胞内において、一定の時間内における該特定の遺伝子の転写事象発生回数は、散発的に発生する事象の出現頻度の記述に適する、ポワソン分布に類する確率分布に従い、また、この確率分布と反しない範囲で、前後する二回の転写事象発生間における時間間隔は、実質的にランダムとなっている。
この３条件を満足するモデルを採用している。
【００３９】
先ず、第一の条件（ａ）を満足する上では、各細胞に関して、その時間的変動を表記する函数の表現形式は、同一であるが、選択される確率変数自体は、個々の細胞で独立に選択することで達成できる。
【００４０】
次に、第二の条件（ｂ）と第三の条件（ｃ）を満たす手法として、本発明では、下記の手法を利用している。
【００４１】
先ず、ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化は、時間（ｔ）を極く短い時間間隔（Δｔ）に分割し、この時間間隔（Δｔ）内におけるｍＲＮＡ量（Ｘi）の変化量ΔＸiを考慮すると、この変化量ΔＸiは、該特定の遺伝子の転写事象発生に相当する、短い時間では、ΔＸiが「０」よりも有意に大きくなるが、その後は、酵素的分解反応に伴って、ΔＸiは「０」よりも有意に小さくなるような、確率函数を用いて表現可能である。
【００４２】
すなわち、スパイク状のｍＲＮＡ量変化を生じさせるメカニズムは、細胞内において、遺伝子発現促進に付随して、転写の促進が起こると、対応するｍＲＮＡの量（Ｘi）が急速に上昇し、その後、転写抑制がなされると、ｍＲＮＡ量の増加は停止し、一方、細胞内に存在するＲＮＡ分子の酵素的分解反応に伴って、ｍＲＮＡの残量は指数関数的に減衰していくことに相当すると仮定した。従って、この過渡的プロセスを極く短い時間間隔（Δｔ）に分割し、この時間間隔（Δｔ）内におけるｍＲＮＡ量（Ｘi）の変化量ΔＸiを考慮すると、この変化量ΔＸiは、転写過程に起因するｍＲＮＡの生成過程と、ＲＮＡ分子の酵素的分解反応に伴う、ｍＲＮＡの分解過程を反映する、二つの主要な成分（項）を含む形式で表現すべきことに想到した。
【００４３】
その際、ＲＮＡ分子の酵素的分解反応の速度は、この酵素的分解反応に関与するＲＮＡ分解酵素タンパク質の濃度、また、分解を受けるＲＮＡ分子の濃度に依存する。一般に、分解を受けるＲＮＡ分子は、対象となるｍＲＮＡ以外にも、複数種のＲＮＡ分子が存在し、その間で競争が起こっていることが想定される。また、ＲＮＡ分解酵素タンパク質（ＲＮａｓｅ）の濃度（ＣRNase）自体も、ある範囲内で経時的に変動していることが予測される。この二つの変動要素を考慮すると、かかる酵素的分解反応に因る、ｍＲＮＡの分解過程の速度；−ｋdegrade＊Ｘiは、ある範囲内で時間的に変動するが、その時間的な平均値は、特定の値に集束するような、時間的に変動する速度定数ｋdegradeを含む形式で表記すべきことが判った。すなわち、この時間的に変動する速度定数ｋdegradeは、時間的な平均値＜ｋdegrade＞は、特定の値を示し、この平均値を中心として、なんらかの確率に従って変動する可変的な函数として、表現することが可能であることを見出した。
【００４４】
一方、転写過程に起因するｍＲＮＡの生成過程に関しては、スパイク状のｍＲＮＡ量変化の前段部分では、ｍＲＮＡ量の増加が急速に進行する形態を反映するような生成の速度として、ｋgenerate＊Ｘiの形式を利用することが可能であることを見出した。この速度定数ｋgenerateも、ある範囲内で時間的に変動するが、その時間的な平均値は、特定の値に集束するような、時間的に変動する速度定数である。すなわち、速度定数ｋgenerateに関しても、時間的な平均値＜ｋgenerate＞は、特定の値を示し、この平均値を中心として、なんらかの確率に従って変動する可変的な函数として、表現することが可能であることを見出した。
【００４５】
ｍＲＮＡ量（Ｘi）の変化量ΔＸiを表現する函数中、ｍＲＮＡ量（Ｘi）に依存する成分として、（ｋgenerate＊Ｘi−ｋdegrade＊Ｘi）を含む形式とする。スパイク状のｍＲＮＡ量変化が生じていない、大半の時間では、（ｋgenerate−ｋdegrade）＊Ｘiは、「０」または、「０」以下（負の値）とすると、ｍＲＮＡ量（Ｘi）の急激な増加が生じてない状態を表現できる。すなわち、ｍＲＮＡ量（Ｘi）には、当然にＸi≧０の条件を付しているので、ｋgenerateとｋdegradeは、それぞれ独立に、なんらかの確率に従って変動するが、大半の時間では、（ｋgenerate−ｋdegrade）は、「０」または、「０」以下（負の値）となるように、確率的に変動する可変的な函数を選択する。
【００４６】
例えば、時間的な平均値＜ｋdegrade＞と時間的な平均値＜ｋgenerate＞とを比較した際、＜ｋdegrade＞≧＜ｋgenerate＞の条件を満足すると、大半の時間では、（ｋgenerate−ｋdegrade）は、「０」または、「０」以下（負の値）となるという条件を満たす。
【００４７】
図１は、先に説明した、本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化に対応する「ｍＲＮＡの生成・分解」過程のモデル化の概念を模式化して示したものである。ＲＮＡ分子の酵素的分解反応に伴う、ｍＲＮＡの分解過程１０２と、転写過程に起因するｍＲＮＡの生成過程１０１とを、時間間隔（Δｔ）内におけるｍＲＮＡ量（Ｘi）の変化量ΔＸiを支配する主要な要素とする点では、セントラルドグマで記述される、タンパク質の生合成における「ｍＲＮＡの生成」過程と、その後の酵素的分解反応に伴う「ｍＲＮＡの分解」過程を考慮するという観点で、従来の手法と共通点を有している。
【００４８】
一方、転写過程に起因するｍＲＮＡの生成過程１０１は、実際には、細胞内に残存しているｍＲＮＡ量（Ｘi）に依存するものではないが、本発明では、ｍＲＮＡ量の増加が急速に進行する形態を反映するような生成の速度として、ｋgenerate＊Ｘiの形式を利用することで、数値演算に適合するモデル化を行っている。細胞内に残存しているｍＲＮＡ量（Ｘi）に依存する、ｍＲＮＡの生成が生じる可能性を有する生体内反応としては、例えば、「ｍＲＮＡによるＲＮＡポリメラーゼの活性化」のような自己触媒的な反応を仮定する必要がある。但し、ｍＲＮＡ分子の種類に依らず、この種の自己触媒的な反応を仮定することは困難である。しかしながら、ｍＲＮＡ量の増加が急速に進行する過程と、その後のｍＲＮＡ量の減少過程とを、簡便に数値演算可能な形式で、近似的に表現する手法として、前述のモデル化を本発明では利用している。
【００４９】
すなわち、本発明において、ｍＲＮＡ量（Ｘi）の変化量ΔＸiを表現する函数中、ｍＲＮＡ量（Ｘi）に依存する成分として、（ｋgenerate＊Ｘi−ｋdegrade＊Ｘi）を含む形式とすることの技術的な意義は、互いに独立した二つの確率に従って変動する、確率変数ｋgenerateとｋdegradeの差違（ｋgenerate−ｋdegrade）が「０」を超える頻度は、限られたものとなり、特に、「０」を超える時間間隔（Δｔ）が複数、一時期に集中する確率は更に限られたものとなり、「離散」的にシャープなピーク状のｍＲＮＡ量の時間的な変化が生じる現象と、形態的に類似する状態を表現可能な「モデル」となる点である。
【００５０】
時間的な平均値＜ｋdegrade＞と＜ｋgenerate＞に関して、＜ｋdegrade＞≧＜ｋgenerate＞の条件に加えて、「負」の範囲で、差違｛＜ｋdegrade＞−＜ｋgenerate＞｝の絶対値が増すとともに、差違（ｋgenerate−ｋdegrade）が「０」を超える頻度は、益々低下していく。その際、「０」を超える時間間隔（Δｔ）が複数、一時期に集中する確率は、極く限られたものとなり、「散発」的にシャープなピーク状のｍＲＮＡ量の時間的な変化が生じる現象と、形態的に類似する状態を表現可能な「モデル」となる。
【００５１】
確率変数ｋgenerateとｋdegradeとは、互いに独立の確率に従うため、差違（ｋgenerate−ｋdegrade）が「０」を超える頻度は、その「独立の確率」の選択の相違に起因して、原理的には、大きな揺らぎを示すものとなり、「０」を超える時間間隔（Δｔ）が複数、一時期に集中する時点の間隔、すなわち、「離散」的ピーク間の時間間隔も実質的にランダムとなる。あるいは、「０」を超える時間間隔（Δｔ）が複数、一時期に集中する際、その集中が生じている期間の長さも、大きな揺らぎを示すものとなり、すなわち、「離散」的ピークのピーク高さも、大きな揺らぎを示す。
【００５２】
以下に、本発明にかかるモデル化方法に基づく、実施の形態の代表的な例を示し、本発明が提案する手法をより詳しく説明する。
【００５３】
本発明にかかるモデル化方法では、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を表記する際、時間（ｔ）を極く短い時間間隔（Δｔ）に分割し、この時間間隔（Δｔ）内におけるｍＲＮＡ量（Ｘi）の変化量ΔＸiを利用して表記する。
【００５４】
その際、本発明の第一の形態にかかるモデル化の方法では、
個々の細胞内における、前記ｍＲＮＡ量（Ｘi）の変化量：ΔＸiを、下記式（Ｉ）：
式（Ｉ）  ：  ΔＸi＝ｋgenerate＊Ｘi−ｋdegrade＊Ｘi
で示し、ｍＲＮＡ量（Ｘi）に依存する項である前記二つの速度項の和：ｋgenerate＊Ｘi−ｋdegrade＊Ｘiのみ含むものとする。
【００５５】
つまり、前記式（Ｉ）に基づき、時刻ｔ＝０における初期値Ｘi（０）から、
Ｘi（ｔ＋Δｔ）＝Ｘi（ｔ）＋ΔＸi
＝Ｘi（ｔ）＋（ｋgenerate−ｋdegrade）＊Ｘi（ｔ）
と、逐次的に数値計算を進めることで、ｍＲＮＡ量（Ｘi）の時間的な挙動が表記される。
【００５６】
この表記法では、ある時点で、Ｘi（ｔ）が「０」となると、その後、（ｋgenerate−ｋdegrade）が変化しても、ｍＲＮＡ量（Ｘi）は、「０」に留まったままとなるため、ｍＲＮＡ量（Ｘi）は、Ｘi（ｔ）＞０の条件を満たすという付加条件を設ける。実際には、上記の式に基づき数値計算により、Ｘi（ｔ）を算出していく過程で、Ｘi（ｔ）≦０となった時点で、「０」と見なせる「最低値」へと置き換える。また、初期値Ｘi（０）も、通常、この「最低値」を設定する。
【００５７】
また、本発明の第二の形態にかかるモデル化の方法では、個々の細胞内における前記ｍＲＮＡ量（Ｘi）を該ｍＲＮＡ分子数と定義する際、つまり、Ｘiは、自然数；１、２、…の離散的な値とする場合、
個々の細胞内における、前記ｍＲＮＡ量（Ｘi）の変化量：ΔＸiを、下記式（II）：
式（II）  ：  ΔＸi＝ｋgenerate＊Ｘi−ｋdegrade＊Ｘi＋Ｃ
で示し、ｍＲＮＡ量（Ｘi）に依存する項である前記二つの速度項の和：ｋgenerate＊Ｘi−ｋdegrade＊Ｘiと、微小なノイズ項に相当する、時間的に変動する定数項とを含むものとする。この微小なノイズ項に相当するＣは、ある範囲内で時間的に変動するが、その時間的な平均値＜Ｃ＞は、特定の値に集束するような、時間的に変動する定数項である。
【００５８】
一方、ｍＲＮＡ量（Ｘi）は、前記式（Ｉ）に基づき、時刻ｔ＝０における初期値Ｘi（０）から、
Ｘi（ｔ＋Δｔ）＝Ｘi（ｔ）＋ΔＸi
＝Ｘi（ｔ）＋（ｋgenerate−ｋdegrade）＊Ｘi（ｔ）＋Ｃ
と、逐次的に数値計算を進めることで、ｍＲＮＡ量（Ｘi）の時間的な挙動が表記される。
【００５９】
この表記法でも、上記の式に基づき数値計算により、Ｘi（ｔ）を算出していく過程で、Ｘi（ｔ）≦０となることもある。その際、Ｘi（ｔ）≧０の条件を満たすように、「０」と見なせる「最低値」へと置き換える。また、初期値Ｘi（０）も、通常、この「最低値」を設定する。Ｘiは、自然数と定義しており、「最低値」として、１を選択する。
【００６０】
大半の時間は、ｍＲＮＡ量（Ｘi）は、前記「最低値」、またはその近辺に留まるという要件は、前記時間的な平均値＜ｋgenerate＞と＜ｋdegrade＞、ならびに＜Ｃ＞とを、
＜ｋdegrade＞≧＜ｋgenerate＞の条件を満足し、
一方、時間的に変動する定数項Ｃは、
０≧＜ｋgenerate＞−＜ｋdegrade＞＋＜Ｃ＞の条件を満足するように、独立な確率函数に従って、時間的に変動する定数項とすることで達成される。
【００６１】
すなわち、大半の時間では、０≧（ｋgenerate−ｋdegrade）となっており、また、Ｘi（ｔ）が、「最低値」の１である状態では、ΔＸi＝（ｋgenerate−ｋdegrade）＋Ｃは、「０」または、「０」以下となる条件も満足されている。
【００６２】
（第一の実施態様）
本発明の第一の形態にかかるモデル化の方法に基づく、一つの実施態様を説明する。
【００６３】
個々の細胞内におけるｍＲＮＡ量（Ｘi）の変化量：ΔＸiを前記式（Ｉ）により表記する際、時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝Ａ＋α＋Ｒ1
ｋdegrade＝Ａ＋Ｒ2
（但し、
Ａは、Ａ＞０  の条件を満たす定数、
αは、（Ａ＋α）＞０≧α  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝０であり、標準偏差σ（Ｒ1）は、（Ａ＋α）＞３σ（Ｒ1）の条件を満たす正規分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝０であり、標準偏差σ（Ｒ2）は、（Ａ＋α）＞３σ（Ｒ2）の条件を満たす正規分布に従う確率変数である）
を選択する。その際、前記式（Ｉ）は、下記式（１）：
式（１）：  ΔＸi＝（Ａ＋α＋Ｒ1）＊Ｘi−（Ａ＋Ｒ2）＊Ｘi
で示される。
【００６４】
この式（１）は、
ΔＸi＝（α＋Ｒ1−Ｒ2）＊Ｘi
と変形でき、二つの確率変数Ｒ1とＲ2と、閾値パラメータ（−α）とで規定される函数となっている。二つの確率変数Ｒ1とＲ2の差違（Ｒ1−Ｒ2）は、平均値は０であるが、単一の確率変数Ｒ（Ｒは、平均値＜Ｒ＞は、＜Ｒ＞＝０であり、標準偏差σ（Ｒ）は、（Ａ＋α）＞３σ（Ｒ）の条件を満たす正規分布に従う確率変数）を用いる場合と比較すると、通常、より大きな揺らぎを示す。差違（Ｒ1−Ｒ2）が、閾値パラメータ（−α）を超えると、ΔＸi≧０となり、ｍＲＮＡ量（Ｘi）の増加につながるが、この閾値パラメータ（−α）を超える頻度の揺らぎは、二つの確率変数Ｒ1とＲ2の標準偏差σ（Ｒ1）と標準偏差σ（Ｒ2）の選択に依存して、大きな揺らぎを示すものとなる。
【００６５】
表現すべき個々の細胞内におけるｍＲＮＡ量（Ｘi）の変化は、スパイク状のピークを示す頻度、ならびに、各スパイク状のピークのピーク値の双方とも、大きな揺らぎを示すものであり、少なくとも「独立」した確率変数２を含む式（１）の表記は、かかる揺らぎ現象のモデル化に適合している。
【００６６】
実際のＲＮＡ分解酵素（ＲＮａｓｅ）が関与する酵素的分解反応に因る、ｍＲＮＡの分解過程の速度；−ｋdecomposition＊ＣRNase＊Ｘiを規定する速度定数：ｋdecompositionは、例えば、酵素的分解過程のみが進行する際には、（ΔＸi／Δｔ）＝−ｋdecomposition＊ＣRNase＊Ｘiのような指数関数的な減少を記述する。外見上は、（ｋdecomposition＊ＣRNase）＊Δｔと、ｋdegrade＝Ａ＋Ｒ2とが対応するように思えるが、上記式（１）に示すモデルでは、定数Ａは、最終的に式中には含まれないものとなり、実際には、両者間に対応関係は無いものである。
【００６７】
すなわち、定数Ａは、任意に選択することが可能である。従って、上記式（１）の一形態として、
式（１’）：  ΔＸi＝（Ａ＋α＋Ｒ1）＊Ｘi−（Ａ＋Ｒ2）＊Ｘi
（但し、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝０であり、標準偏差σ（Ｒ1）は、σ（Ｒ1）＝１の条件を満たす正規分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝０であり、標準偏差σ（Ｒ2）は、σ（Ｒ2）＝１の条件を満たす正規分布に従う確率変数、
Ａは、Ａ＞３  の条件を満たす定数、
αは、（Ａ＋α）＞０≧αの条件を満たす定数である）
の表記を利用することができる。
【００６８】
前述する、本発明にかかるモデル化の方法に対応する、式（１’）の表現形式のモデル函数を利用して、個々の細胞内におけるｍＲＮＡ量（Ｘi）の時間推移をシミュレーションする際には、数値演算により行われる。
【００６９】
例えば、式（１’）の表現形式のモデル函数として、シミュレーションする時間を、極く短い時間間隔（Δｔ）×Ｎ区分に分割し、時間（ｔ）を、ｔk＝ｋ＊Δｔ（ｋ＝０、１、２、…、Ｎ）の数列とし、各細胞（ｑ）に対して、時刻（ｔk）における確率変数Ｒ1q（ｔk）、Ｒ2q（ｔk）を設定する。任意定数Ａ（Ａ＞３）と、閾値パラメータ（−α）とは、生細胞群が均一と仮定しているので、各細胞（ｑ）に対して、同じ値を選択する。
【００７０】
従って、各細胞について、その細胞内におけるｍＲＮＡ量（Ｘiq）の変化量：ΔＸiqを表記するモデル函数として、
ΔＸiq(tk)＝（Ａ＋α＋Ｒ1q(tk)）＊Ｘiq(tk)−（Ａ＋Ｒ2q(tk)）＊Ｘiq(tk)
が設定される。
【００７１】
なお、Ｘiqの最低値：Ｘiq0は、（Ｘiq0×Δｔ×Ｎ）の値が、シミュレーションされたｍＲＮＡ量（Ｘiq）の時間変化を、この時間幅（Δｔ×Ｎ）に関して数値積分した際、離散的に出現しているスパイク状ピーク部分の積分値合計の１％を超えない範囲に選択する。また、シミュレーションの際、初期値として採用するＸiq(t0)は、通常、この最低値：Ｘiq0を選択する。すなわち、時刻（ｔ0）では、生細胞群の全ての細胞では、対象とするｍＲＮＡ量は、実質的に「０」の状態とする。
【００７２】
数値演算によるシミュレーションでは、図２に例示するような情報処理装置が利用される。基本的には、外部記憶装置２０１、中央処理装置（ＣＰＵ）２０２、メモリ２０３、入出力装置２０４からなる装置構成とされる。外部記憶装置２０１は、上記のモデル函数の生成に必要な確率変数（Ｒ1、Ｒ2）の発生、そのモデル函数に基づく数値演算などの実行に利用する各種プログラムや、シミュレーションの結果を保持する。中央処理装置（ＣＰＵ）２０２は、主に、前記の各種プログラムの実行、周辺装置を含む、システム全体の制御を行う。メモリ２０３は、中央処理装置（ＣＰＵ）２０２が使用する実行プログラム、サブルーチン、あるいは、データを一時的に記録する。入出力装置２０４は、ユーザーとのインタラクションを行う。多くの場合、本発明に基づく、モデル函数を利用するシミュレーションなどを実施する際、各種のプログラム実行のトリガーは、この入出力装置を介してユーザーが出す。また、ユーザーが、シミュレーション結果を見たり、上述の各種のパラメータ、定数の入力をしたり、プログラムの実行制御などを、この入出力装置を介して行う。
【００７３】
本発明にかかるモデル化の方法に基づく、モデル函数系を適用した数値演算によるシミュレーションは、例えば、図２８に示すフローに従って実施することが可能である。シミュレーション工程は、シミュレーション開始ステップ、ｍＲＮＡ量（Ｘi）の時間推移を逐次的に数値演算する計算ループのステップ、シミュレーション終了ステップの三段階のステップで構成される。
【００７４】
先ず、シミュレーション開始ステップは、シミュレーション条件を設定した後、計算ループへと移行させる工程である。すなわち、適用するモデル函数系を定義する各種パラメータ、逐次的な数値演算の結果得られる、各時刻におけるｍＲＮＡ量（Ｘi）の発散を回避する目的で設定される、上下の境界条件（境界値）、適用するモデル函数系に用いる、複数の確率変数に関する条件設定（対応する確率函数種類の選択）、複数の確率変数それぞれの生成に利用する乱数系の選択（乱数系の種類と、初期値の設定）、計算ループのループ段数の設定（シミュレーションする時間幅の選択）など、個々のシミュレーション計算に必要な、外部から入力される条件パラメータ・セットを用意し、計算ループを開始する工程である。
【００７５】
一方、シミュレーション終了ステップは、計算ループにおける実施するループ段数が完了した後、計算ループを終了させ、そのシミュレーション結果を、そのシミュレーション計算において用いた、前記条件パラメータ・セットとともに、外部へ出力（転送・記録）する工程である。
【００７６】
シミュレーション開始ステップにおける、計算ループへの移行は、ユーザーのトリガーにより行う。一方、シミュレーション終了ステップにおける、計算ループの終了は、通常、予め設定したループ段数に到達したか、否かという条件判断が、計算ループに設けられており、その条件判断に従って、自動的に実施される。なお、前記の実行済ループ段数の「カウント・アップ」に基づく条件判断に従った、自動的終了に加えて、予め設定したループ段数に到達していない段階でも、ユーザーのトリガーにより、強制的に計算ループの終了を行わせる「強制終了」モードも、通常設けられている。この「強制終了」モードでは、例えば、その時点で進行中のループを中断し、その時点で既に完了しているループ段数までのシミュレーション結果と、そのシミュレーション計算において用いた、前記条件パラメータ・セットとが、外部へ出力（転送・記録）される形態とすることができる。
【００７７】
計算ループでは、各ループ段において、各時刻のｍＲＮＡ量（Ｘi）を逐次的に数値演算し、予め設定したループ段数に到達した時点で、所定の時間幅におけるｍＲＮＡ量（Ｘi）の時間推移のシミュレーション結果が得られる。この逐次的な数値演算過程では、上述するように、シミュレーション条件を規定するパラメータ・セットに基づき、演算操作が実行される。なお、その演算操作のプログラムは、通常、シミュレーション条件を規定するパラメータ・セットに関して、デフォルト設定のパラメータ値を含む形式で作成されており、初期ステップにおいて、ユーザーが外部から入力する、特定の条件パラメータ・セットで、デフォルト設定のパラメータ値を置き換えた上で、その後の一連のループが実行される。その際、特定の条件パラメータ・セットのうち、デフォルト設定のパラメータ値と同じ値を選択する際には、その部分に関しては、ユーザーが外部から入力する、特定の条件パラメータ・セットにおいて、割愛し、デフォルト設定のパラメータ値を利用することが可能な形態とされる。
【００７８】
計算ループの各ループは、乱数発生、四則演算、境界値設定の３つのステップで構成される。乱数発生ステップでは、シミュレーションにおいて採用するモデル函数系において定義されている、各確率変数それぞれについて、各ループで用いる値の割り付けに利用する乱数の発生がなされる。すなわち、各確率変数は、対応する確率函数種類、例えば、正規分布やポワソン分布などの確率分布に基づき、作成される「数値集合」から、発生された乱数に基づき、選択される。
【００７９】
四則演算ステップでは、先ず、発生された乱数に基づき、各確率変数それぞれについて、そのループで用いる値の割り付けを行い、次いで、採用されているモデル函数系に従って、前ループ段（直前の時刻）のｍＲＮＡ量（Ｘi）から、そのループ段（時刻）のｍＲＮＡ量（Ｘi）を算出する。
【００８０】
境界値設定ステップでは、各時刻におけるｍＲＮＡ量（Ｘi）の発散を回避する目的で設定される、上下の境界条件（境界値）に基づき、算出されたｍＲＮＡ量（Ｘi）が前記上下の境界条件を満たすか否かを判定する。仮に、算出されたｍＲＮＡ量（Ｘi）が上下の境界条件を逸脱している場合には、上下の境界値で置き換え、そのループ段（時刻）のｍＲＮＡ量（Ｘi）とする。
【００８１】
計算ループでは、以上の３つのステップが、予め設定したループ段数繰り返される。
【００８２】
図３に、上記式（１’）に示すモデル函数系を用いたシミュレーション結果の一例を示す。図３に示す例では、シミュレーションする時間幅（Δｔ）×Ｎは、Ｎ＝１０００、Δｔ＝０．１秒、任意定数Ａは、Ａ＝１０、閾値パラメータ（−α）は、α＝−０．０２５をそれぞれ設定している。Ｘiの最低値：Ｘi0は、Ｘi0＝１．０＊１０-10を設定し、初期値：Ｘi(t0)として、Ｘi(t0)＝１．０＊１０-10としている。なお、数値演算の過程で、Ｘi(t)が、Ｘi(t)＜１．０＊１０-10となった場合、最低値Ｘi0：１．０＊１０-10の値に置き換えるという境界条件を選択している。
【００８３】
図３に示すシミュレーション結果の一例において、離散的なスパイク状のピークが示され、その時間間隔には、規則性は見出されない。また、各ピークの高さも大きな揺らぎを示している。従って、上記式（１’）に示すモデル函数系は、本発明が目標とするモデル化の条件を満足する結果を与えていると判断される。なお、各ピークの高さの水準、１．０＊１０-8のオーダーと比較し、最低値Ｘi0＝１．０＊１０-10は、１％以下となっており、この境界条件は、概ね妥当な値であることも確認される。
【００８４】
上記式（１’）に示すモデル函数系において、図３に示すシミュレーションに用いた、条件；Ｎ＝１０００、Δｔ＝０．１秒、Ａ＝１０、α＝−０．０２５、Ｘi0＝１．０＊１０-10、初期値：Ｘi(t0)＝１．０＊１０-10は固定し、確率変数（Ｒ1、Ｒ2）に関して、合計１００００個の独立な確率変数の組を用いて、計１００００個の独立なシミュレーション結果を得た。このシミュレーション結果において、“ある瞬間（tk）”の値：Ｘiq(tk)の集合｛Ｘiq(tk)：ｑ＝１〜１００００｝を作成し、この集合に含まれる１００００個の値を、大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を、図４に示す。
【００８５】
すなわち、細胞数１００００個の生細胞群に関して、“ある瞬間（tk）”における、各細胞（ｑ）中のｍＲＮＡ量（Ｘiq(tk)）の分布をシミュレーション結果に相当している。図４に示すプロット結果は、概ね直線で近似できる分布を示し、遺伝子の発現量ｋ（ｍＲＮＡ量）の確率分布Ｐ（ｋ）は、Ｐ（ｋ）∝ｋ-rで示される「べき乗分布」を示すことに相当している。
【００８６】
さらに、上記式（１’）に示すモデル函数系において、条件；Ｎ＝１０００、Δｔ＝０．１秒、Ａ＝１０、Ｘi0＝１．０＊１０-10、初期値：Ｘi(t0)＝１．０＊１０-10は固定し、閾値パラメータ（−α）は、α＝０．０、α＝−０．００２５、α＝−０．０５の３条件を選択し、確率変数（Ｒ1、Ｒ2）は、互いに独立とした上で、３種のシミュレーションを行った。図５に、α＝０．０、図６に、α＝−０．００２５、図７に、α＝−０．０５を選択した際のシミュレーション結果をそれぞれ示す。
【００８７】
まず、図５〜図７に示すシミュレーション結果では、いずれも、離散的なスパイク状のピークが示され、その時間間隔には、規則性は見出されない。また、各ピークの高さも大きな揺らぎを示している。この二つの点では、共通している。但し、各ピークの高さの平均的な高さは、図５〜図７に示す結果の間で、明確な相違があり、αの値が小さくなる（閾値パラメータ（−α）が大きくなる）と伴に、減少する傾向を示している。従って、ピークの高さが、１．０＊１０-8の水準を超えるピーク数も、αの値が小さくなる（閾値パラメータ（−α）が大きくなる）と伴に、減少する傾向を示している。
【００８８】
このシミュレーションされたｍＲＮＡ量（Ｘiq）の時間変化を、時間幅（Δｔ×Ｎ）に関して、数値積分した結果を、図８に対比して示す。この図８に示す積分値は、αの値が小さくなる（閾値パラメータ（−α）が大きくなる）と伴に、減少する傾向をより明確に示している。すなわち、上記式（１’）に示すモデル函数系において、閾値パラメータ（−α）は、遺伝子からの転写過程の抑制に相当する、発現制御の大小を反映するパラメータとなっていることが検証される。
【００８９】
（第二の実施態様）
本発明の第二の形態にかかるモデル化の方法に基づく、一つの実施態様を説明する。
【００９０】
個々の細胞内における前記ｍＲＮＡ量（Ｘi）を該ｍＲＮＡ分子数、すなわち、ｍＲＮＡ量（Ｘi）は、自然数で示される場合、
個々の細胞内におけるｍＲＮＡ量（Ｘi）の変化量：ΔＸiを前記式（II）により表記する際、時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝α＋Ｒ1
ｋdegrade＝β＋Ｒ2
また、時間的に変動する定数Ｃとして、
Ｃ＝Ｒ3
（但し、
αは、α≧０  の条件を満たす定数、
βは、β≧０  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝λ1であり、ポワソン分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝λ2であり、ポワソン分布に従う確率変数、
Ｒ3は、平均値＜Ｒ3＞は、＜Ｒ3＞＝λ3であり、ポワソン分布に従う確率変数であって、
（β＋λ2）≧（α＋λ1）  ならびに  ０≧（α＋λ1）−（β＋λ2）＋λ3  の条件を満たす）
を選択する。その際、上記式（II）は、下記式（２）：
式（２）：    ΔＸi＝（α＋Ｒ1）＊Ｘi−（β＋Ｒ2）＊Ｘi＋Ｒ3
で示される。この式（２）は、
ΔＸi＝｛（α−β）＋（Ｒ1−Ｒ2）｝＊Ｘi＋Ｒ3
と変形できる。すなわち、ランダム変数Ｒ1、Ｒ2、Ｒ3は、ポワソン分布の確率密度関数：ｆ（Ｒ）＝ｅ-R＊Ｒλ／（λ！）（λは、自然数を表す）に従う確率変数である。また、α、βは、整数値をとる定数である。
【００９１】
その際、ランダム変数Ｒ1、Ｒ2の差違（Ｒ1−Ｒ2）が、閾値パラメータ（−（α−β））の値を超える時間のみで、ΔＸiのｍＲＮＡ量（Ｘi）に依存する増加が起こり、閾値パラメータ（−（α−β））の値を下回ると、ｍＲＮＡ量（Ｘi）依存的な減少に転じる。従って、上記式（１’）に示すモデル函数系と同様に、この式（２）で示すモデル函数系でも、ランダム変数Ｒ1、Ｒ2の差違（Ｒ1−Ｒ2）に依存して、離散的なシャープなピークを示す時間変動を示すものとなる。一方、大半の時間では、｛（α−β）＋（Ｒ1−Ｒ2）｝は、負の値を示すため、ΔＸiは、「０」あるいは、負の値となり、ｍＲＮＡ量（Ｘi）は、低水準で僅かに変動している状態となる。
【００９２】
式（２）に示すモデル函数系を用いる場合、前記離散的なシャープなピークを示す際、そのピーク付近の時刻においては、｛（α−β）＋（Ｒ1−Ｒ2）｝が正負の間で振動すると、｛（α−β）＋（Ｒ1−Ｒ2）｝＊Ｘiの値が正負の間を大きく変動する。一方、前記離散的なシャープなピークを除く、多くの時間帯では、｛（α−β）＋（Ｒ1−Ｒ2）｝＊Ｘiの値は、負の値となっているが、ランダム変数Ｒ3を加えることで、ΔＸiは、「０」の前後で僅かに変動する状態とでき、下限領域に対する「緩衝成分」項に相当している。
【００９３】
ランダム変数（Ｒ1、Ｒ2、Ｒ3）のセットは、各細胞において、Ｒ1、Ｒ2、Ｒ3相互を独立に、また、細胞間においても、（Ｒ1、Ｒ2、Ｒ3）のセットは、相互に独立に選択する。
【００９４】
図９に、上記式（２）に示すモデル函数系を用いたシミュレーション結果の一例を示す。図９に示す例では、シミュレーションする時間幅（Δｔ）×Ｎは、Ｎ＝１０００、Δｔ＝０．１秒、定数αは、α＝０、βは、β＝０をそれぞれ設定している。また、ランダム変数（Ｒ1、Ｒ2、Ｒ3）に関して、利用されるポワソン分布の確率密度関数：ｆ（Ｒ）として、Ｒ1には、その平均λ1がλ1＝１、Ｒ2には、その平均λ2がλ2＝２、Ｒ3には、その平均λ3がλ3＝１を選択している。また、Ｘiの最低値：Ｘi0は、Ｘi0＝１を設定し、初期値：Ｘi(t0)として、Ｘi(t0)＝１としている。なお、数値演算の過程で、Ｘi(t)が、Ｘi(t)≦０となった場合、最低値Ｘi0：１の値に置き換えるという境界条件を選択している。
【００９５】
なお、上記式（２）に示すモデル函数系では、閾値パラメータ（−（α−β））の値によっては、Ｘiの値が、一定水準を超えて高くなると、その後、不必要に大きな振幅で振動する状態が出現することがある。この巨大振幅の振動現象を回避するため、Ｘiの値に上限を設けている。図９に示す例では、このＸiの上限値：ＸiLIMITとして、ＸiLIMIT＝１００，０００を設定している。
【００９６】
図９に示す結果には、離散的なシャープなピークを示す、時間的な変動が得られている。一方、離散的なシャープなピークのピーク高さは、１００を超えてなく、予定した最低値：Ｘi0＝１と比較して、ピーク高さは有意に高いとは言えないものであった。このような低い水準で推移する、ｍＲＮＡ量（Ｘi）の時間変動は、遺伝子の転写が抑制されている状況に相当している。
【００９７】
図９に示す結果に関して、各時刻（tk）における、Ｘi(tk)の値を大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を、図１０に示す。各時刻（tk）の半分以上では、Ｘi(tk)の値は、１となっており、概ね半分の時間では、遺伝子の転写はなく、Ｘi≧１０とっている延べ時間は、３０／１０００程度の比率となっている。図１０に示すプロット結果も、概ね直線で近似できる分布を示し、遺伝子の発現量ｋ（ｍＲＮＡ量）の確率分布Ｐ（ｋ）は、Ｐ（ｋ）∝ｋ-rで示される「べき乗分布」を示すことに相当している。
【００９８】
値ｋをとる確率分布ｐ（ｋ）が、Ｐ（ｋ）∝ｋ-rで示される「べき乗分布」のうち、べき乗数；ｒ＝２であるものは、「Ｚｉｐｆの法則」と称される現象に相当するものである。群に含まれる要素数がＮTOTALの時、その中に、値ｋを示すものの個数は、ｎ（ｋ）＝ＮTOTAL＊Ｐ（ｋ）となる。従って、値ｋi以上を示す要素数の総和；Ｎ（ｋi）は、
Ｎ（ｋ）＝∫kＮTOTAL＊Ｐ（ｋ）ｄｋ
＝［ＮTOTAL＊Π（ｋ）］k＝ＮTOTAL＊｛Π（∞）−Π（ｋ）｝
となる。その際、Ｐ（ｋ）∝ｋ-2であれば、Π（ｋ）∝ｋ-1となり、Ｎ（ｋ）∝ｋ-1の依存性を示すものとなる。従って、両対数プロット上、ｌｏｇｋ＝Ｃ＋（−１）ｌｏｇＮ（ｋ）で示される関係となる。
【００９９】
一方、図１０に示すプロット結果は、ｌｏｇＸi＝Ｃ＋（−１／２）ｌｏｇＮ（Ｘi）で示される関係と近似でき、順位Ｎ（Ｘi）とＸiとの関係は、Ｎ（Ｘi）∝Ｘi-2、換言すると、確率分布Ｐ（Ｘi）としては、Ｐ（Ｘi）∝Ｘi-3程度の「べき乗分布」を示すことに相当している。
【０１００】
図１１に示す例では、シミュレーションする時間幅（Δｔ）×Ｎは、Ｎ＝１０００、Δｔ＝０．１秒、定数αは、α＝１、βは、β＝０をそれぞれ設定している。また、ランダム変数（Ｒ1、Ｒ2、Ｒ3）に関して、利用されるポワソン分布の確率密度関数：ｆ（Ｒ）として、Ｒ1には、その平均λ1がλ1＝１、Ｒ2には、その平均λ2がλ2＝２、Ｒ3には、その平均λ3がλ3＝１を選択している。また、Ｘiの最低値：Ｘi0は、Ｘi0＝１を設定し、初期値：Ｘi(t0)として、Ｘi(t0)＝１としている。なお、数値演算の過程で、Ｘi(t)が、Ｘi(t)≦０となった場合、最低値Ｘi0：１の値に置き換えるという境界条件を選択している。図１１に示す例では、このＸiの上限値：ＸiLIMITとして、ＸiLIMIT＝１００，０００を設定している。
【０１０１】
図１１に示す結果でも、離散的なシャープなピークを示す、時間的な変動が得られている。一方、離散的なシャープなピークのピーク高さは、１０，０００を超え、予定した最低値：Ｘi0＝１と比較して、ピーク高さが有意に高いものも存在している。従って、図９に示す遺伝子の転写が抑制されている状況と比較すると、ある時期には、遺伝子の転写が顕著に生じている状況、全体として、遺伝子の発現はある程度活発化していることに相当する。
【０１０２】
図１１に示す結果に関して、各時刻（tk）における、Ｘi(tk)の値を大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を、図１２に示す。Ｘi(tk)の値が、１００以上である延べ時間は、５０／１０００程度の比率でしかなく、殆どの時間では、低水準となっている。図１２に示すプロット結果も、概ね直線で近似できる分布を示し、遺伝子の発現量ｋ（ｍＲＮＡ量）の確率分布Ｐ（ｋ）は、Ｐ（ｋ）∝ｋ-rで示される「べき乗分布」を示すことに相当している。
【０１０３】
図１２に示すプロット結果は、ｌｏｇＸi＝Ｃ＋（−２）ｌｏｇＮ（Ｘi）で示される関係と近似でき、順位Ｎ（Ｘi）とＸiとの関係は、Ｎ（Ｘi）∝Ｘi-1/2、換言すると、確率分布Ｐ（Ｘi）としては、Ｐ（Ｘi）∝Ｘi-3/2程度の「べき乗分布」を示すことに相当している。
【０１０４】
図１３に示す例では、シミュレーションする時間幅（Δｔ）×Ｎは、Ｎ＝１０００、Δｔ＝０．１秒、定数αは、α＝２、βは、β＝０をそれぞれ設定している。また、ランダム変数（Ｒ1、Ｒ2、Ｒ3）に関して、利用されるポワソン分布の確率密度関数：ｆ（Ｒ）として、Ｒ1には、その平均λ1がλ1＝１、Ｒ2には、その平均λ2がλ2＝２、Ｒ3には、その平均λ3がλ3＝１を選択している。また、Ｘiの最低値：Ｘi0は、Ｘi0＝１を設定し、初期値：Ｘi(t0)として、Ｘi(t0)＝１としている。なお、数値演算の過程で、Ｘi(t)が、Ｘi(t)≦０となった場合、最低値Ｘi0：１の値に置き換えるという境界条件を選択している。図１１に示す例では、このＸiの上限値：ＸiLIMITとして、ＸiLIMIT＝１００，０００を設定している。
【０１０５】
図１３に示す結果では、離散的なシャープなピークを示す、時間的な変動を示すが、離散的なシャープなピークのピーク高さは、１０００，０００を超え、予定した上限値：ＸiLIMITを超えている時間が相当に存在している。従って、図１３に示す結果は、頻繁に、遺伝子の転写が顕著に生じている状況、全体として、遺伝子の発現は高い水準まで活発化していることに相当する。
【０１０６】
図１３に示す結果に関して、各時刻（tk）における、Ｘi(tk)の値を大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を、図１４に示す。Ｘi(tk)の値が、１００以上である延べ時間は、３００／１０００程度の比率まで増しており、全体的に高い水準となっている。一方、Ｘi(tk)の値が、最低値Ｘi0：１の値を示す比率は、２００／１０００程度の比率で残っている。図１４に示すプロット結果でも、概ね直線で近似できる分布が示されており、遺伝子の発現量ｋ（ｍＲＮＡ量）の確率分布Ｐ（ｋ）は、Ｐ（ｋ）∝ｋ-rで示される「べき乗分布」を示すことに相当している。
【０１０７】
図９、図１１、図１３に示すシミュレーション結果では、いずれも、離散的なスパイク状のピークが示され、その時間間隔には、規則性は見出されない。また、各ピークの高さも大きな揺らぎを示している。この二つの点では、共通している。但し、各ピークの高さの平均的な高さは、図９、図１１、図１３に示す結果の間で、明確な相違があり、（α−β）の値が大きくなる（閾値パラメータ（−（α−β））が小さくなる）と伴に、増大する傾向を示している。すなわち、上記式（２）に示すモデル函数系において、閾値パラメータ（α−β）は、遺伝子からの転写過程の抑制・促進に相当する、発現制御の大小を反映するパラメータとなっていることが検証される。さらには、パラメータαは、遺伝子からの転写過程の促進、パラメータβは、遺伝子からの転写過程の抑制にそれぞれ相当する、二種の発現制御機構による寄与の大小を反映するパラメータとなっていることが検証される。
【０１０８】
図１０、図１２、図１４に示す結果は、遺伝子の発現量ｋ（ｍＲＮＡ量）の確率分布Ｐ（ｋ）は、Ｐ（ｋ）∝ｋ-rで示される「べき乗分布」を示す系に対するモデル化が、上記式（２）に示すモデル函数系を用いることでも達成されることの査証となっている。
【０１０９】
上で説明した通り、上記式（２）に示すモデル函数系は、「独立」した確率変数として、ランダム変数Ｒ1、Ｒ2、Ｒ3を有しているが、表現すべき個々の細胞内におけるｍＲＮＡ量（Ｘi）の変化は、特に、スパイク状のピークを示す頻度、ならびに、各スパイク状のピークのピーク値に関しては、ランダム変数Ｒ1、Ｒ2によって、その大きな揺らぎがもたらされている。すなわち、表現すべき個々の細胞内におけるｍＲＮＡ量（Ｘi）の変化は、スパイク状のピークを示す頻度、ならびに、各スパイク状のピークのピーク値の双方とも、確率分布に従って揺らぎを示すものであり、少なくとも、ランダム変数Ｒ1、Ｒ2という「独立」した確率変数２を含む式（２）の表記は、目標とするモデル化に適合している。
【０１１０】
（第三の実施態様）
本発明の第二の形態にかかるモデル化の方法に基づく、別の実施態様を説明する。
【０１１１】
個々の細胞内における前記ｍＲＮＡ量（Ｘi）を該ｍＲＮＡ分子数、すなわち、ｍＲＮＡ量（Ｘi）は、自然数で示される場合、
個々の細胞内におけるｍＲＮＡ量（Ｘi）の変化量：ΔＸiを前記式（II）により表記する際、時間的に変動する速度定数ｋgenerateとｋdegradeとして、
ｋgenerate＝α1＋α2＊Ｒ1
ｋdegrade＝β1＋β2＊Ｒ2
前記時間的に変動する定数Ｃとして、
Ｃ＝Ｒ3
（但し、
α1は、α1≧０  の条件を満たす定数、α2は、α2＞０  の条件を満たす定数、
β1は、β1≧０  の条件を満たす定数、β2は、β2＞０  の条件を満たす定数、
Ｒ1は、平均値＜Ｒ1＞は、＜Ｒ1＞＝λ1であり、ポワソン分布に従う確率変数、
Ｒ2は、平均値＜Ｒ2＞は、＜Ｒ2＞＝λ2であり、ポワソン分布に従う確率変数、
Ｒ3は、平均値＜Ｒ3＞は、＜Ｒ3＞＝λ3であり、ポワソン分布に従う確率変数であって、
（β1＋β2＊λ2）≧（α1＋α2＊λ1）  ならびに  ０≧（α1＋α2＊λ1）−（β1＋β2＊λ2）＋λ3  の条件を満たす）
を選択する。その際、上記式（II）は、下記式（３）：
式（３）：    ΔＸi＝（α1＋α2＊Ｒ1）＊Ｘi−（β1＋β2＊Ｒ2）＊Ｘi＋Ｒ3
で示される。この式（３）は、
ΔＸi＝｛（α1−β1）＋（α2＊Ｒ1−β2＊Ｒ2）｝＊Ｘi＋Ｒ3
と変形できる。
【０１１２】
この式（３）に示すモデル函数系は、上記式（２）に示すモデル函数系を、より一般化したものに相当している。
【０１１３】
上で説明した通り、式（２）に示すモデル函数系と同様に、この式（３）に示すモデル函数系も、「独立」した確率変数として、ランダム変数Ｒ1、Ｒ2、Ｒ3を有しているが、表現すべき個々の細胞内におけるｍＲＮＡ量（Ｘi）の変化は、特に、スパイク状のピークを示す頻度、ならびに、各スパイク状のピークのピーク値に関しては、ランダム変数Ｒ1、Ｒ2によって、その大きな揺らぎがもたらされている。すなわち、表現すべき個々の細胞内におけるｍＲＮＡ量（Ｘi）の変化は、スパイク状のピークを示す頻度、ならびに、各スパイク状ピークのピーク値の双方とも、大きな揺らぎを示すものであり、少なくとも、ランダム変数Ｒ1、Ｒ2という「独立」した確率変数２を含む式（３）の表記も、目標とするモデル化に適合している。
【０１１４】
また、式（３）に示すモデル函数系においても、閾値パラメータ（α−β）は、遺伝子からの転写過程の抑制・促進に相当する、発現制御の大小を反映する主パラメータとして機能する。加えて、（α2＊Ｒ1−β2＊Ｒ2）という確率変数部に含まれる、副パラメータα2は、遺伝子からの転写過程の促進、副パラメータβ2は、遺伝子からの転写過程の抑制に、それぞれ対応する副次的な寄与を持っている。従って、遺伝子からの転写過程の促進は、主パラメータα1と、副パラメータα2による、二種の制御メカニズムが、また、遺伝子からの転写過程の抑制は、主パラメータβ1と、副パラメータβ2による、二種の制御メカニズムが、それぞれ独立に機能するような、より適用範囲が広いものとなっている。
【０１１５】
例えば、上記式（２）に示すモデル函数系では、閾値パラメータ（−（α−β））の値によっては、Ｘiの値が、一定水準を超えて高くなると、その後、不必要に大きな振幅で振動する状態が出現することがある。この巨大振幅の振動現象を回避するため、Ｘiの値に上限を設けている。図１３に示す例では、このＸiの上限値：ＸiLIMITとして、ＸiLIMIT＝１００，０００を設定しているが、相当数の時刻では、この上限値：ＸiLIMITによって、制限を受けている。式（３）に示すモデル函数系では、（α2＊Ｒ1−β2＊Ｒ2）という確率変数部に含まれる、副パラメータα2と副パラメータβ2を適正に選択することで、不必要に大きな振幅で振動する状態が出現するＸiの値の上限をより高くすることが可能となる。従って、上記式（２）に示すモデル函数系における、図１３に示す例と同等のシミュレーションを行う際、式（３）に示すモデル函数系を用い、副パラメータα2と副パラメータβ2を適正に選択することで、Ｘiの上限値：ＸiLIMITとして、より大きな値を選択することが可能となる。
【０１１６】
（第四の実施態様）
本発明にかかるモデル化方法では、上述する第一の形態にかかるモデル化方法、あるいは、第二の形態にかかるモデル化方法のように、ｍＲＮＡ量（Ｘi）の変化量ΔＸiを表現する函数中、ｍＲＮＡ量（Ｘi）に依存する成分として、（ｋgenerate＊Ｘi−ｋdegrade＊Ｘi）を含む形式とする、ｍＲＮＡ量（Ｘi）の変化量ΔＸiの時間変動を「モデル化する」モデル函数系を利用する手法以外に、ｍＲＮＡ量（Ｘi）の時間変動を直接、確率変数を利用して、表記する手法も適用できる。
【０１１７】
この本発明にかかる第三の形態にかかるモデル化方法では、個々の細胞内における前記ｍＲＮＡ量の時間的な変動を確率論的に推定可能な形態で表記するモデルとして、
（ａ）  一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生する。
（ｂ）  一つの細胞内における該特定の遺伝子の転写事象発生一回に付随して、ｍＲＮＡ量の過渡的な増加とその後の減少を伴う、シャープなピーク状のｍＲＮＡ量の時間的な変化が生じている。
（ｃ’）  一つの細胞内において、前後する二回の転写事象発生間における時間間隔は、実質的にランダムとなっており、また、各シャープなピーク状のｍＲＮＡ量の時間的な変化（ピーク高さ）も大きく揺らいでいる。
この３条件と、実質的に同等な条件として、下記の条件を満足するモデルとしている。
【０１１８】
先ず、シミュレーションする時間を、極く短い時間間隔（Δｔ）×Ｎ区分に分割し、時間（ｔ）を、ｔk＝ｋ＊Δｔ（ｋ＝０、１、２、…、Ｎ）の数列とし、各細胞（ｑ）に対して、時刻（ｔk）における、その細胞内におけるｍＲＮＡ量（Ｘiq(tk)）の数列を定義し、各細胞（ｑ）におけるｍＲＮＡ量（Ｘi）の時間変動とする。
【０１１９】
その際、前記のシャープなピークとして認識可能な、一定値以上のピーク高さを有する変化は、非常に散発的に発生している。すなわち、Ｎ区分のうち、その極く一部のみで、ｍＲＮＡ量（Ｘiq(tk)）は、一定値以上の大きさを示すが、残る大半は、ｍＲＮＡ量（Ｘiq(tk)）は、極く小さな値である。
【０１２０】
図１５に、この一定値以上のピーク高さを有するピークが、非常に散発的に発生している状態を模式的に示す。その際、ピーク間の時間間隔Δｔp-pは、Δｔp-p＝ｘである、確率分布は、ｆ（ｘ）∝ｅ-2|x-τ|/τ（τは、ピーク総数の期待値ｎ1で、全期間（Δｔ×Ｎ）を除した値である）で記述される指数分布に従うと仮定する。一方、各ピークのピーク高さ（Ｘiq(tk)）は、ｆ（ｘ）∝１／（１＋ｘ）2で記述される第２種のパレート分布（べき乗分布）に従うと仮定する。
【０１２１】
この条件に従って、Ｎ区分のうち、その極く一部のみに、一定値以上の大きさ（有意なピーク値：１００以上）を示すｍＲＮＡ量（Ｘiq(tk)）を、前述の確率函数に従って設定した。図１６に、この手法に基づき、生成されたｍＲＮＡ量（Ｘiq(tk)）の数列（時間変化）を図示した結果を示す。図１７には、図１６に記載されるｍＲＮＡ量（Ｘiq(tk)）の数列に関して、大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を示す。
【０１２２】
図１７におけるプロットにおいて、Ｘiq(tk)≧１０の領域は、概ね直線で近似できる分布が示されており、遺伝子の発現量ｋ（ｍＲＮＡ量）の確率分布Ｐ（ｋ）は、Ｐ（ｋ）∝ｋ-rで示される「べき乗分布」と実質的に等価なものとなっている。
【０１２３】
なお、各ピークのピーク高さ（Ｘiq(tk)）として、上述のパレート分布（べき乗分布）と同様の裾広がりを有する分布を選択する限り、ピーク間の時間間隔Δｔp-pは、Δｔp-p＝ｘである、確率分布は、例えば、平均値τ（τは、ピーク総数の期待値ｎ1で、全期間（Δｔ×Ｎ）を除した値である）の正規分布などを用いても、同等の結果となる。その際、各ピークのピーク高さ（Ｘiq(tk)）の分布は、一定値以上の大きさ（例えば、有意なピーク値：１００以上）の確率分布が、上述のパレート分布（べき乗分布）と同様であれば、小さな値の領域（例えば、１０以下）の分布の形状は、全体の傾向に大きな影響を与えない。
【０１２４】
この第三の形態にかかるモデル化方法では、遺伝子発現の大小は、一定値以上のピーク高さを有するピークの総数の期待値ｎ１と、さらには、各ピークのピーク高さ（Ｘiq(tk)）に関する確率分布の係数を制御パラメータとして利用することで、モデル化がなされる。
【０１２５】
次いで、本発明にかかるｍＲＮＡの定量方法について、より詳しく説明する。
【０１２６】
同一種類の細胞複数で構成される生細胞群において、
（ａ）  一つの細胞内における該特定の遺伝子の転写事象は、他の細胞内における該特定の遺伝子の転写事象とは独立して発生する。
という条件が成り立つと、全体の系が均質であっても、各時刻における、個々の細胞における該特定の遺伝子の転写に起因するｍＲＮＡ量は、図４に示すように、細胞間において、大きな揺らぎを示す。一方、図５〜図７、図８の結果を参照すると、仮に、ある時間幅内において、離散的に起こる特定の遺伝子の転写に対応する、スパイク状のピークを示すｍＲＮＡ量（Ｘi）の変動があっても、各細胞内におけるｍＲＮＡ量（Ｘi）の時間変化を追跡し、一定の時間幅において、観測されるｍＲＮＡ量（Ｘi）の時間積分値を算定すると、各細胞において、特定の遺伝子の転写が活性に進行しているか、否かを定量的に評価することが可能であることが判る。
【０１２７】
しかし、個々の細胞において、おけるｍＲＮＡ量（Ｘi）の時間変化を追跡することは不可能であり、それに代わり、同一種類の細胞複数で構成される生細胞群から、その一部をサンプリングして、該サンプル群を構成する複数の細胞について、ある時刻におけるｍＲＮＡ量（Ｘi）の小計を用いて、生細胞群全体が、特定の遺伝子の転写が活性に進行しているか、否かを、高い確度で評価する必要がある。
【０１２８】
生細胞群全体について、ある時刻におけるｍＲＮＡ量（Ｘi）の合計を、この生細胞群から、その一部をサンプリングした、サンプル群におけるｍＲＮＡ量（Ｘi）の小計を用いて、高い確度で推定する方法について、以下に説明する。
【０１２９】
（第五の実施形態）
最近、ＤＮＡマイクロアレイやＳＡＧＥ  （ｓｅｒｉａｌ  ａｎａｌｙｓｉｓ  ｏｆ  ｇｅｎｅ  ｅｘｐｒｅｓｓｉｏｎ）などの方法が確立され、均質な状況下に置かれている細胞複数個からなる一つのサンプル群に付いて、その細胞内に存在する複数の遺伝子に関して、各遺伝子から発現されるｍＲＮＡ量の合計を、包括的に測定することが可能となった。しかし、同一の状況下にある多数個の生細胞群から、その一部を採取したサンプル群における測定結果を比較した際、測定結果に大きな差違が見出され、データの再現性に疑問が生じる場合がある。また、その測定結果に大きな差違が生じた場合、その差違を引き起こす要因の解釈は非常に難しいのが現実である。
【０１３０】
この第５の実施形態においては、ＤＮＡマイクロアレイやＳＡＧＥのような複数種のｍＲＮＡ量の同時測定に適する手法を適用し、複数個の細胞で構成されるサンプル群を採取し、このサンプル群の細胞中から抽出される、非常に多くの種類のｍＲＮＡに付いて、各ｍＲＮＡ量を測定する実験を行う場合を想定する。その際、測定される各ｍＲＮＡ量は、“ある時刻”において、該サンプル群を構成している複数個の細胞内に存在している各ｍＲＮＡ量を合計したものに相当する。この状況をシミュレーションするため、先ず、該サンプル群を構成しているＮcell個の細胞について、各細胞内に存在しているＮ種類の遺伝子から発現される各ｍＲＮＡ量の時間的変化を個々の細胞について、本発明にかかるモデル化を利用して、シミュレーションする。次いで、Ｎ種類の遺伝子から発現される各ｍＲＮＡ量に関して、“ある時刻”における各細胞内の存在量を前記Ｎcell個の細胞について合計し、該サンプル群全体のｍＲＮＡ量のシミュレーション値とする。
【０１３１】
前記の手法を利用して、独立にシミュレーションされた、二つのサンプル群全体のｍＲＮＡ量のシミュレーション値間で、どの程度の差違があるかを比較する。さらに、該サンプル群を構成している細胞数Ｎcellを変えた際、かかる二つのサンプル群全体のｍＲＮＡ量のシミュレーション値間で差違の程度が変化するかを比較する。
【０１３２】
図４に示す結果では、細胞数１０，０００個について、各時刻における、個々の細胞における特定の遺伝子の転写に起因するｍＲＮＡ量のシミュレーションに、本発明の第一の形態にかかる、式（１’）のモデル函数系を適用している。一方、本第５の実施態様では、本発明の第二の形態にかかる、式（２）のモデル函数系を適用して、個々の細胞における特定の遺伝子の転写に起因するｍＲＮＡ量のシミュレーションを行う事例を示す。
【０１３３】
その際、対象とする細胞は、多数種の遺伝子が存在しており、生細胞群全体の平均をとると、各遺伝子の発現の活発さには、差違が存在していると仮定する。この事例では、多数種の遺伝子のうち、１，０００種類の遺伝子を取り上げ、
式（２）：    ΔＸi＝（α＋Ｒ1）＊Ｘi−（β＋Ｒ2）＊Ｘi＋Ｒ3
のモデル函数系における、遺伝子からの転写過程の促進程度を反映するパラメータαに関して、前記１，０００種類の遺伝子について、２≧α＞０の範囲とし、図１８に示すような、第Ｎgene番目の遺伝子に対するパラメータαgene（Ｎgene）を、
αgene（Ｎgene）＝２＊Ｎgene／１，０００（Ｎgene＝１〜１，０００）
と選択してシミュレーションを行う。
【０１３４】
このパラメータα以外の条件、シミュレーションする時間幅（Δｔ）×Ｎ；Ｎ＝１０００、Δｔ＝０．１秒、パラメータβ；β＝０、ランダム変数（Ｒ1、Ｒ2、Ｒ3）に関して、利用されるポワソン分布の確率密度関数：ｆ（Ｒ）として、Ｒ1の平均λ1；λ1＝１、Ｒ2の平均λ2；λ2＝２、Ｒ3の平均λ3；λ3＝１を選択している。また、Ｘiの最低値：Ｘi0は、Ｘi0＝１を設定し、初期値：Ｘi(t0)として、Ｘi(t0)＝１としている。なお、数値演算の過程で、Ｘi(t)が、Ｘi(t)≦０となった場合、最低値Ｘi0：１の値に置き換えるという境界条件を選択している。一方、Ｘiの上限値：ＸiLIMITとして、ＸiLIMIT＝１００，０００を設定している。
【０１３５】
遺伝子からの転写過程の促進程度に相違を有する、前記１，０００種類の遺伝子に相当する、１，０００種のパラメータαgene（Ｎgene）：αgene（Ｎgene）＝２＊Ｎgene／１，０００（Ｎgene＝１〜１，０００）について、該サンプル群を構成している細胞の個数：Ｎcellに相当する、互いに独立なランダム変数（Ｒ1、Ｒ2、Ｒ3）の組；（Ｎcell）×（１，０００）組を２セットを生成する。各セットに関して、１，０００種のパラメータαgene（Ｎgene）について、それぞれＮcell個のｍＲＮＡ量Ｘiの時間的変化のシミュレーション結果を得る。二つのセットに関して、そのＮcell個のｍＲＮＡ量の時間的変化のシミュレーション結果における、ある時刻ｔkにおけるｍＲＮＡ量Ｘi（ｔk）の小計を算出し、各パラメータαgene（Ｎgene）に関して、第一のセットにおけるｍＲＮＡ量Ｘi（ｔk）を横軸として、第二のセットにおけるｍＲＮＡ量Ｘi（ｔk）を縦軸として、プロットする。
【０１３６】
図１９〜図２２は、上記の条件でシミュレーションされた、該サンプル群を構成している細胞の個数：Ｎcellが、Ｎcell＝１００、１，０００、１０，０００、１００，０００の４種のシミュレーション結果を示す。
【０１３７】
図１９に示すＮcell＝１００の結果では、サンプル群全体のｍＲＮＡ量のシミュレーション値は、第一のセットと第二のセット間において、大きな相違（揺らぎ）を示しているが、該サンプル群を構成している細胞の個数：Ｎcellが増すに従って、第一のセットと第二のセット間における相違は低減されている。特に、サンプル群全体のｍＲＮＡ量のシミュレーション値が大きくなる領域における、二つのセット間での相違（揺らぎ）は、該サンプル群を構成している細胞の個数：Ｎcellの増加とともに、急速に低減している。
【０１３８】
例えば、パラメータαgene＝２に相当する遺伝子は、上で示した図１４のシミュレーション結果のように、各細胞において、かかるシミュレーションする時間幅（Δｔ）×Ｎ；Ｎ＝１０００、Δｔ＝０．１秒の約１０％程度の時間において、高いｍＲＮＡ量（ＸiLIMIT＝１００，０００）を示している。換言すると、同一の状況下にある多数個の生細胞群（例えば、細胞総数Ｎcell＝１，０００，０００）中では、該パラメータαgene＝２に相当する遺伝子に関して、ある時刻ｔkにおいて、高いｍＲＮＡ量（ＸiLIMIT＝１００，０００）を示す細胞の占める比率は、約１０％程度となっている。一方、同一の状況下にある多数個の生細胞群（例えば、細胞総数Ｎcell＝１，０００，０００）中から、サンプル群として、細胞数Ｎcell＝１００程度の小さな部分集合を選別すると、かかる高いｍＲＮＡ量（ＸiLIMIT＝１００，０００）を示す細胞の個数は、期待値は、Ｎcellの１０％程度（Ｎcell／１０）であるは、その前後に、少なくとも（Ｎcell／１０）1/2程度の変動幅を示す可能性を示唆している。一方、サンプル群の細胞数Ｎcellの桁が増すとともに、期待値（Ｎcell／１０）に対して、その前後の変動幅は、（Ｎcell／１０）1/2程度の大きさである場合、期待値（Ｎcell／１０）：（Ｎcell／１０）1/2程度の変動幅の相対的な比率は、急速に減少していく。図１９〜図２２の結果は、例えば、ポワソン分布における期待値：Ｘi-expectedに対する、変動幅は、（Ｘi-expected）1/2に比例して変化する傾向と同等の傾向を反映している結果となっている。
【０１３９】
実際に、上記の条件において、パラメータβをβ＝０からβ＝１に変化させ、式（２）のモデル函数系における、閾値パラメータ（α−β）をより小さくした上で、Ｎcell＝１０，０００において、シミュレーションした結果を図２３に示す。パラメータβを増して、全体的に、各遺伝子からの発現を抑制した状況に変化させると、図２１に示すβ＝０における結果と比較し、図２３に示すβ＝１における結果では、二つのセット間での相違（揺らぎ）は、相対的に大きなものとなっている。この変化も、期待値：Ｘi-expectedに対する、変動幅は、（Ｘi-expected）1/2に比例して変化する傾向を反映しており、例えば、各遺伝子からの発現を抑制した状況では、サンプル群において測定されるｍＲＮＡ量の期待値：Ｘi-expectedは約１．５桁減少する際、（Ｘi-expected）1/2は、約０．８桁しか減少しないとすると、相対的に、二つのセット間での相違（揺らぎ）は大きくなる。
【０１４０】
図１９に示す、Ｎcell＝１００のサンプル群における第一のセットのシミュレーション結果、すなわち、前記１，０００種類の遺伝子に相当する、１，０００種のパラメータαgene（Ｎgene）に対する、サンプル群全体のｍＲＮＡ量のシミュレーション値を、そのシミュレーション値の大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を、図２４に示す。同様に、図２２に示す、Ｎcell＝１００，０００のサンプル群における第一のセットのシミュレーション結果、すなわち、前記１，０００種類の遺伝子に相当する、１，０００種のパラメータαgene（Ｎgene）に対する、サンプル群全体のｍＲＮＡ量のシミュレーション値を、そのシミュレーション値の大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を、図２５に示す。
【０１４１】
図２４において、第１００番目のｍＲＮＡ量のシミュレーション値は、４×１０5程度であり、その値は、Ｘiの上限値：ＸiLIMIT＝１×１０5を示す細胞数は２〜３個程度であることに相当する。一方、図２５において、第１００番目のｍＲＮＡ量のシミュレーション値は、６×１０8程度であり、その値は、Ｘiの上限値：ＸiLIMIT＝１×１０5を示す細胞数は、４〜５×１０3個程度であることに相当する。第１００番目のｍＲＮＡ量のシミュレーション値に寄与する、このＸiの上限値：ＸiLIMIT＝１×１０5を示す細胞数Ｎmajorに関する変動幅として、（Ｎmajor）1/2程度の幅を想定すると、Ｎcell＝１００のサンプル群に関しては、Ｎmajor＝２〜３個に対して、（Ｎmajor）1/2＝１．４〜１．８個程度の幅であり、Ｎcell＝１００，０００のサンプル群に関しては、Ｎmajor＝４〜５×１０3個に対して、（Ｎmajor）1/2＝６．５〜７．５×１０1個程度の幅であり、その相対的に変動幅が顕著に減少している傾向を反映する指標となっている。一方、図２５に示す結果においては、Ｎcell＝１００，０００のサンプル群に関しては、Ｎmajor＝２〜３となるｍＲＮＡ量のシミュレーション値が、６×１０8程度となる順位は、べき乗分布を示す領域となっている。
【０１４２】
図２４においては、第１００番目のｍＲＮＡ量のシミュレーション値でさえ、既にべき乗分布を示す領域の末端部に含まれるが、図２５においては、第２００番目のｍＲＮＡ量のシミュレーション値に達して、ようやくべき乗分布を示す領域の末端部に含まれる状態となっている。実際の測定においても、サンプル群を構成する細胞数Ｎcellに関して、そのサンプル群において観測される各ｍＲＮＡ量について、その観測量の多い順にソートした上で、横軸のソート順位に対して、縦軸にその観測値を、両対数プロットした結果において、べき乗分布を示す領域の末端と見なせる順位を選択する。その選択される順位数が一定の水準値を超えるように、サンプル群を構成する細胞数Ｎcellを増すことで、二つのサンプル群相互間における、対応する遺伝子の発現に起因するｍＲＮＡ量の測定値の相違（揺らぎ）を低減することが可能であると、以上の結果から推断される。
【０１４３】
（第六の実施形態）
本実施形態においては、それぞれ同一の状況下にある多数個の生細胞群二群から、それぞれサンプルを採取し、この二つのサンプルについて、ターゲットとする遺伝子の発現に由来するｍＲＮＡ量を測定し、二つサンプルにおけるｍＲＮＡ量の測定結果に基づき、母体群におけるターゲット遺伝子の発現の大小を検定する手法を示す。
【０１４４】
二つの母体群Ａ、Ｂからそれぞれ採取された、サンプルＡとサンプルＢとで測定されるｍＲＮＡ量を比較して、対象とするターゲット遺伝子の発現の大小関係；（母体群Ａの発現量）＞（母体群Ｂの発現量）を主張したいとする。その際、棄無仮説は、（母体群Ａの発現量）≦（母体群Ｂの発現量）となる。すなわち、サンプルＡとサンプルＢとで測定されるｍＲＮＡ量を比較した際、（サンプルＡのｍＲＮＡ量測定値）／（サンプルＢのｍＲＮＡ量測定値）の比が何倍を超えたときに、有意水準ｐで（母体群Ａの発現量）＞（母体群Ｂの発現量）が主張できるかを調べるためには、（母体群Ａの発現量）＝（母体群Ｂの発現量）とした時、つまり、サンプルＡとサンプルＢに含まれるｍＲＮＡ量が全く同じ確率分布に従っていると仮定した際に、（サンプルＡのｍＲＮＡ量測定値）／（サンプルＢのｍＲＮＡ量測定値）の比が何倍以上となる現象が、確率ｐで起こるか否かを調べればよい。
【０１４５】
二つの母体群Ａ，Ｂにおいて、ある指標値ｋは、全く同じ確率分布ｐ（ｋ）を有している際に、例えば、各母体群からそれぞれ無作為にサンプリングした、１００組サンプルにおける該指標値（ｋi（Ａ）、ｋi（Ｂ）；ｉ＝１〜１００）について、その比の集合｛ｋi（Ａ）／ｋi（Ｂ）；ｉ＝１〜１００｝｝を調べる。その際、集合｛ｋi（Ａ）／ｋi（Ｂ）；ｉ＝１〜１００｝｝のうち、（ｋi（Ａ）／ｋi（Ｂ））＞２を満足する要素数は、５個未満であったとすると、（ｋi（Ａ）／ｋi（Ｂ））＞２という事象が起こる確率は５％以下であると推断される。逆に、実際に、二つの母体群Ａ，Ｂについて、各母体群からそれぞれ無作為にサンプリングしたサンプルにおける該指標値（ｋ（Ａ）、ｋ（Ｂ））を比較した場合、（ｋ（Ａ）／ｋ（Ｂ））＞２であったとすると、有意水準５％で、「二つの母体群Ａ，Ｂは、指標値ｋに関して、全く同じ確率分布ｐ（ｋ）を有している」という棄無仮説を否定することができる。換言すると、有意水準５％で、「二つの母体群Ａ，Ｂは、指標値ｋに関して、異なる確率分布ｐ（ｋ）を有している」と推定することができる。
【０１４６】
以下、２種類の確率分布を想定して、前述と同様に、二つの母体群Ａ，Ｂにおいて、ある指標値ｋは、全く同じ確率分布ｐ（ｋ）を有している際に、例えば、各母体群からそれぞれ無作為にサンプリングした、１０００組サンプルにおける該指標値（ｋi（Ａ）、ｋi（Ｂ）；ｉ＝１〜１０００）について、その比の集合｛ｋi（Ａ）／ｋi（Ｂ）；ｉ＝１〜１０００｝｝をシミュレーションした結果を示す。
【０１４７】
第一の確率分布は、平均１０、標準偏差１の正規分布（Ｎ（１０，１））である。この正規分布に従う母体群は、（標準偏差／平均）が０．１のデータと同一視できる。第二の確率分布は、指標値ｋを示す確率密度関数は、ｆ（ｋ）＝１／（１＋ｋ）2で記述される第２種のパレート分布とする。
【０１４８】
各確率分布について、個々の要素における指標値ｋが該確率分布に従って変化している二つの母体群Ａ，Ｂを作成し、母体群Ａ，Ｂのそれぞれから無作為に指標値ｋをサンプリングして、１０００個のペア（ｋi（Ａ）、ｋi（Ｂ）；ｉ＝１〜１０００）を選択する。その後、その比の集合｛ｋi（Ａ）／ｋi（Ｂ）；ｉ＝１〜１０００｝｝を高い値の順にソートし、上位５０位、１０位、５位の比の値を決定する。
【０１４９】
上で述べたように、総数１０００個の要素からなる、比の集合｛ｋi（Ａ）／ｋi（Ｂ）；ｉ＝１〜１０００｝｝において、ｋi（Ａ）／ｋi（Ｂ）の値が上位５０位、１０位、５位の比の値を超える確率は、それぞれ５％、１％、０．５％以下であると推断される。
【０１５０】
上記の二種の確率分布について、シミュレーションした結果における、比ｋi（Ａ）／ｋi（Ｂ）の値に関して決定された、上位５０位、１０位、５位の値を図２６に表形式で示す。すなわち、図２６に示す比の値は、その値を超えている際には、５％、１％、０．５％の有意水準で、「二つの母体群Ａ，Ｂは、指標値ｋに関して、異なる確率分布ｐ（ｋ）を有している」と推定する上での、閾値と見なせる。
【０１５１】
指標値ｋが、ターゲットとする遺伝子の発現に由来するｍＲＮＡ量であり、対照とする母体群Ｂが、（標準偏差／平均）が０．１の正規分布と見なせるｍＲＮＡ量を示す場合、検定対象の母体群Ａと、対照とする母体群からそれぞれサンプリングした、二つのサンプルＡ，Ｂにおいて測定されるｍＲＮＡ量を比較し、（サンプルＡのｍＲＮＡ量測定値）／（サンプルＢのｍＲＮＡ量測定値）の比が、１．５を超えている際には、有意水準０．５％で、「対象とするターゲット遺伝子の発現の大小関係は、（母体群Ａの発現量）＞（母体群Ｂの発現量）である」と主張できる。一方、対照とする母体群Ｂが、上述するパレート分布のような裾の広い分布をしめしている場合、検定対象の母体群Ａと、対照とする母体群からそれぞれサンプリングした、二つのサンプルＡ，Ｂにおいて測定されるｍＲＮＡ量を比較し、（サンプルＡのｍＲＮＡ量測定値）／（サンプルＢのｍＲＮＡ量測定値）の比に基づき、有意水準５％で、「対象とするターゲット遺伝子の発現の大小関係は、（母体群Ａの発現量）＞（母体群Ｂの発現量）である」と主張するために、（サンプルＡのｍＲＮＡ量測定値）／（サンプルＢのｍＲＮＡ量測定値）＞９０という結果が必要となる。
【０１５２】
以上に例示するように、検定の対象となる「ターゲットとする遺伝子の発現に由来するｍＲＮＡ量」が裾の広い分布を示す場合には、従来、正規分布型の分布を対象とする際に利用されてきた、有意な差を示す「比（倍率）」と較べると、遥かに大きな「比（倍率）」を検定の基準（閾値）とすることが必要となる。
【０１５３】
（第七の実施形態）
本実施形態においては、同一の状況下にある多数個の生細胞群から、サンプルを採取し、そのサンプルについて、二つ遺伝子の発現に由来する、二つのｍＲＮＡ量を測定し、サンプルにおける二つのｍＲＮＡ量の測定結果に基づき、母体群における、この二つ遺伝子の発現に関して、その大小関係を検定する手法を示す。
【０１５４】
先に、第二の実施形態に示すように、同一の状況下にある多数個の細胞で構成される生細胞群において、個々の細胞中におけるターゲットとする遺伝子の発現に由来するｍＲＮＡ量の時間的変化が式（２）で示すモデル函数系で表記可能な場合、該遺伝子発現の活発さの相違は、式（２）のパラメータβ、λ1、λ2、λ3、Ｘiの上限値（ＸiLIMIT）を固定し、αを変化させることでシミュレーションが可能である。その際、図９、図１１、図１３に例示されるｍＲＮＡ量Ｘi（ｔk）の時間的変化に基づき、図１０、図１２、図１４に示されるＸi（ｔk）の分布プロットが得られ、そのべき乗分布に従う領域の傾きにも相違が見出される。つまり、ｍＲＮＡ量Ｘi（ｔk）の時間積分値のみでなく、ｍＲＮＡ量Ｘi（ｔk）の分布プロットの形状にも、該遺伝子発現の活発さの差違が反映すると予測される。
【０１５５】
二つ遺伝子の発現に由来する、二つのｍＲＮＡ量；Ｘi（ｔk）とＸj（ｔk）とが、ともに、式（２）で示すモデル函数系で表記可能であり、式（２）のパラメータβ、λ1、λ2、λ3と、Ｘi、Ｘjの上限値（ＸiLIMIT、ＸjLIMIT）は同じとし、両者の間でαのみを変化することでシミュレーションが可能である場合にも、両者の分布プロットの形状を比較すると、二つの遺伝子間における、発現の活発さの差違が反映されると予測される。
【０１５６】
同一の状況下にある多数個の細胞からなる生細胞群について、該細胞内において、二つの遺伝子Ａと遺伝子Ｂの発現がなされている際、この生細胞群全体として、この二つの遺伝子Ａと遺伝子Ｂの発現の活発さに有意な差違があるか否かは、下記の手法を利用することでも解析可能である。
【０１５７】
この生細胞群全体は、一定の期間、同一の状況下にあり、その期間内では、二つの遺伝子Ａと遺伝子Ｂの発現は、それぞれ生細胞群全体を平均すると、時間的な変動を示さないと仮定する。この生細胞群を構成する多数個の細胞から、“ある時刻（ｔk）”において、ある程度多数の細胞サンプルを採取し、各細胞サンプル中の遺伝子Ａの発現に由来するｍＲＮＡ量と、遺伝子Ａの発現に由来するｍＲＮＡ量とを測定する。例えば、細胞サンプル数として、１０００個のサンプルを採取した際、各細胞サンプルにおける遺伝子ＡのｍＲＮＡ量測定値の集合｛ＸA（ｍ）；ｍ＝１〜１０００｝と、遺伝子ＢのｍＲＮＡ量測定値の集合｛ＸB（ｍ）；ｍ＝１〜１０００｝が得られる。この二つの集合について、それぞれ、ｍＲＮＡ量測定値の大きい順にソートし、順位を横軸に、そのｍＲＮＡ量測定値を縦軸にして、両対数プロットを行う。この遺伝子ＡのｍＲＮＡ量測定値の集合｛ＸA（ｍ）；ｍ＝１〜１０００｝と、遺伝子ＢのｍＲＮＡ量測定値の集合｛ＸB（ｍ）；ｍ＝１〜１０００｝について、それぞれ作成される分布プロットは、図１０、図１２、図１４に示す分布プロットと同様に、少なくとも、ソート順位の後半部分は、べき乗分布に類する分布の形状を示すものとなる。この分布プロットにおいて、べき乗分布に従う領域における「傾き」、すなわち、べき乗分布型の確率分布  ｐ（ｋ）∝ｋ-r  における指数（−ｒ）は、上記第二の実施形態で説明したように、各遺伝子の発現が活発か、否かを反映する指標となっている。具体的には、分布プロットにおいて、べき乗分布に従う領域における「傾き」が急なものほど、その遺伝子の発現がより活発であると推定される。
【０１５８】
仮に、遺伝子ＡのｍＲＮＡ量測定値の集合｛ＸA（ｍ）；ｍ＝１〜１０００｝の分布プロットが、図１２に示すような「傾き」が急なもので、一方、遺伝子ＢのｍＲＮＡ量測定値の集合｛ＸB（ｍ）；ｍ＝１〜１０００｝の分布プロットが、図１０に示すような「傾き」がより緩やかなものであれば、この生細胞群全体として、遺伝子Ａと遺伝子Ｂについえ、その発現の活発さの大小関係は、有意に（遺伝子Ａの発現）＞（遺伝子Ｂの発現）であると推定される。
【０１５９】
また、二つの生細胞群において、特定の遺伝子に関して、その発現の活発さの大小関係を比較する際にも、各生細胞群から、それぞれ“ある時刻（ｔk）”において、ある程度多数の細胞サンプルを採取し、各細胞サンプル中の遺伝子Ａの発現に由来するｍＲＮＡ量を測定した上で、上記の手法に準じで、分布プロットにおける分布形状の対比を行うことで、発現の活発さの大小関係を推定することが可能である。
【０１６０】
（第八の実施形態）
同一の状況下にある多数個の細胞で構成される生細胞群について、“ある時刻（ｔk）”における、個々の細胞中におけるターゲットとする遺伝子の発現に由来するｍＲＮＡ量を測定し、各細胞のｍＲＮＡ量測定値の集合を作成できたならば、そのｍＲＮＡ量測定値の大きい順にソートし、順位を横軸に、そのｍＲＮＡ量測定値を縦軸にして、両対数プロットを行うと、分布の形状を推定することは可能である。例えば、各細胞のｍＲＮＡ量測定値の集合が、正規分布を示すものである場合、ｍＲＮＡ量測定値は、平均値を中心にし、その前後に、標準偏差σの３倍程度の範囲に分布するため、裾広がりの小さな分布となる。一方、各細胞のｍＲＮＡ量測定値の集合が、べき乗分布に従うと、ｍＲＮＡ量測定値は、裾広がりの大きな分布となる。
【０１６１】
但し、べき乗分布に従うと裾広がりの大きな分布を対象とする際には、各細胞のｍＲＮＡ量測定値の集合の要素数（細胞サンプル数）が、１００個程度では、必ずしも、その分布形状を反映する分布プロットが得られるとは限らない。一方、集合の要素数（細胞サンプル数）が１０，０００個程度となると、その分布形状を反映する分布プロットが得られるものの、各細胞サンプルについて、同じ確度で含有するｍＲＮＡ量を測定することは、測定するｍＲＮＡ量が少なくなるとともに、技術的困難さが増す。
【０１６２】
本実施形態では、同一の状況下にある多数個の細胞で構成される生細胞群が、“ある時刻（ｔk）”における、各細胞のｍＲＮＡ量の確率分布として、べき乗分布でよく近似できる分布を有するか否かを推定する手法を示す。
【０１６３】
まず、同一の状況下にある多数個の細胞で構成される生細胞群から、ある程度多量の細胞数（２×Ｎcell）を有するサンプル群を採取する。図２７に示すように、このサンプル群を無作為に二分割して、細胞数（（１／２）0×Ｎcell＝Ｎcell）の二つの部分サンプル群Ｎ（（１／２）0）とする。次いで、一方の部分サンプル群Ｎ1/1を、無作為に二分割して、細胞数（（１／２）1×Ｎcell）の二つの部分サンプル群Ｎ（（１／２）1）とする。その後、この二分割操作を繰り返すことで、一連の部分サンプル群；Ｎ（（１／２）0）、Ｎ（（１／２）1）、…、Ｎ（（１／２）m）が調製される。この一連の部分サンプル群について、それを構成する複数の細胞から、全ｍＲＮＡを抽出した上で、各部分サンプル群全体における、特定の遺伝子の発現に由来するｍＲＮＡ量を測定する。その際、細胞からのｍＲＮＡ抽出、増幅、ラベル化などのｍＲＮＡ量の定量的な測定を行う際の一連の生化学的操作は、各部分サンプル群で同一条件とする。
【０１６４】
当初採取した細胞数（２×Ｎcell）を有するサンプル群を構成する、各細胞のｍＲＮＡ量測定値の集合が、正規分布を示す場合、一連の部分サンプル群；Ｎ（（１／２）0）、Ｎ（（１／２）1）、…、Ｎ（（１／２）m）において測定されるｍＲＮＡ量測定値は、一般に、各部分サンプル群の細胞数（（１／２）m×Ｎcell；ｍ＝０、１、…）が次第に減少しているのと対応して、各部分サンプル群全体のｍＲＮＡ量合計の測定値は、次第に減少していく。
【０１６５】
一方、当初採取した細胞数（２×Ｎcell）を有するサンプル群を構成する、各細胞のｍＲＮＡ量測定値の集合が、べき乗分布を示す場合、各部分サンプル群の細胞数（（１／２）m×Ｎcell；ｍ＝０、１、…）が順次減少している場合であっても、一連の部分サンプル群において測定されるｍＲＮＡ量測定値は、場合によっては、「逆転」現象が起こることもある。
【０１６６】
この「逆転」現象が起こる確率は、各細胞のｍＲＮＡ量測定値の集合が、べき乗分布を示す場合と、各細胞のｍＲＮＡ量測定値の集合が、正規分布を示す場合とでは顕著に異なることを利用すると、仮に、「逆転」現象が起こっていれば、「各細胞のｍＲＮＡ量測定値の集合が、べき乗分布を示す」と推定することが可能となる。
【０１６７】
一方、「各細胞のｍＲＮＡ量測定値の集合が、べき乗分布を示す」場合であっても、一連の部分サンプル群；Ｎ（（１／２）0）、Ｎ（（１／２）1）、…、Ｎ（（１／２）m）において、その部分サンプル群を構成する細胞数が十分に大きな範囲では、一連の部分サンプル群において測定されるｍＲＮＡ量測定値は、ほぼその細胞数の減少に対応して減少する。この細胞数の減少と、部分サンプル群において測定されるｍＲＮＡ量測定値の減少とが、ほぼ比例関係を維持する範囲であれば、部分サンプル群において測定されるｍＲＮＡ量測定値を、細胞数で除した、「細胞あたりのｍＲＮＡ量平均」は、同一の状況下にある多数個の細胞で構成される生細胞群全体における、「細胞あたりのｍＲＮＡ量平均」をほぼ反映した値を示す。
【０１６８】
換言するならば、上述する一連の部分サンプル群；Ｎ（（１／２）0）、Ｎ（（１／２）1）、…、Ｎ（（１／２）m）において測定されるｍＲＮＡ量測定値において、少なくとも、部分サンプル群；Ｎ（（１／２）0）、…Ｎ（（１／２）4）の範囲では、細胞数の減少と、部分サンプル群において測定されるｍＲＮＡ量測定値の減少とが、ほぼ比例関係を維持するように、当初採取するサンプル群の細胞数（２×Ｎcell）を選択することで、部分サンプル群において測定されるｍＲＮＡ量測定値を、細胞数で除した、「細胞あたりのｍＲＮＡ量平均」に基づき、生細胞群全体における、「細胞あたりのｍＲＮＡ量平均」を高い確度で推定することが可能である。
【０１６９】
同一の状況下にある多数個の細胞で構成される生細胞群において、この群を構成する各細胞のｍＲＮＡ量の確率分布が、べき乗分布を示す場合、上述する第五の実施形態において説明したように、サンプル群の細胞数のうち、特定の遺伝子を活発に発現し、そのｍＲＮＡ量が高い水準となっている細胞数Ｎhighの期待値；Ｎhigh-expectedに対する、その変動幅は、例えば、ポワソン分布における期待値：Ｘi-expectedに対する、変動幅は、（Ｘi-expected）1/2に比例して変化する傾向と同等の傾向を反映すると判断される。
【０１７０】
上述する「逆転」現象が起こる確率が顕著化する条件は、部分サンプル群の細胞数（（１／２）m×Ｎcell）中における、ｍＲＮＡ量が高い水準となっている細胞数Ｎhighの期待値；Ｎhigh-expectedに対して、（Ｎhigh-expected）1/2程度の変動幅が想定される際、２（Ｎhigh-expected）1/2≧Ｎhigh-expectedとなる状態に相当する。すなわち、４≧Ｎhigh-expectedとなる状態である。十分に多数個の細胞で構成される生細胞群において、特定の遺伝子を活発に発現し、そのｍＲＮＡ量が高い水準となっている細胞数の比率をｐ（high）とすると、部分サンプル群の細胞数（（１／２）m×Ｎcell）に対する、期待値；Ｎhigh-expectedは、Ｎhigh-expected＝（（１／２）m×Ｎcell）×ｐ（high）と表すことが可能である。
【０１７１】
従って、（（１／２）m×Ｎcell）×ｐ（high）≧４が満たされる範囲では、「逆転」現象が起こる確率は顕著化しないと推定される。この「逆転」現象の抑制に対する、境界条件は、Ｎcell≧２m+2×（１／ｐ（high））と表現される。先に示すように、部分サンプル群；Ｎ（（１／２）0）、…Ｎ（（１／２）4）の範囲では、細胞数の減少と、部分サンプル群において測定されるｍＲＮＡ量測定値の減少とが、ほぼ比例関係を維持するためには、少なくとも、ｍ＝５において、前記の境界条件が満足されていること、すんわち、Ｎcell≧２7×（１／ｐ（high））を満足するように、当初採取するサンプル群の細胞数（２×Ｎcell）を選択することが望ましいことを意味している。
【０１７２】
従って、同一の状況下にある多数個の細胞で構成される生細胞群において、この群を構成する各細胞のｍＲＮＡ量の確率分布が、べき乗分布のように、大きな揺らぎを示す確率分布である場合、特定の遺伝子が活発に発現され、そのｍＲＮＡ量が高い水準となっている細胞数の比率をｐ（high）とすると、当初採取するサンプル群の細胞数（２×Ｎcell）を、（２×Ｎcell）≧２8×（１／ｐ（high））を満足するように選択することが、このサンプル群において測定される「細胞あたりのｍＲＮＡ量平均」に基づき、生細胞群全体における「細胞あたりのｍＲＮＡ量平均」を高い確度で推定する上でより好ましいと判断される。
【０１７３】
例えば、対象とする遺伝子が、ｐ（high）≧１／１００の高い発現頻度を示す場合、Ｎcell≧１２６×（１００）≒１０，０００の細胞数を示す、部分サンプル群における「細胞あたりのｍＲＮＡ量平均」に基づき、生細胞群全体における、「細胞あたりのｍＲＮＡ量平均」を高い確度で推定することが可能である。対象とする遺伝子毎に、ｐ（high）≧１／１００〜１／１０，０００，０００と変化するとともに、部分サンプル群における「細胞あたりのｍＲＮＡ量平均」に基づき、生細胞群全体における、「細胞あたりのｍＲＮＡ量平均」を高い確度で推定する上で望ましい、部分サンプル群の細胞数の基準は、Ｎcell≧１０，０００〜１，０００，０００，０００と推定される。
【０１７４】
実際には、細胞内においては、同時に複数種の遺伝子が発現されており、また、発現されている各遺伝子間で、その遺伝子の発現が活発化した際に、生成されるｍＲＮＡ量は本質的な差違を示すことも多い。すなわち、細胞内における、該遺伝子の発現に伴って、生合成されるタンパク質量が、全タンパク質の１０％程度と多量となるもののある一方で、生合成されるタンパク質量が元来全タンパク質の０．０１％程度と少量となるもののある。このように、遺伝子が活発に発現されている際に、該遺伝子の発現に由来する「細胞内のｍＲＮＡ量が高い水準にある」状態に質的な差違がある場合もある。上で述べた「特定の遺伝子が活発に発現され、そのｍＲＮＡ量が高い水準となっている細胞数の比率」は、対象とする遺伝子の種類に応じて、例えば、式（１）のモデル函数系を利用するシミュレーション結果を利用して、適宜推定することも可能である。
【０１７５】
仮に、単一の生細胞に関して、一定の期間、その細胞内に存在している、特定の遺伝子に由来するｍＲＮＡの量（Ｘi）を追跡したとすると、該ｍＲＮＡ量Ｘiは、定常的に一定水準の値の周辺で小さな変動幅で推移する状態ではなく、その期間内の限られた、離散的な複数の時期では、シャープなピークを示すｍＲＮＡ量Ｘiの上昇が見出されるが、残る期間においては、ｍＲＮＡ量Ｘiは、実質的に「０」の水準となる、散発的に遺伝子発現が起こる状態が観測される。なお、この散発的な遺伝子発現に起因する、シャープなピークを示すｍＲＮＡ量Ｘiの上昇量を、長期間にわたり積分し、その積分値を積分期間で除した「ｍＲＮＡ量Ｘiの時間平均値」は、散発的な遺伝子発現の頻度と、各散発的な遺伝子発現に伴い、生成されるｍＲＮＡ量の平均的な上限との積を反映している。
【０１７６】
十分に多数個の細胞で構成される生細胞群において、特定の遺伝子を活発に発現し、そのｍＲＮＡ量が高い水準となっている細胞数の比率；ｐ（high）は、前述する、一定の期間内において、散発的に該遺伝子の発現が起こっている頻度と、その発現が維持されている時間幅の積が占める比率に相当している。
【０１７７】
同一の状況下にある多数個の細胞で構成される生細胞群から、ある程度多量の細胞数（２×Ｎcell）を有するサンプル群を採取し、上述する一連の部分サンプル群を調製した上で、各部分サンプル群から全ＲＮＡを抽出した上で、各ｍＲＮＡ量を定量する場合、この測定される各ｍＲＮＡ量は、“ある時刻（ｔk）”（短時間幅）に各細胞中に存在するｍＲＮＡ量を合計した値に相当する。
【０１７８】
例えば、図５に示されるように、細胞内において、特定の遺伝子からの発現が活発に（高い頻度）で起こっている場合、ある程度多量の細胞数（２×Ｎcell）を有するサンプル群を採取すると、“ある時刻（ｔk）”（短時間幅）に、「細胞内のｍＲＮＡ量が高い水準となっている状態」の細胞数の比率ｐ（high）は、相当に高いものとなっている。従って、サンプル群の細胞数（２×Ｎcell）が比較的少ない場合でも、部分サンプル群における「細胞あたりのｍＲＮＡ量平均」に基づき、生細胞群全体における、「細胞あたりのｍＲＮＡ量平均」を高い確度で推定することが可能である。但し、細胞内においては、同時に複数種の遺伝子が発現されており、また、発現されている各遺伝子間で、その遺伝子の発現が活発化した際に、生成されるｍＲＮＡ量は本質的な差違を示すことも多いため、生細胞群全体における、「細胞あたりのｍＲＮＡ量平均」の絶対値は、対象とする遺伝子に依存して、基本的にバラツキを示す。しかし、“ある時刻（ｔk）”（短時間幅）における、散発的な遺伝子発現の頻度が同程度に高い場合には、二つの遺伝子間で、それぞれの遺伝子の発現に起因する、細胞内のｍＲＮＡ量について、「細胞あたりのｍＲＮＡ量平均」の比較を、サンプル群の細胞数（２×Ｎcell）が比較的に少ない場合でも、高い確度で実施することが可能である。
【０１７９】
一方、図６や図７に示されるように、細胞内において、特定の遺伝子からの発現が活発でない（低い頻度でしか起こっていない）場合には、サンプル群の細胞数（２×Ｎcell）が比較的に少ない際には、このサンプル群において測定されるｍＲＮＡ量自体、絶対的にも相対的にも大きなバラツキを示す恐れがある。つまり、ハウスキーピング遺伝子のような「どのような細胞でも、そこそこ発現していると考えられる遺伝子」であっても、“ある時刻（ｔk）”（短時間幅）に、「細胞内のｍＲＮＡ量が高い水準となっている状態」の細胞数の比率ｐ（high）が低い場合には、サンプル群の細胞数（２×Ｎcell）を非常に多くしないと、該サンプル群において測定される「細胞あたりのｍＲＮＡ量平均」は、個々のサンプル群の間で非常に大きなバラツキを示すことが予想される。従来、細胞数が少なくサンプルを用いた際、「全ての細胞で発現しているはずのハウスキーピング遺伝子のｍＲＮＡが、必ずしも全てのサンプルから検出されない」という実験結果は、この現象を反映していると考えられる。
【０１８０】
以下、“ある時刻（ｔk）”（短時間幅）における、散発的な遺伝子発現の頻度を考慮した上で、具体的に、どの程度の細胞数を有するサンプルを採取して、「細胞あたりのｍＲＮＡ量平均」の測定を行う必要があるかということを試算してみる。
【０１８１】
測定対象の遺伝子からの発現が、外部環境の変化に応じて、１００秒単位で応答するような応答速度の高い遺伝子である場合には、“ある時刻（ｔk）”（短時間幅）における、散発的な遺伝子発現の頻度は、極めて高くなり、サンプル群を構成する細胞数が１０００個程度の、非常に小さなサンプルを用いても、「細胞あたりのｍＲＮＡ量平均」の測定値は相対的に安定したものとなる。
【０１８２】
このように“ある時刻（ｔk）”（短時間幅）における、散発的な遺伝子発現の頻度は、極めて高くなる、応答性の高い遺伝子の例としては、例えば、大腸菌の代謝酵素の例がある。先ず、大腸菌自体、対数増殖期においては、セル・サイクル時間が短い。加えて、通常、大腸菌では、グルコース代謝回路が活性化しているが、ラクトースしかエネルギー源がない環境におかれると、速やかに、その代謝回路が切り替わる。その際、今まで発現していなかった、ラクトースを利用する代謝回路に関与する酵素タンパク質の遺伝子が、短期間内に高い水準で発現される。
【０１８３】
また、所謂、シグナル伝達系を介する、非常に速い細胞応答に対応して、細胞内で生合成されるタンパク質の遺伝子も、短期間内に高い水準で発現される。例えば、各種のホルモンの分泌量によって、組織内の細胞の状態に変化が生じる系では、ホルモン・レセプターを介するシグナル伝達系、あるいは、ホルモン結合型の転写因子タンパク質の活性化に付随し、これら転写因子タンパク質による転写促進を受ける、タンパク質の遺伝子の発現が急速に促進される。また、形態形成の過程において、積極的な細胞死（アポトーシス）が誘起されることもあるが、その際、かかるアポトーシスに関わる調整タンパク質の遺伝子も急速に発現する。これらの特定の時期に、外的な要因に起因して、発現の促進がなされる遺伝子は、“ある時刻（ｔk）”（短時間幅）における、散発的な遺伝子発現の頻度は、極めて高くなる、応答性の高い遺伝子の代表例に含まれる。
【０１８４】
それに対して、哺乳動物における、癌細胞のマーカー遺伝子の発現は、動物細胞自体、その増殖は、微生物細胞と較べて、格段に遅い。非癌細胞が癌化する過程では、数日をかけてマーカー遺伝子の発現量が増加していくことが考えられる。例えば、１１日ぐらいをかけて非癌細胞から癌細胞へ変化する過程では、“ある時刻（ｔk）”（短時間幅）、例えば、１００秒単位における、マーカー遺伝子の発現頻度は、極めて、小さな頻度となる。換言すると、癌化が進行している組織から、採取されるサンプルに含まれる細胞数は、１．０＊１０7個程度、つまり、１０ｍｇ以上の癌細胞を集めないと、マーカー遺伝子の発現を反映する「細胞あたりのｍＲＮＡ量平均」を相対的に安定して測定できないということになる。なお、非癌細胞が癌化する過程が、さらに緩やかに進行する場合には、マーカー遺伝子の発現を反映する「細胞あたりのｍＲＮＡ量平均」を相対的に安定して測定する上で必要となる、サンプルに含まれる細胞数は、さらに多くなる。
【０１８５】
同様に、ハウスキーピング遺伝子などの全ての細胞で発現しているが、全体の生細胞群のうち、“ある時刻（ｔk）”（短時間幅）において、実際に発現を行っている細胞数の比率が低い場合には、採取したサンプル群における「細胞あたりのｍＲＮＡ量平均」を相対的に安定的に測定するためには、非常に大量のサンプルが必要となる。例えば、細胞分裂の一環として、発現がなされるようなハウスキーピング遺伝子は、細胞増殖がほぼ完了し、安定期に達した生細胞群では、細胞分裂を起こしている細胞の比率は低減しており、この細胞分裂と同期して起こるハウスキーピング遺伝子の発現の頻度も低減する。その際、よく行われる「ハウスキーピング遺伝子のｍＲＮＡ量を基準として、外的要因によって、発現の促進がなされる、各種の遺伝子のｍＲＮＡ量の正規化を行う」という処理は、採取されるサンプル群の細胞数によっては、非常に大きな「系統的な誤差」を含む可能性がある。
【０１８６】
（その他の実施形態）
なお、本発明にかかるモデル化の方法に基づく、シミュレーションは、複数の機器（例えば、ホストコンピュータ、インタフェイス機器、リーダ、プリンタなど）から構成されるシステムにおいて、かかる確率論的なシミュレーションを、コンピュータ上において実現するための演算処理用制御プログラムを適用して、実行することができる。さらには、中心的なシステムに対して、情報ネットワーク網を介して、情報データの遣り取りが可能なネットワーク接続されている、個別的な装置（例えば、パーソナル・コンピュータなど）複数の余剰タスク時間を利用して、並列処理形式によって、全体の処理を実行することもできる。なお、シミュレーション規模によっては、単位時間当たりの処理能力が大きくない、個別的な装置（例えば、パーソナル・コンピュータなど）のみを用いて、全体の処理を実行することもできる。
【０１８７】
その際、本発明にかかるモデル化の方法に基づく、確率論的なシミュレーションは、例えば、上述する実施形態に説明する、各種のモデル函数系を適用して構築されるソフトウェアのプログラムコードを記録した記憶媒体を、個々のシステムあるいは装置に供給し、そのシステムあるいは装置中のコンピュータ（またはＣＰＵやＭＰＵ）が、該記憶媒体に格納されたプログラムコードを読出し、実行する形態とすることができる。あるいは、ネットワーク接続を介して、プログラム・ソースやデータを収納するサーバーにアクセスし、大容量の記憶媒体上の一部に収納されるプログラムコードを読出し、実行する形態とすることができる。
【０１８８】
通常、各種のモデル函数系を適用して構築されるソフトウェアのプログラムコードを記録した記憶媒体として供給する際、その記憶媒体として、例えば、フレキシブル・ディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、磁気テープ、不揮発性のメモリカード、ＲＯＭなどを利用することができる。
【０１８９】
また、コンピュータが読出したプログラムコードは、その実行に際して、利用される全ての情報処理・演算操作機能のプログラムを含み、他のソースを利用することなく、例えば、前述した各実施形態に記載する各種のシミュレーション機能の実現がなされる、「スタンド・アローン」形態とすることができる。あるいは、そのプログラムコードの指示に基づき、そのコンピュータ上で稼働しているＯＳ（オペレーティングシステム）など、他のソースを利用して、実際の処理の一部または全部が実行される、他のソースをリンクすることで、例えば、前述した実施形態に記載する各種のシミュレーション機能の実現がなされる形態であってもよい。
【０１９０】
すなわち、記憶媒体から読出されたプログラムコードが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書込まれた後、そのプログラムコードの指示に基づき、その機能拡張ボードや機能拡張ユニットに備わるＣＰＵなどで、付属されているＲＯＭ中に提供されている、各種のソースを利用して、実際の処理の一部または全部を行い、その処理操作の制御・情報の遣り取りの管理のみに、コンピュータ本体が関与する形態をとって、前述した実施形態に記載する各種のシミュレーション機能の実現がなされる場合も含まれることは言うまでもない。
【産業上の利用可能性】
【０１９１】
本発明は、全体が均質な条件下に置かれている、同一種類の生細胞群を対象とする際、各細胞内における、特定の遺伝子の転写に起因して生成されるｍＲＮＡ量が時間的に大きな揺らぎを示し、同時に、各細胞内で転写に起因して生成されるｍＲＮＡ量は、細胞相互においても、大きな揺らぎを示している系において、かかる生細胞群全体から採取されるサンプルを利用し、相対的に安定的に遺伝子発現量の定量を行う上で必要となるサンプル量を推定する際に利用できる。
【図面の簡単な説明】
【０１９２】
【図１】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化に対応する「ｍＲＮの生成・分解」過程のモデル化の概念を模式化して表す図である。
【図２】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数を利用する解析に利用される情報処理装置の基本構成を模式的に示すブロック図である。
【図３】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数として、第一の実施形態に例示する式（１）の函数型を利用して、シミュレーションされた、一細胞内における特定のｍＲＮＡ量の時間的な変化の一例を示す図である。
【図４】第一の実施形態に例示する式（１）の函数型を利用して、確率変数（Ｒ1、Ｒ2）として、独立な１００００種を選択して、シミュレーションされた、一細胞内における特定のｍＲＮＡ量（Ｘ）の時間的な変化、計１００００通りの計算結果において、ある特定の時刻におけるｍＲＮＡ量（Ｘ）の値について、値の大きな順にソートし、横軸をソート順位、縦軸のそのｍＲＮＡ量（Ｘ）の値として、両対数プロットした結果を示す図である。
【図５】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数として、第一の実施形態に例示する式（１）の函数型を利用し、パラメータαに関して、α＝０と設定した上でシミュレーションされた、一細胞内における特定のｍＲＮＡ量の時間的な変化の一例を示す図である。
【図６】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数として、第一の実施形態に例示する式（１）の函数型を利用し、パラメータαに関して、α＝−０．０２５と設定した上でシミュレーションされた、一細胞内における特定のｍＲＮＡ量の時間的な変化の一例を示す図である。
【図７】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数として、第一の実施形態に例示する式（１）の函数型を利用し、パラメータαに関して、α＝−０．０５０と設定した上でシミュレーションされた、一細胞内における特定のｍＲＮＡ量の時間的な変化の一例を示す図である。
【図８】図５〜図６に示す、一細胞内における特定のｍＲＮＡ量の時間的な変化のシミュレーション結果に関して、表示時間範囲のｍＲＮＡ量（Ｘ）の値を時間積分した値を対比する図である。
【図９】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数として、第二の実施形態に例示する式（２）の函数型を利用し、パラメータαに関して、α＝０と設定した上でシミュレーションされた、一細胞内における特定のｍＲＮＡ量の時間的な変化の一例を示す図である。
【図１０】図９に示す、一細胞内における特定のｍＲＮＡ量の時間的な変化に基づき、各時刻におけるｍＲＮＡ量（Ｘ）の値について、値の大きな順にソートし、横軸をソート順位、縦軸のそのｍＲＮＡ量（Ｘ）の値として、両対数プロットした結果を示す図である。
【図１１】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数として、第二の実施形態に例示する式（２）の函数型を利用し、パラメータαに関して、α＝１と設定した上でシミュレーションされた、一細胞内における特定のｍＲＮＡ量の時間的な変化の一例を示す図である。
【図１２】図１１に示す、一細胞内における特定のｍＲＮＡ量の時間的な変化に基づき、各時刻におけるｍＲＮＡ量（Ｘ）の値について、値の大きな順にソートし、横軸をソート順位、縦軸のそのｍＲＮＡ量（Ｘ）の値として、両対数プロットした結果を示す図である。
【図１３】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数として、第二の実施形態に例示する式（２）の函数型を利用し、パラメータαに関して、α＝２と設定した上でシミュレーションされた、一細胞内における特定のｍＲＮＡ量の時間的な変化の一例を示す図である。
【図１４】図１３に示す、一細胞内における特定のｍＲＮＡ量の時間的な変化に基づき、各時刻におけるｍＲＮＡ量（Ｘ）の値について、値の大きな順にソートし、横軸をソート順位、縦軸のそのｍＲＮＡ量（Ｘ）の値として、両対数プロットした結果を示す図である。
【図１５】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル化の第三の形態における、モデル化の概念を模式的に示す図である。
【図１６】本発明にかかる転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル化の第三の形態に基づき、シミュレーションされた、一細胞内における特定のｍＲＮＡ量の時間的な変化の一例を示す図である。
【図１７】図１６に示す、一細胞内における特定のｍＲＮＡ量の時間的な変化に基づき、各時刻におけるｍＲＮＡ量（Ｘ）の値について、値の大きな順にソートし、横軸をソート順位、縦軸のそのｍＲＮＡ量（Ｘ）の値として、両対数プロットした結果を示す図である。
【図１８】第五の実施形態において、シミュレーションに利用した、遺伝子からの転写過程の促進程度を反映するパラメータαに関して、１，０００種類の遺伝子について、２≧α＞０の範囲とし、図１８に示すような、第Ｎgene番目の遺伝子に対するパラメータαgene（Ｎgene）を、  αgene（Ｎgene）＝２＊Ｎgene／１，０００（Ｎgene＝１〜１，０００）と選択したパラメータ・セットを示す図である。
【図１９】第五の実施形態において、シミュレーションされた、各細胞数Ｎcell＝１００の第一のセットと第二のセット間において、パラメータαの一致するシミュレーション値相互を、第一のセットにおけるｍＲＮＡ量を横軸として、第二のセットにおけるｍＲＮＡ量を縦軸として、プロットした結果を示す図である。
【図２０】第五の実施形態において、シミュレーションされた、各細胞数Ｎcell＝１０００の第一のセットと第二のセット間において、パラメータαの一致するシミュレーション値相互を、第一のセットにおけるｍＲＮＡ量を横軸として、第二のセットにおけるｍＲＮＡ量を縦軸として、プロットした結果を示す図である。
【図２１】第五の実施形態において、シミュレーションされた、各細胞数Ｎcell＝１０，０００の第一のセットと第二のセット間において、パラメータαの一致するシミュレーション値相互を、第一のセットにおけるｍＲＮＡ量を横軸として、第二のセットにおけるｍＲＮＡ量を縦軸として、プロットした結果を示す図である。
【図２２】第五の実施形態において、シミュレーションされた、各細胞数Ｎcell＝１００，０００の第一のセットと第二のセット間において、パラメータαの一致するシミュレーション値相互を、第一のセットにおけるｍＲＮＡ量を横軸として、第二のセットにおけるｍＲＮＡ量を縦軸として、プロットした結果を示す図である。
【図２３】図２１の作製に際して用いた、パラメータβをβ＝０からβ＝１に変化させ、式（２）のモデル函数系における、閾値パラメータ（α−β）をより小さくした上で、第五の実施形態において、シミュレーションされた、各細胞数Ｎcell＝１０，０００の第一のセットと第二のセット間において、パラメータαの一致するシミュレーション値相互を、第一のセットにおけるｍＲＮＡ量を横軸として、第二のセットにおけるｍＲＮＡ量を縦軸として、プロットした結果を示す図である。
【図２４】図１９に示す、Ｎcell＝１００のサンプル群における第一のセットのシミュレーション結果、すなわち、前記１，０００種類の遺伝子に相当する、１，０００種のパラメータαgene（Ｎgene）に対する、サンプル群全体のｍＲＮＡ量のシミュレーション値を、そのシミュレーション値の大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を、
【図２５】図２２に示す、Ｎcell＝１００，０００のサンプル群における第一のセットのシミュレーション結果、すなわち、前記１，０００種類の遺伝子に相当する、１，０００種のパラメータαgene（Ｎgene）に対する、サンプル群全体のｍＲＮＡ量のシミュレーション値を、そのシミュレーション値の大きい順にソートした上で、横軸のソート順位に対して、縦軸にその値を、両対数プロットした結果を、
【図２６】第六の実施形態において、二種の確率分布を想定した上で、その値を超えている際には、５％、１％、０．５％の有意水準で、「二つの母体群Ａ，Ｂは、指標値ｋに関して、異なる確率分布ｐ（ｋ）を有している」と推定する上での、閾値をシミュレーションした結果を表形式で示す図である。
【図２７】第八の実施形態において説明する、確率分布の形態を特定する際に利用される手法を模式的に説明する図である。
【図２８】本発明にかかるモデル化方法に従って、例えば、転写産物（ｍＲＮＡ）量の時間的な変化を近似的に表記するモデル函数を利用する解析する際に利用される、情報処理操作の基本的構成の一例を模式的に示すブロック図である。
【符号の説明】
【０１９３】
１０１  転写過程
１０２  酵素的分解過程
１０３  翻訳過程
２０１  外部記憶
２０２  演算・制御部
２０３  内部メモリー
２０４  入出力装置
【図１】
2006101733.tif 000002
【図２】
2006101733.tif 000003
【図３】
2006101733.tif 000004
【図４】
2006101733.tif 000005
【図５】
2006101733.tif 000006
【図６】
2006101733.tif 000007
【図７】
2006101733.tif 000008
【図８】
2006101733.tif 000009
【図９】
2006101733.tif 000010
【図１０】
2006101733.tif 000011
【図１１】
2006101733.tif 000012
【図１２】
2006101733.tif 000013
【図１３】
2006101733.tif 000014
【図１４】
2006101733.tif 000015
【図１５】
2006101733.tif 000016
【図１６】
2006101733.tif 000017
【図１７】
2006101733.tif 000018
【図１８】
2006101733.tif 000019
【図１９】
2006101733.tif 000020
【図２０】
2006101733.tif 000021
【図２１】
2006101733.tif 000022
【図２２】
2006101733.tif 000023
【図２３】
2006101733.tif 000024
【図２４】
2006101733.tif 000025
【図２５】
2006101733.tif 000026
【図２６】
2006101733.tif 000027
【図２７】
2006101733.tif 000028
【図２８】
2006101733.tif 000029