(57)【要約】【課題】膜電位変化の時系列データといった間接情報を用いて、細胞が生きている状態下での細胞内に存在する分子群の生化学反応経路の推定方法及び装置を提供する。【解決手段】誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測し、膜電位変化を生成する生化学反応経路の全候補を数理モデルで表し、数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、事前確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出し、事前確率分布で指定するパラメータの範囲で膜電位変化と同様な変化を示す数理モデルを抽出し、生化学反応経路を推定する。【選択図】図１
2019154299.tif 000002
【特許請求の範囲】
【請求項１】
観測対象の細胞内に存在する分子群の生化学反応経路の推定方法であって、
膜電位変化を生成する前記経路の数理モデルに対して、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとの少なくとも２種類のパラメータを設定し、誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データと同様な変化を示す前記数理モデルをベイズ推定により抽出し、前記経路を推定することを特徴とする生化学反応経路の推定方法。
【請求項２】
１）誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測する計測ステップと、
２）膜電位変化を生成する生化学反応経路の全候補を数理モデルで表すモデル化ステップと、
３）前記数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定するパラメータ設定ステップと、
４）前記膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、前記事前確率分布とから、前記細胞共通パラメータ及び前記細胞固有パラメータに関する事後確率分布を算出し、前記事前確率分布で指定するパラメータの範囲で前記膜電位変化と同様な変化を示す前記数理モデルを抽出し、前記経路を推定する推定ステップ、
を備えたことを特徴とする請求項１に記載の生化学反応経路の推定方法。
【請求項３】
前記細胞固有パラメータは、
セカンドメッセンジャーの細胞内基質内での拡散時定数、及び、
各イオンチャネルによる膜電位の変化量である、ことを特徴とする請求項１又は２に記載の生化学反応経路の推定方法。
【請求項４】
前記細胞共通パラメータは、
Hill係数、各分子の正規化反応速度、各イオンチャネル密度の反応速度、及び、各イオンチャネルの膜電位への寄与度の平均値である、ことを特徴とする請求項１〜３の何れかに記載の生化学反応経路の推定方法。
【請求項５】
前記経路の全候補の最大数は、
前記分子群の内で経路を構成する分子間相互作用の経路数と、各分子からイオンチャネルへ向かう径路数との合計数をべき数とし、
分子間相互作用の態様数を基数として算出される、ことを特徴とする請求項１〜４の何れかに記載の生化学反応経路の推定方法。
【請求項６】
分子間相互作用の態様は、
相互作用無し、活性作用有り、及び、抑制作用有り、の３態様であることを特徴とする請求項５に記載の生化学反応経路の推定方法。
【請求項７】
前記パラメータ設定ステップにおける事前確率分布は、
膜電位変化へのイオンチャネルの作用を抑制するか否かの制御パラメータと、
分子群の何れかの分子の作用を抑制するか否かの制御パラメータとの少なくとも何れかが、
更に含まれることを特徴とする請求項１〜６の何れかに記載の生化学反応経路の推定方法。
【請求項８】
請求項１〜７の何れかの生化学反応経路の推定方法の各ステップを、コンピュータに実行させるためのプログラム。
【請求項９】
観測対象の細胞内に存在する分子群の生化学反応経路の推定装置であって、
１）誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測する膜電位変化時系列データ計測部と、
２）膜電位変化を生成する前記経路の全候補を数理モデルで表し、前記数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、前記膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、前記事前確率分布とから、前記細胞共通パラメータ及び前記細胞固有パラメータに関する事後確率分布を算出し、前記事前確率分布で指定するパラメータの範囲で前記膜電位変化と同様な変化を示す前記数理モデルを抽出し、前記分子群の生化学反応経路を推定する生化学反応経路推定部、
を備えたことを特徴とする生化学反応経路の推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、細胞内の生化学反応経路の推定に関するものであり、特に、非直接的観測データである膜電位の時系列データを用いた生化学反応経路の推定に関するものである。
【背景技術】
【０００２】
細胞の様々な機能は、細胞内の様々な分子群からなる生化学反応によって実現される。細胞の機能解明のためには、細胞内分子群の生化学反応経路を知ることが重要である。細胞内の生化学反応は試験管内の反応と大きく異なるため、細胞内という環境を維持した状態で反応経路を明らかにしなければならない。しかしながら、現在の生物学の観測技術では、細胞内分子群の生化学反応経路を直接観測することは非常に困難である。
【０００３】
従来、細胞内に存在する分子を同定し、分子と分子が結合する割合を計測して、細胞内分子群の生化学反応経路を推定する方法として、主に次の２つの方法が知られている。１つ目の従来方法は、図３６（１）に示すように、細胞１内の標的分子に抗体を結合させた複数の細胞をつぶし、電気泳動により分子間の結合割合を調べるものである（非特許文献１を参照）。この場合、細胞内の環境で、かつ細胞が生きている状態における分子間の親和性を計測することができないといった問題がある。また、細胞が生きている状態で時間変化する分子間相互反応を観測することも困難である。
２つ目の従来方法は、図３６（２）に示すように、細胞１内の標的分子に蛍光タンパク（プローブ）を発現させ、細胞内の分子間相互反応を観測するものである（非特許文献２，３を参照）。この場合、分子に結合させる蛍光タンパクの種類が最大で３種程度であるため、１〜２種類の分子間相互反応しか観測することができないといった問題がある。
【０００４】
一方、神経細胞間の情報伝達は、軸索と樹状突起で構築されるシナプスを介して行われるが、神経系細胞の発達段階において、神経ネットワークを確立するために軸索誘導は不可欠である。軸索誘導は、神経細胞から軸策が誘導因子に従って伸長して標的細胞に向かう。神経突起の先端にある成長円錐は、細胞外の様々なシグナルに対して応答することにより、神経線維を正確な標的細胞まで導くことが知られている。成長円錐の誘導方向の決定の主な要因として誘導因子の存在がある。成長円錐は、細胞外のシグナルを細胞膜上の受容体により受取り、軸索が標的細胞へと誘導される。また、誘導因子のセカンドメッセンジャーの濃度によって誘導方向がスイッチングする現象が知られ、細胞内の分子群による分子ネットワーク（以下、分子システム）と誘導因子によるシグナルとの関連性など成長円錐の誘導メカニズムの解明が要望されている。西山らは、成長円錐の膜電位を人工的に変化させることにより、成長円錐の誘導方向を決定できることを示し（非特許文献４を参照）、成長円錐の膜電位が膜電位依存のイオンチャネルを制御して成長円錐の誘導方向を決定するといったモデルを提唱している。
【０００５】
しかしながら、成長円錐の膜電位を制御する細胞内分子システムについては部分的にしか解明されていないのが実情である。神経突起先端の成長円錐のように、構造的に小さく、応力に弱い細胞では、実験が困難であり、限られた観測データに基づいて、細胞内分子システムの分子間相互作用のメカニズムを解明していかなければならない。
【０００６】
本発明者の一人である山田は、西山らが示した成長円錐の生物学的知見を用いて、成長円錐の膜電位の時系列データを用いて、分子システムの同定を試みている（非特許文献５を参照）。非特許文献５では、誘導因子(Semaphorin 3A；Sema 3A)のセカンドメッセンジャー(細胞外のシグナルを細胞内へ伝達する際に変換される分子)であるcGMPの細胞内濃度を人工的に操作した条件下のシグナルによる膜電位変化の時系列データを用い、西山らが提唱する分子システムにおいて、未だ明らかでない分子間相互作用の解明を行っている。具体的には、西山らが計測した膜電位変化の時系列データに対し、生化学反応方程式を用いて数理モデル化する。そして、同一条件においても、細胞により膜電位変化の時系列データがばらつくことから、このデータのばらつきを数理モデルのモデルパラメータのばらつきとして扱い、ベイズ推定により、分子システムの同定を行っている。
【０００７】
山田が行った分子システムの同定では、モデルパラメータとして、cGMPの拡散時定数、カリウムチャネルの膜電位寄与度、Hill係数、イオンチャネル濃度電位変換係数、定常速度係数、解離定数などを用い、これらのモデルパラメータを一括りとして全ての細胞に共通するモデルパラメータとして扱っていた。また、数理モデルの候補については、１〜３種類のモデルパラメータを用いた限定されたモデル候補から尤もらしい分子システムを示す数理モデルを選択していた。このため、モデル候補に挙がっていない数理モデルや、細胞毎に存在する個性を考慮できず、細胞内分子システムの分子間相互作用のメカニズムを解明するには不十分であった。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】Shingo Takeuchi, et al., Eph/ephrin reverse signaling induces axonal retraction through RhoA/ROCK pathway, The Journal of Biochemistry volume 158, pages 245-252 (2015).
【非特許文献２】Sadaiea W, Haradab Y, Matsuda M, and Aoki K, Quantitative In Vivo Fluorescence Cross-Correlation Analyses Highlight the Importance of Competitive Effects in the Regulation of Protein-Protein Interactions, Molecular and Cellular Biology, 34, 3272-3290 (2014).
【非特許文献３】Nina K. Thiede-Stan, et al., Attractive and replusive factors act through multi-subunit receptor complexes to regulate nerve fiber growth, Journal of Cell Science colume 128, pages 2403-2414 (2015).
【非特許文献４】Makoto Nishiyama, et al., Membrane potential shifts caused by diffusible guidance signals direct growth-cone turning, Nature Neuroscience volume 11, pages 762-771 (2008).
【非特許文献５】山田達也，奈良先端科学技術大学院大学  情報科学研究科  情報生命科学専攻  修正論文，NAIST-IS-MT1051204（2012年2月）．
【発明の概要】
【発明が解決しようとする課題】
【０００９】
上述の如く、非特許文献４に開示された従来技術では、モデル候補に挙がっていない数理モデルや、細胞毎に存在する個性を考慮できず、細胞内分子システムの分子間相互作用のメカニズムを解明するには不十分であるという問題がある。
また、細胞内分子システムの生化学反応については部分的にしか解明されておらず、構造的に小さく、応力に弱い細胞では、限られた観測データに基づいて、分子間相互作用のメカニズムを解明していく必要がある。
【００１０】
そこで、本発明は、膜電位変化の時系列データといった間接情報を用いて、細胞が生きている状態で発生している分子間相互反応の推定、すなわち、直接観測が困難な細胞内に存在する分子群の生化学反応経路の推定方法及び装置を提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記課題を達成すべく、本発明の生化学反応経路の推定方法は、観測対象の細胞内に存在する分子群の生化学反応経路の推定方法であって、膜電位変化を生成する生化学反応経路の数理モデルに対して、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとの少なくとも２種類のパラメータを設定する。そして、誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データと同様な変化を示す数理モデルをベイズ推定により抽出し、生化学反応経路を推定する。
細胞内の生化学反応経路の結果として表れる膜電位変化の計測は、従来から神経科学の技術として確立されている。本発明では、計測可能な膜電位を用い、膜電位変化の時系列データと同様な変化を示す数理モデルをベイズ推定により抽出して生化学反応経路を推定する。その際、数理モデルに対して、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとの少なくとも２種類のパラメータを設定することで、推定精度を向上させる。
生化学反応の伝搬をトリガーする化学的刺激を細胞に与えると、反応経路の出力として膜電位変化が起きる。この膜電位変化を生成する経路の候補を全て数理モデルで定量的に表現し、観測される膜電位変化と同様の変化を示す数理モデルを抽出する。個々の細胞が持つ個性により、観測される膜電位変化は様々な系列を持つため、抽出した数理モデルにおいて、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータの少なくとも２種類のパラメータを確率分布として表現し、ベイズ推定により最適な反応経路を推定する。
【００１２】
本発明の生化学反応経路の推定方法では、具体的には、下記１）〜４）のステップを備える。
１）誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測する計測ステップ。
２）膜電位変化を生成する生化学反応経路の全候補を数理モデルで表すモデル化ステップ。
３）上記２）の数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定するパラメータ設定ステップ。
４）上記１）の膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、上記３）の事前確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出し、事前確率分布で指定するパラメータの範囲で膜電位変化と同様な変化を示す数理モデルを抽出し、生化学反応経路を推定する推定ステップ。
【００１３】
ここで、細胞固有パラメータは、具体的には、セカンドメッセンジャーの細胞内基質内での拡散時定数、及び、各イオンチャネルによる膜電位の変化量である。
また、細胞共通パラメータは、具体的には、Hill係数、各分子の正規化反応速度、各イオンチャネル密度の反応速度、及び、各イオンチャネルの膜電位への寄与度の平均値である。
【００１４】
本発明の生化学反応経路の推定方法において、生化学反応経路の全候補の最大数は、分子群の内で経路を構成する分子間相互作用の経路数と、各分子からイオンチャネルへ向かう径路数との合計数をべき数とし、分子間相互作用の態様数を基数として算出される。また、分子間相互作用の態様は、具体的には、相互作用無し、活性作用有り、抑制作用有り、の３態様である。
例えば、細胞内分子システムに２つの分子が存在して相互作用がある場合は、分子群の内で経路を構成する分子間相互作用の経路数が２となる。また、細胞内に２つのイオンチャネルがあるとし、２つの分子からそれぞれ特定のイオンチャネルへ向かう場合には、各分子からイオンチャネルへ向かう径路数が２となる。また、分子間相互作用の態様数が、相互作用無し、活性作用有り、抑制作用有りの３態様の場合は、基数が３となる。その場合、生化学反応経路の全候補の最大数は、３４＝８１となる。詳細については後述の実施例で説明する。
【００１５】
本発明の生化学反応経路の推定方法において、パラメータ設定ステップにおける事前確率分布は、膜電位変化へのイオンチャネルの作用を抑制するか否かの制御パラメータと、分子群の何れかの分子の作用を抑制するか否かの制御パラメータとの少なくとも何れかが、更に含まれることでもよい。細胞外の環境に応じて、イオンチャネルの作用、分子の作用を抑制した場合の生化学反応経路を推定できる。
【００１６】
本発明の生化学反応経路推定プログラムは、上記の生化学反応経路の推定方法の各ステップを、コンピュータに実行させるためのプログラムである。
【００１７】
本発明の生化学反応経路の推定装置は、観測対象の細胞内に存在する分子群の生化学反応経路の推定装置であって、下記１）の膜電位変化時系列データ計測部と、下記２）の生化学反応経路推定部を備える。
１）誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測する膜電位変化時系列データ計測部。
２）膜電位変化を生成する生化学反応経路の全候補を数理モデルで表し、数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、パラメータの事前確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出し、事前確率分布で指定するパラメータの範囲で膜電位変化と同様な変化を示す数理モデルを抽出し、分子群の生化学反応経路を推定する生化学反応経路推定部。
【発明の効果】
【００１８】
本発明によれば、観測困難な細胞内の生化学反応を直接観測することなく、計測容易な非直接的観測データである膜電位変化の時系列データといった間接情報を用いて、細胞が生きている状態で、細胞内に存在する分子群の生化学反応経路の推定できるといった効果がある。
【図面の簡単な説明】
【００１９】
【図１】本発明の生化学反応経路の推定方法の概念図
【図２】本発明の生化学反応経路の推定方法の説明図（１）
【図３】本発明の生化学反応経路の推定方法の説明図（２）
【図４】細胞内の分子システムの情報伝達の説明図（１）
【図５】細胞内の分子システムの情報伝達の説明図（２）
【図６】細胞内の分子システムの情報伝達の説明図（３）
【図７】細胞内の分子システムの情報伝達の説明図（４）
【図８】数理モデル候補の説明図
【図９】本発明の生化学反応経路の推定方法のフロー図
【図１０】他の細胞内の分子システムの情報伝達の説明図（１）
【図１１】他の細胞内の分子システムの情報伝達の説明図（２）
【図１２】実施例１の細胞内の分子システムの情報伝達の説明図（１）
【図１３】実施例１の細胞内の分子システムの情報伝達の説明図（２）
【図１４】実施例１の膜電位の計測についての説明図
【図１５】膜電位時系列データの説明図
【図１６】実施例１の細胞内の分子システムの数理モデル候補の説明図
【図１７】cGMPの拡散シミュレーションの説明図
【図１８】膜電位時系列データと数理モデルの尤度についての説明図
【図１９】膜電位時系列データセットの説明図
【図２０】数理モデル候補の尤度の推定結果の説明図（１）
【図２１】数理モデル候補の尤度の推定結果の説明図（２）
【図２２】実施例１の交差検証手順の説明図
【図２３】平均適合度の評価の説明図
【図２４】実施例１の細胞内の分子システムの情報伝達の説明図（３）
【図２５】細胞共通パラメータの説明図
【図２６】実施例２の交差検証手順の説明図
【図２７】実施例２の膜電位の計測についての説明図
【図２８】実施例２の数理モデル候補の尤度の推定結果の説明図（１）
【図２９】実施例２の数理モデル候補の尤度の推定結果の説明図（２）
【図３０】実施例３の予測可能性検証手順の説明図
【図３１】実施例３の数理モデル候補の尤度の推定結果の説明図
【図３２】実施例４の予測可能性検証手順の説明図
【図３３】実施例４の数理モデル候補の尤度の推定結果の説明図
【図３４】実施例５の定常膜電位変化の説明図（１）
【図３５】実施例５の定常膜電位変化の説明図（２）
【図３６】従来の生化学反応経路の推定方法の説明図
【発明を実施するための最良の形態】
【００２０】
本発明の生化学反応経路の推定方法について、図１，２を参照して説明する。図１の概念図に示すように、本発明の生化学反応経路の推定方法では、観測対象の細胞１に対して、細胞膜２上にある受容体３と細胞外の誘導因子（図示せず）が反応し結合すると、二次的に情報伝達物質であるセカンドメッセンジャーが産生される。そして、セカンドメッセンジャーのシグナルにより、細胞内に存在する分子群に細胞周辺の情報を伝達していき、最終的に、細胞膜にある膜貫通タンパク質の一種であるイオンチャネル（７ａ，７ｂ）を介して、情報を隣接細胞に伝える。イオンチャネル（７ａ，７ｂ）は、細胞膜表面に多種類のイオンチャネル分子を備え、これらのイオンチャネル分子が細胞内外のイオンを出し入れする。
【００２１】
本発明の生化学反応経路の推定方法は、上記のように観測対象の細胞内に存在する分子群の情報伝達の生化学反応経路の推定を行うものであり、観測データとしてイオンチャネル（７ａ，７ｂ）の膜電位時系列データを計測し、存在する反応経路を推定する。
図２（１）に示すように、実験系において、生きている細胞の膜電位変化の時系列データを観測する。例えば、塩素（Ｃｌ）イオンチャネル（ＤＣｌ）７ａとナトリウム（Ｎａ）イオンチャネル（ＤＮａ）７ｂの膜電位時系列データを計測する。また、細胞内に存在する分子群の情報伝達経路に関して、候補の反応経路の数理モデルに対して、パラメータを設定し、観測した膜電位時系列データと誤差が少ない反応経路を抽出する。細胞により膜電位変化の時系列データがばらつくことから、このデータのばらつきを数理モデルのモデルパラメータのばらつきとして扱う。すなわち、本発明の生化学反応経路の推定方法は、図３に示すように、候補に挙げた複数の数理モデル（１，２，・・・，ｎ）から、膜電位変化の時系列データと同様な変化を示す数理モデルｉを、ベイズ推定により抽出し、生化学反応経路を推定する。
【００２２】
ここで、細胞内に存在する分子群の情報伝達の生化学反応経路について、図４〜７の細胞内の模式図の一例を参照して説明する。
図４に示すように、観測対象の細胞１に対して、誘導因子を注入（８ａ）し、細胞膜２上にある受容体３と誘導因子を結合させると、細胞内の分子群（図示せず）が活性化され、細胞膜上のイオンチャネル（７ａ，７ｂ）の活性密度が調整され、その結果として、細胞外からのイオンの流入量が変化し膜電位変化を生じる。その膜電位変化の観測１０して、細胞が生きている状態で、細胞内に存在する分子群の生化学反応経路を推定する。
【００２３】
図５に示すように、誘導因子を注入（８ａ）すると、二次的に情報伝達物質であるセカンドメッセンジャー４が産生される（８ｂ）。そして、セカンドメッセンジャー４のシグナルにより、例えば、細胞内に存在する第１／第２の分子（５，６）が活性化され（８ｃ，８ｄ）、イオンチャネル（７ａ，７ｂ）の活性密度、すなわち活性状態の各イオンチャネル分子の濃度が調整される。
【００２４】
図６に示すように、例えば、第１／第２の分子（５，６）からイオンチャネル（７ａ，７ｂ）までに至る分子間の反応経路は未知な部分が多いとする。過去の実験などから、第１の分子５の活性が起因してイオンチャネル７ａの活性密度が調整される反応経路（８ｅ，９ａ）が存在すること、及び、第２の分子６の活性が起因してイオンチャネル７ｂの活性密度が調整される反応経路（８ｆ，９ｂ）が存在すること、これらが既知としても、それ以外の分子間相互作用については未知であるとする。このような場合に、本発明の推定方法は、未知の分子間相互作用の反応経路を推定できる。すなわち、図７に示すように、本発明の推定方法は、例えば、第１／第２の分子（５，６）からイオンチャネル（７ａ，７ｂ）までの未知の分子間相互作用の反応経路（１１ａ，１１ｂ，１２，１３）の有無を推定できる。
【００２５】
本発明の推定方法では、第１／第２の分子（５，６）からイオンチャネル（７ａ，７ｂ）までの未知の分子間相互作用の反応経路（１１ａ，１１ｂ，１２，１３）について、全ての候補を数理モデルとして表現する。分子間相互作用として、相互作用無し（No interaction）、活性作用（Activation）、活性抑制作用（inhibition）の３通りがあるとする。第１の分子５から第２の分子６への反応経路１１ａ、第２の分子６から第１の分子５への反応経路１１ｂ、第１の分子５からイオンチャネル７ｂへの反応経路１２、及び、第２の分子６からイオンチャネル７ａへの反応経路１３の４通りがあり、それぞれの反応経路において、相互作用無し（No interaction）、活性作用（Activation）、活性抑制作用（inhibition）の３通りがある。そのため、反応経路（１１ａ，１１ｂ，１２，１３）の全ての候補は、３４＝８１通り存在することになる。
【００２６】
図８に示すように、８１通りの反応経路（１１ａ，１１ｂ，１２，１３）をマトリックス状に整理する。図８において、右上に示す通り、マトリックス内に示す記号“−”，“→”，“⇒”は、それぞれ相互作用無し（No interaction）、活性作用（Activation）、活性抑制作用（inhibition）を意味する（矢印の向きは反応経路の向きに合せて示す）。
図８のマトリックスの上辺は、第１の分子５から第２の分子６への反応経路１１ａを示し、左辺は、第２の分子６から第１の分子５への反応経路１１ｂを示し、右辺は、第１の分子５からイオンチャネル７ｂへの反応経路１２を示し、下辺は、第２の分子６からイオンチャネル７ａへの反応経路１３を示している。例えば、Ｍ１は、反応経路１１ａ，１１ｂ，１２，１３の全てが相互作用無しであることを意味し、Ｍ９は、反応経路１１ａ，１３の２つが活性抑制作用として働き，反応経路１１ｂ，１２は相互作用無しであることを意味している。同様に、Ｍ７３は、反応経路１１ｂ，１２の２つが活性抑制作用として働き，反応経路１１ａ，１３は相互作用無しであることを意味し、Ｍ８１は、反応経路１１ａ，１１ｂ，１２，１３の全てが活性抑制作用として働くことを意味している。例えば、Ｍ９の下のＭ１８（図示せず）では、反応経路１１ｂが活性作用として働き、反応経路１１ａ，１３の２つが活性抑制作用として働き、反応経路１２は相互作用無しであることを意味している。
【００２７】
本発明の推定方法では、上述の未知の反応経路（１１ａ，１１ｂ，１２，１３）の全ての候補を数理モデルとして表現し、誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測結果と突き合わせ、８１通りの候補の数理モデルに対する尤もらしさから、膜電位変化と同様な変化を示す数理モデルを抽出し、生化学反応経路を推定する。
本発明の推定方法では、数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定する。本発明の推定方法では、細胞内の分子群の相互反応は、個々の細胞に依存せず全ての細胞で共通しているプロセスと、個々の細胞の個性に依存しているプロセスが存在することを前提とする。細胞で共通しているプロセスを数理モデルに表現するものが、個々の細胞に依存しない細胞共通パラメータであり、細胞の個性に依存しているプロセスを数理モデルに表現するものが、個々の細胞に依存する細胞固有パラメータである。
図７に細胞１内には、細胞で共通しているプロセスと、細胞の個性に依存しているプロセスとが共存している概念を符号１ａ，１ｂで示している。
【００２８】
本発明の推定方法では、数理モデルにおいて、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、その事前確率分布と、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出して、８１通りの候補の数理モデルから、膜電位変化と同様な変化を示す数理モデルを抽出する。生化学反応経路を推定する。
【００２９】
本発明の推定方法のフローを図９に示す。観測対象の細胞に誘導因子を注入し（Ｓ０１）、注入した誘導因子のセカンドメッセンジャーの濃度変化によって誘導方向が変わり、細胞内の分子群の相互作用によって、膜電位変化が生じる（Ｓ０２）。公知の計測方法を用いて、膜電位時系列データを計測する（計測ステップ；Ｓ０３）。そして、膜電位変化を生成する生化学反応経路の全候補を数理モデルで表現する（モデル化ステップ；Ｓ０４）。表現した数理モデルにおいて、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定する（パラメータ設定ステップ；Ｓ０５）。なお、モデル化ステップ（Ｓ０４）、パラメータ設定ステップ（Ｓ０５）は、計測ステップ（Ｓ０３）の以前、又は、併行して行っても構わない。
そして、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、個細胞共通パラメータと細胞固有パラメータとに関する事前確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出し、計測した膜電位変化の時系列データを用いて、膜電位変化と同様な変化を示す数理モデルを抽出し、反応経路を推定する（推定ステップ；Ｓ０６）。
【００３０】
上述の説明では、誘導因子のセカンドメッセンジャーのシグナルによって生じる反応経路、すなわち、第１／第２の分子（５，６）からイオンチャネル（７ａ，７ｂ）までの未知の分子間相互作用の反応経路（１１ａ，１１ｂ，１２，１３）について、数理モデルが８１通りの存在することを示した。
本発明の推定方法は、様々なケースに応用できる。例えば、図１０に示すように、誘導因子のセカンドメッセンジャー４のシグナルにより、３つの分子（１５〜１７）が活性化され（８ｃ，８ｄ，８ｅ）、３つのイオンチャネル（７ａ，７ｂ，７ｃ）の活性状態の各イオンチャネル分子の濃度が調整され、膜電位変化が生じるとする。上述と同様に、分子間相互作用として、相互作用無し（No interaction）、活性作用（Activation）、活性抑制作用（inhibition）の３通りがあるとする。
【００３１】
図１１に示すように、３つの分子（１５〜１７）の相互反応経路が６つ、３つの分子（１５〜１７）それぞれから３つのイオンチャネル（７ａ，７ｂ，７ｃ）への反応経路が９つ存在する。そのため、反応経路の全ての候補は、３１５通り存在することになる。但し、既に何れかの反応経路が確認できており、既知の反応経路（図７の反応経路９ａ，９ｂのように既知の経路）が存在するのであれば、候補数を減らせることができる。また、イオンチャネルは、透過するイオンの選択性により、ナトリウム（Ｎａ）イオンチャネル、塩素（Ｃｌ）イオンチャネル、カルシウム（Ｃａ）イオンチャネル、カリウム（Ｋ）イオンチャネルが存在するが、誘導因子のセカンドメッセンジャーのシグナルによって膜電位変化が生じるイオンチャネルに絞りこんで反応経路の候補を決定できる。
【００３２】
以下、本発明の実施形態の一例を、図面を参照しながら詳細に説明していく。なお、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更及び変形が可能である。
【実施例１】
【００３３】
上述したように、西山らは、成長円錐の膜電位を人工的に変化させることにより、成長円錐の誘導方向を決定できることを示し、成長円錐の膜電位が膜電位依存のイオンチャネルを制御して成長円錐の誘導方向を決定するといった細胞内分子システムのモデルを提唱している。
本実施例では、西山らが提唱する分子システムにおいて、未知の反応経路の推定を行った結果を示す。
【００３４】
先ず、誘導因子Semaphorin 3A（以下、Sema 3A）のセカンドメッセンジャーであるcGMPの細胞内濃度を人工的に操作した。実際には、cGMPと同等の効果を持つアナログである8-Br-cGMPを用いた。図１２に示すように、誘導因子のSema 3Aは、成長円錐の細胞膜２上の受容体３に受容されると、cGMPへとシグナル変換される。cGMPはリン酸基を有するエネルギー運搬体であり、PKG（protein kinase G）やCNGC（cyclic nucleotide gated ion channel）の分子と結合する（８ｃ，８ｄ）。PKGやCNGCの分子が直接的または間接的に塩素イオンチャネル（ClC）７ａ及びナトリウムイオンチャネル（NaC）７ｂへ作用し、膜電位変化を生じさせる。ここで、過去の西山らの実験によって、PKGはナトリウムイオンチャネル（NaC）７ｂを介して膜電位の脱分極（静止膜電位から膜電位を上げる働き）を引き起こし、また、CNGCは塩素イオンチャネル（ClC）７ａを介して膜電位の過分極（静止膜電位から膜電位を下げる働き）を引き起こすことが既知である。しかし、それ以外の分子間相互作用については判明しておらず、図１３に示すような反応経路（９ａ，９ｂ）の存在は既知であるが、反応経路（１１ａ，１１ｂ，１２，１３）が存在するか否かは未知である。
本実施例では、図１３に示す反応経路（１１ａ，１１ｂ，１２，１３）の存在について、膜電位時系列データを用いて推定する。
【００３５】
図１４に示すように、cGMPのアナログである8-Br-cGMPを、マイクロピペット３０を用いて観測対象の成長円錐細胞３１に注入し、その後の膜電位変化を計測した。実験は生物にとって通常に近い条件下で実施した。8-Br-cGMPの濃度を変えて、幾つかの濃度で実験を行った。注入後１０分程度の膜電位時系列データを計測した。なお、成長円錐の膜電位の計測は、マイクロピペット３０を固定し、Whole Cell Recordingによって計測した。図１５に、8-Br-cGMPの濃度を変えて膜電位変化を計測した様子を示す。図１５のグラフにおいて、横軸は時間（秒）、縦軸は成長円錐膜電位（ｍＶ）を示す。8-Br-cGMPの濃度を変えると、膜電位時系列データは異なる。
【００３６】
図１６は、PKGとCNGCの２つの分子からClCとNaCの２つのイオンチャネルまでの未知の分子間相互作用の反応経路（１１ａ，１１ｂ，１２，１３）について、全ての候補をマトリックスに表現したものである。分子間相互作用として、相互作用無し（No interaction）、活性作用（Activation）、活性抑制作用（inhibition）の３通りがある。CNGCからPKGへの反応経路、PKGからCNGCへの反応経路、CNGCからNaCへの反応経路、及び、PKGからClCへの反応経路の４通りがあり、それぞれの反応経路において、相互作用無し、活性作用、活性抑制作用の３通りがあり、反応経路の全ての候補は、３４＝８１通り存在する。
【００３７】
図１６のマトリックスの上辺は、CNGCからPKGへの反応経路１１ａを示し、左辺は、PKGからCNGCへの反応経路１１ｂを示し、右辺は、CNGCからNaCへの反応経路１２を示し、下辺は、PKGからClCへの反応経路１３を示している。Ｍ１は、反応経路１１ａ，１１ｂ，１２，１３の全てが相互作用無しであることを意味し、Ｍ９は、反応経路１１ａ，１３の２つが活性抑制作用として働き，反応経路１１ｂ，１２は相互作用無しであることを意味し、Ｍ７３は、反応経路１１ｂ，１２の２つが活性抑制作用として働き、反応経路１１ａ，１３は相互作用無しであることを意味し、Ｍ８１は、反応経路１１ａ，１１ｂ，１２，１３の全てが活性抑制作用として働くことを意味している。
【００３８】
未知の反応経路（１１ａ，１１ｂ，１２，１３）の全ての候補を数理モデルとして表現し、膜電位変化の時系列データを計測結果と突き合わせ、８１通りの数理モデルに対する尤もらしさから、膜電位変化と同様な変化を示す数理モデルを抽出する方法について、以下に詳細に説明する。
まず、数理モデルにおいて、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定する。具体的な細胞共通パラメータと細胞固有パラメータを、それぞれ下記表１，２に示す。
【００３９】
【表１】
2019154299.tif 000003
【００４０】
【表２】
2019154299.tif 000004
【００４１】
細胞内の分子群の相互反応において、個々の細胞に依存せず全ての細胞で共通しているプロセスと、個々の細胞の個性に依存しているプロセスが存在するとし、細胞で共通しているプロセスを細胞共通パラメータ（上記表１の１０個のパラメータ）で表現し、細胞の個性に依存しているプロセスを細胞固有パラメータ（上記表２の４個のパラメータ）で表現する。
【００４２】
そして、ベイズ推定を用いて、８１通りの候補の数理モデルから、膜電位変化と同様な変化を示す数理モデルを抽出する。具体的には、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、その事前確率分布と、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出する。以下、具体的な数式を示しながら説明する。
【００４３】
実験では、アフリカツメガエルの脊髄ニューロンの成長円錐細胞を用い、成長円錐細胞の中心付近にcGMPのアナログ（8-Br-cGMP）の溶液で満たされた記録ピペットを取り付け、8-Br-cGMPを成長円錐全体に拡散させた。実験における8-Br-cGMPの刺激濃度は数μＭであり、成長円錐の内在的なcGMP（数１０ｎＭ）との間に大きな濃度差があるため、ピペットから成長円錐細胞内への8-Br-cGMPの流入には拡散過程を考慮した。成長円錐Ｓにおける8-Br-cGMP濃度の時系列は、時間ｔの単純指数関数として近似され、下記式（１）によって拡散過程として表すことができる。
ここで、Ｓｍａｘはピペット中の8-Br-cGMP濃度であり、これが入力濃度刺激の最大値になる。τsは拡散時定数である。 拡散時定数τsは、成長円錐細胞の体積に依存し、上記表２に示す細胞固有パラメータである。
【００４４】
【数１】
2019154299.tif 000005
【００４５】
拡散時定数τsの変化を表現するために、図１７に示す３つの第１〜第３区画（１，２，３）からなるモデル成長円錐における数値計算から推定された確率分布としてモデル化した。内側の２つの第１，２区画（１，２）は１μｍの厚さを有し、外側の第３区画（３）はガウスノイズε〜Ｎ（０，１２）によって修正された１〜２μｍの可変厚さを有する。第ｎ区画の8-Br-cGMP濃度Ｓｎ（μＭ）はフィックの法則に従うべきであり、下記式（２）として表すことができる。
【００４６】
【数２】
2019154299.tif 000006
【００４７】
ここで、Ｄ、Ａn-1,n、ｄn-1,n、Ｖｎは、8-Br-cGMP（１μm２／ｓ）２の拡散係数、ｎ−１番目とｎ番目の区画中点間の距離、および第ｎ区画の容積にそれぞれ対応するものである。上記式（２）の境界条件Ｓ０＝Ｓmaxを用いて数値的に計算し、第２区画（３）の時間経過が上記数式１の指数関数で近似した。種々の成長円錐サイズにおけるモンテカルロシミュレーションによって、拡散時定数τsは、平均が４０秒である非負のガウス分布とした。
【００４８】
上述の如く、反応経路の数理モデルでは、8-Br-cGMPの結合時のCNGCやPKGの活性化の反応経路の存在を考慮した。以下、Ｘ，Ｙとして表されるCNGCおよびPKGの下流因子（ＤＦ）の活性レベルが、8-Br-cGMPの濃度レベルが増加するにつれて増加し、ＤＦの活性レベルは、正常な状態（または飽和）に達するまでの時間で観測した。上記数式ＳのHill型関数として、下記式（３），（４）に示すように、ＸとＹを表現する準定常近似を適用した。ここで、Hill係数ｎ，ｍは、下流プロセスに対して効果的なトリガーを行う。 Ｘ及びＹは、[０，１]の範囲で正規化され、実際の活動レベルはその最大値を掛けて表現される。
【００４９】
【数３】
2019154299.tif 000007
【００５０】
上記式（３），（４）における解離定数ＫｘおよびＫｙは、細胞内の8-Br-cGMP濃度とＤＦ活性との間の反応の順方向速度（ＫｘｂおよびＫｙｂ）に逆方向速度（ＫｘｆおよびＫｙｆ）の比によって定義される。 順方向速度および逆方向速度に直線近似を適用すると、有効解離定数は、下記式（５），（６）で表すことができる。ここで、分母と分子はそれぞれ順反応速度と逆反応速度を表す。 数式中の８つのパラメータの値は、細胞共通パラメータである。
【００５１】
【数４】
2019154299.tif 000008
【００５２】
膜電位（MP）変化は、イオンチャネル分子の密度に依存する。このため、数理モデルは、CNGCおよびPKGによるClCおよびNaCのイオンチャネル分子の密度の調節を考慮している。イオンチャネル分子の密度のダイナミクスは、それぞれＤＣｌ、ＤＮａとして表されるイオンチャネル密度の比の変化によって定義される。ＤＣｌ、ＤＮａの調整はＸやＹに依存するため、細胞共通パラメータは、数理モデルにおける信号処理の速度（微小時間スケール）を制限することになる。 したがって、生化学反応のダイナミクスは、常微分方程式（ODE：Ordinary Differential  Equation）によって記述され、下記数式で表すことができる。ここで、添え字付きのｋは、正規化された変数Ｘ，Ｙの定数または関数のいずれかであり、有効な順方向反応速度または逆方向反応速度である。
【００５３】
【数５】
2019154299.tif 000009
【００５４】
線形近似によって、生化学反応のダイナミクスは、下記式（９）〜（１２）として表すことができる。下記式において、反応速度ｋ**0は一定であり（ＸまたはＹとは無関係）、他の反応速度はＸまたはＹによって与えられる最大速度である。下記式（９）〜（１２）の全てのパラメータは、細胞共通パラメータであり、上記式（５）および（６）のパラメータを用いて、数理モデルを表す。
【００５５】
【数６】
2019154299.tif 000010
【００５６】
上記式（９）〜（１２）において、イオンチャネルからCNGCおよびPKGへのフィードバックの反応経路を仮定しなかった。いくつかの状態を伴うフィードバック反応経路は、発散、振動、双安定性などのシステム状態の不安定さをもたらすが、そのような不安定な応答は、実験では膜電位時系列データでは観察されなかった。
【００５７】
イオンチャネル密度（ＤＣｌおよびＤＮａ）が膜電位（MP）を調整し、イオンチャネル密度からMPへの信号変換については、HodgkinおよびHuxleyの式によって十分に確立されていることから、膜電位は下記式（１３）の線形表現として表すことができる。
【００５８】
【数７】
2019154299.tif 000011
【００５９】
ここで、Ｃｍはイオンチャネル容量と膜電位に依存する膜容量、Ｖは膜電位、ＶＮａ、ＶＫ、ＶＣｌは反転電位、ｇＮａ、ｇＫ、ｇＣｌはチャネルコンダクタンスである。静止状態でのカリウムコンダクタンスｇＫは、支配的な透過性（ｇＫ>>ｇＣｌ、ｇＮａ）を有する。上記式（１３）をゼロに設定することによって、定常状態および静止状態での膜電位（MP）は、下記式（１４）で表すことができる。
【００６０】
【数８】
2019154299.tif 000012
【００６１】
ｇＮａ、ｇＫ、ｇＣｌ（上付きバー）は、定常状態におけるチャネルコンダクタンスである。数理モデルでは、ＤＣｌおよびＤＮａは8-Br-cGMPの下流因子によってイオンチャネルの分子密度が調節される。 静止電位における膜電位Ｖの小さな変化のため、8-Br-cGMPのシグナルによる膜電位変化は、以下の直線近似式に従うことになる。
【００６２】
【数９】
2019154299.tif 000013
【００６３】
ＤＣｌおよびＤＮａは、それぞれ8-Br-cGMPのシグナルによって誘導されたClCおよびNaCの有効濃度の正規化された変化であり、それらの変化が最大に達した場合には1つになる。また、パラメータＶＫ、ＡＣｌおよびＡＮａは個々の細胞特性に依存し、その値は各細胞に特異的であり、上記表２に示す細胞固有パラメータである。
【００６４】
成長円錐体積およびイオンチャネル密度のパラメータが膜電位変化に影響を及ぼすことから、ベイズ推定を用いる。数理モデルのパラメータは、細胞共通パラメータと細胞固有パラメータ  と実験環境として強制的に導入された実験条件パラメータの３つに分類される。
生化学反応速度、Hill係数およびＡＣｌおよびＡＮａの平均（ＡＣｌ０およびＡＮａ０）を含む細胞共通パラメータセット（θ）は、すべての細胞に共通であるため、全ての膜電位時系列データセットから推定した。一方、個々の細胞の特性（例えば、成長円錐体積、チャネル密度）に大きく依存する細胞固有パラメータセット（φ）＝ {ＶＫ、ＡＣｌ、ＡＮａ、τＳ}を各膜電位時系列データについて推定した。実験条件パラメータセット（ｃ）＝ {ηＣｌ、ηＮａ、Ｓｍａｘ}は、薬理学的条件および適用された8-Br-cGMP濃度などの実験的に誘導されたパラメータである。
【００６５】
ベイズ推定のアプローチにより、事前分布（上記表１に示すように、殆どガウス分布の非負部分）によって与えられた制約下で、計測した膜電位時系列データからモデルパラメータを推定した。モデルパラメータの事後分布を下記式（１６）で表した。
【００６６】
【数１０】
2019154299.tif 000014
【００６７】
ここで、実験条件パラメータセットｃは、見やすさのために省略されている。エビデンスＥ（Ｖ，Ｍ）は、膜電位時系列データセットＶおよび数理モデルＭによって決定されるモデル評価基準である。事前分布のパラメータｐ（φ、θ｜Ｍ）は、上述の表１，２に示したものである。与えられたデータセットに対する適応レベルを表す尤度は、ガウス分布の積によって与えられる。
【００６８】
【数１１】
2019154299.tif 000015
【００６９】
ΔＶｉ(ｔ)は、観測された膜電位とＶｉ(ｔ)との間の誤差を示す。ここで、iは膜電位時系列データのインデックスであり、Ｉ＝１６、Ｔｉ＝３００〜９００である。ｉ番目の膜電位時系列データのσiはデータ前処理中に計算された観測ノイズの大きさである。
なお、膜電位の観測データは、細胞の自発的興奮等の影響により、観測ノイズ以外の予期しないトレンドが生じるため、これらのノイズをデータの前処理によりデータ補正したものを推定に用いる時系列データとしている。具体的には、明らかにはずれ値であるとわかる区間については、範囲を指定し、線形補間によりデータを内挿した。
【００７０】
次に、推定結果について説明する。
数理モデル間で、計測した膜電位時系列データの尤度を比較し、上述の８１通りの候補の数理モデルのそれぞれのモデル適合度を調べ、モデル適合度を評価してモデル抽出を行った。事前確率分布によって制約された細胞共通パラメータ（θ）と細胞固有パラメータ（φ）を用いて数理モデルの膜電位時系列データを計算した。上述の表１，２に示したように、パラメータの事前確率分布のほとんどにゼロ平均のガウス分布の非負部分を使用した。
【００７１】
そして、単一の実験的膜電位時系列データと数理モデルの膜電位時系列データとの間の誤差を尤度関数として比較した。実験的膜電位時系列データに細胞依存性ノイズが存在するため、各ノイズレベルで誤差を標準化した。サンプリングされた膜電位時系列データと同じ膜電位時系列データの平滑化された時間平均との間の誤差からノイズレベル（σi）を推定した。
このプロセスは、すべての膜電位時系列データの全データ点について反復され、すべてのデータセットのモデル適合度として全尤度を得るために全ての尤度が乗算された。パラメータ値に依存するモデル適合度は、実験的膜電位時系列データと数理モデルの膜電位時系列データとの間の類似性のみを提供し、パラメータの妥当性を無視する。したがって、数理モデルがデータセットにあてはまるときにパラメータが不合理な値を取る可能性があるため、モデルの適合度だけでは、もっともらしい数理モデルを選択するには不十分である。
【００７２】
上述の如く、ベイズ推定を使用して、データセットのオーバーフィットからパラメータを回避し、パラメータの妥当性を組み込むことによって、数理モデルの妥当性を計算した。数理モデルの妥当性は、モデル適合度としての尤度とパラメータの妥当性としての事前確率との積として表される。
ここで、尤度として表されるモデル適合度は、図１８（４）に示す尤度（likelihood）グラフに示すように定義される。定数パラメータを無視して上記式（１８）の自然対数をとることによって得られる。エビデンスの対数を数理モデルの妥当性の指標とし、エビデンスの実際の計算は、マルコフ連鎖モンテカルロ法（MCMC）による事後分布からのサンプリングにより実行した。
【００７３】
実験的な膜電位時系列データセットの数は１６であり、その内訳は、１０μＭ 8-Br-cGMPによる時系列データ数ｎ＝７、膜電位変化へのClCの寄与をブロッキングするDNDSを実験条件として加えた時系列データ数ｎ＝５、膜電位変化へのNaCの寄与をブロッキングするSTXを実験条件として加えた時系列データ数ｎ＝４である。図１９に示すように、全て（１６セット）の実験的な膜電位時系列データセットから得られる細胞共通パラメータをθallとして表す。また、細胞に依存する細胞固有パラメータ（φi）は、モデル適応度を計算するための全ての膜電位時系列データセットからｉ番目の膜電位時系列データからのみ推定されたものを表す。１０μＭ 8-Br-cGMPから誘発される膜電位時系列データセットから細胞内の分子システム内のシグナル伝達ネットワークを同定することを目指して、検証、特異性、および予測可能性のための全データセットから推定された細胞共通パラメータ（θall）を利用したテストを行った。
【００７４】
図２０に示すように、各々の数理モデルの実験的膜電位時系列データセットから計算されたエビデンスは、より小さなエビデンス（グレーの薄い領域）を有するグループａと、より大きなエビデンスを有するグループｂ（グレーの濃い領域）を有していた。
【００７５】
図２１に示すように、各グループ内の数理モデル間のエビデンスの違いは、ａ，ｂの２つのグループ間のものよりもずっと小さく、グループｂの数理モデルに一般的に存在する反応経路が重要な役割を果たす可能性があることを示唆している。 注目すべきことは、グループｂ内に存在する最も単純な分子相互作用を表すＭ７は、膜電位変化中にPKGを介してClCを抑制するcGMPのシグナル伝達に関与することを明らかにしていることである。
【００７６】
細胞共通パラメータがどの程度安定しているかを評価するために、交差検証（CV；Cross-validation）を実行することによって、推定された細胞共通パラメータの特異性を確認した。図２２に示すように、最初に１６個のデータセットを２つのサブセットに分割する。１つのサブセットは、１５のデータセットであり、細胞共通パラメータ（θiMAP、i番目の膜電位時系列データを除く）を訓練する。もう１つのサブセットは、膜電位時系列データセットを１つ（抜き出したデータセット）だけ有し、細胞固有パラメータ（φi; i番目の膜電位時系列データセット）を訓練する。抜き出した膜電位時系列データセットは、訓練されていないデータとしての可能性を計算するために使用される。これは、細胞共通パラメータに組み込まれていないためである。全データセットサンプル（ｉ＝１，２，・・・，１６）について１６回訓練と計算手順を繰り返し、モデル適合度を得るために１６の尤度を乗算した。
【００７７】
細胞共通パラメータの推定値が安定している場合には、Ｍ７の細胞共通パラメータは交差検証によって大きく変更されるべきではなく、訓練されていないデータセットに対するモデル適応度のマトリックス表現は訓練データセットのそれと同様でなければならない。
図２３（１）と（２）を比較すると、データポイントあたりの対数尤度で表されるモデル適合度の平均値（θiMAPとの平均適合度）は、グループａでは小さく、グループｂでは大きく見える。
【００７８】
交差検証において、各々の訓練ごとに細胞共通パラメータ（θiMAP; ｉ＝１，２，・・・，１６）を推定した。図２５に示すように、Ｍ７の細胞共通パラメータ（θallMAP）がθiMAPの±１の標準偏差（ｓ．ｄ．）内に入っている。これらの結果は、本推定方法が、交差検証のデータセットの組み合わせによって著しく影響を受けることなく、細胞共通パラメータが適度に安定していることを示している。
【００７９】
以上述べたように、誘導因子Sema 3AのセカンドメッセンジャーであるcGMPのアナログである8-Br-cGMPを用いて、成長円錐の細胞に刺激を与えることにより、塩素イオンチャネルとナトリウムイオンチャネルの動態により変化する膜電位時系列データを計測し、計測した時系列データ用いて、本発明の推定方法を適用し、未知の生化学反応経路を推定できた。
【実施例２】
【００８０】
図２６は、交差検証の手順を示している。図２６に示す検証手順は、上述の図２２に示す検証手順と異なり、エビデンス評価に全く用いていないデータを用いて細胞固有パラメータφを推定した。最初に１６個のデータセットを用いて、細胞共通パラメータ（θallMAP）を訓練し、その後、エビデンス評価に全く用いていないｎ個のデータセットを用いて、細胞固有パラメータ（φi; i番目の膜電位時系列データセット）を訓練し、候補の数理モデルを評価した。
【００８１】
本実施例では、下記表３の試薬を用いて、特定の分子の働きを抑制した条件下で膜電位時系列データを計測した。図２７に、膜電位時系列データ計測の実験手順を示す。図２７に示すように、時点Ａにおいてチャネルブロッカー（STX／DNDS）やPKG抑制薬（KT5823）を導入し、十分な時間（３０分以上）経過後の時点Ｂにおいて、cGMPのアナログである8-Br-cGMPを導入すると同時に膜電位の観測を行った。
【００８２】
【表３】
2019154299.tif 000016
【００８３】
図２８（１）は、１０μＭ 8-Br-cGMPのデータセットを用い、チャネルブロッカー（STX／DNDS）を導入した条件下で、チャネルブロッカーの効果（ηＣｌ、ηＮａ）の値（存在／非存在）を入れ替えたデータセットを用いた検証結果、図２８（２）は、５μＭ 8-Br-cGMPのデータセットを用い、刺激を弱め、（１）と同様に、チャネルブロッカー（STX／DNDS）の導入条件下で、チャネルブロッカーの効果（ηＣｌ、ηＮａ）の値を入れ替えたデータセットを用いた検証結果を示している。図２８に示すように、共にcGMPによって駆動される細胞内分子システムが機能していることが確認できた。
【００８４】
一方、図２９（１）は、１０μＭ 8-Br-cGMPのデータセットを用い、PKG抑制薬の導入条件下でのデータセットを用いた検証結果、図２９（２）は、Netrin（脊髄交連神経細胞の軸索を誘引する分泌タンパク質として同定された軸索ガイダンス分子）の環境下でのデータセットを用いた検証結果を示している。共に、cGMPによって駆動される本来の分子システムとは異なるシステムが機能していることが確認できた。
【実施例３】
【００８５】
図３０は、予測可能性検証の手順を示している。１０μＭ 8-Br-cGMPの１６個のデータセットを２つのサブセットに分割し、１つのサブセット（１５のデータセット）を用いて、細胞共通パラメータ（θiMAP、i番目の膜電位時系列データを除く）を訓練する。もう１つのサブセット（i番目の膜電位時系列データ）において、時系列の過渡期（０〜２５０秒）の膜電位時系列データを用いて細胞固有パラメータφを推定し、推定に用いていない残りの膜電位時系列データ（２５０秒〜）に対する予測誤差を尤度により評価した。すなわち、図３１（１）に示すグラフにおいて、グレイ領域が時系列の過渡期（０〜２５０秒）の領域であり、この領域内の膜電位時系列データを用いてパラメータφを推定し、残りの領域の膜電位時系列データを尤度評価に用いた。図３１（１）のグラフにおいて、黒線は実験データを表し、グレイの実線および破線はそれぞれ、Ｍ７，Ｍ１による予測した膜電位時系列データを表す。
図３１（２）は、各候補の数理モデルに対し評価した尤度を示している。図３１（２）から、cGMPによって駆動される細胞内分子システムが機能していると言える。
【実施例４】
【００８６】
図３２は、予測可能性検証のその他の手順を示している。１０μＭ 8-Br-cGMPの１６個のデータセットを用いて、細胞共通パラメータ（θallMAP）を訓練する。そして、エビデンス評価に全く用いていない１１個のデータセットであって、上述の実施例と同様に、時系列の過渡期（０〜２５０秒）の膜電位時系列データを用いて細胞固有パラメータφを推定する。エビデンス評価に全く用いていない１１個のデータセットは、５μＭ 8-Br-cGMPにより刺激を弱めたデータである。推定に用いていない残りの膜電位時系列データ（２５０秒〜）に対する予測誤差を尤度により評価した。すなわち、図３３（１）に示すグラフにおいて、グレイ領域が時系列の過渡期（０〜２５０秒）の領域であり、この領域内の膜電位時系列データを用いてパラメータφを推定し、残りの領域の膜電位時系列データを尤度評価に用いた。図３３（１）のグラフにおいて、黒線は実験データを表し、グレイの実線および破線はそれぞれ、Ｍ７，Ｍ１による予測した膜電位時系列データを表す。
図３３（２）は、各候補の数理モデルに対し評価した尤度を示している。図３３（２）に示すように、cGMPによって駆動される細胞内分子システムが機能していることが確認できた。
【実施例５】
【００８７】
膜電位は、定常状態では、ほぼ一定の値を維持する。これは、8-Br-cGMP刺激のレベルに双方向の依存性を有するからである。実施例１における図１６のマトリックスのカラムＭ７において強調されたClCのPKG抑制を含む細胞内分子システムが、cGMP依存性の双方向の膜電位変化の誘発に必要か否かについて確認した。異なる３通りの濃度のcGMPと、１０μＭの8-Br-cGMP（ｎ＝１６）から推定されるθallMAPを用いてＭ７を計算した。
モデルシミュレーションは、（１）低濃度（０．５μＭ）の8-Br-cGMPによって刺激されて生じる過分極、（２）適度な8-Br-cGMPによって刺激され徐々に回復するシャープな過分極、（３）高濃度（２０μＭ）の8-Br-cGMPによって刺激され、脱分極に変換されたシャープな過分極の３通り実施した。
【００８８】
次に、モデルの定常状態の膜電位（図３４の長時間経過後の円点を参照）と実験データを比較した。図３４は、数理モデルＭ７による各cGMP濃度に対する膜電位時系列データの一例を示している。図３５は、数理モデルＭ７によるcGMP濃度依存性の膜電位変化と実験データの比較を示している。図中、●印は、図３４と対応している。シミュレートされたcGMP濃度依存性の膜電位変化は、双方向の膜電位変化を示す（図３５の点線のグラフを参照）。
モデルシミュレーションは双方向の膜電位変化を示したが、脱分極の大きさははるかに大きく、脱分極の発生は実験的に観測されたものよりはるかに低い8-Br-cGMP濃度で現れた。この差は、実験での8-Br-cGMPの浴適用（bath）と比較して、記録ピペットを通した8-Br-cGMPによる直接刺激との差から生じるものである。 
8-Br-cGMPの膜浸透の単純なHill様モデルを導入したとき、8-Br-cGMPは、実験的に観測されたものにはるかに近い双方向の膜電位変化を再現することができた（図３５の黒四角を参照）。
【産業上の利用可能性】
【００８９】
本発明は、神経細胞などの細胞内分子群の生化学反応経路の推定に有用である。
【符号の説明】
【００９０】
１  細胞
２  細胞膜
３  受容体
４  セカンドメッセンジャー
５  第１の分子
６  第２の分子
７ａ〜７ｃ  イオンチャネル
８ａ〜８ｆ，９ａ，９ｂ，１１〜１３    反応経路
１０  膜電位変化の観測
３０  マイクロピペット
３１  成長円錐細胞
【図１】
2019154299.tif 000017
【図２】
2019154299.tif 000018
【図３】
2019154299.tif 000019
【図４】
2019154299.tif 000020
【図５】
2019154299.tif 000021
【図６】
2019154299.tif 000022
【図７】
2019154299.tif 000023
【図８】
2019154299.tif 000024
【図９】
2019154299.tif 000025
【図１０】
2019154299.tif 000026
【図１１】
2019154299.tif 000027
【図１２】
2019154299.tif 000028
【図１３】
2019154299.tif 000029
【図１４】
2019154299.tif 000030
【図１５】
2019154299.tif 000031
【図１６】
2019154299.tif 000032
【図１７】
2019154299.tif 000033
【図１８】
2019154299.tif 000034
【図１９】
2019154299.tif 000035
【図２０】
2019154299.tif 000036
【図２１】
2019154299.tif 000037
【図２２】
2019154299.tif 000038
【図２３】
2019154299.tif 000039
【図２４】
2019154299.tif 000040
【図２５】
2019154299.tif 000041
【図２６】
2019154299.tif 000042
【図２７】
2019154299.tif 000043
【図２８】
2019154299.tif 000044
【図２９】
2019154299.tif 000045
【図３０】
2019154299.tif 000046
【図３１】
2019154299.tif 000047
【図３２】
2019154299.tif 000048
【図３３】
2019154299.tif 000049
【図３４】
2019154299.tif 000050
【図３５】
2019154299.tif 000051
【図３６】
2019154299.tif 000052