TOP > 国内特許検索 > 膜電位時系列を用いた細胞内分子群の生化学反応経路の推定方法及び装置 > 明細書

明細書 :膜電位時系列を用いた細胞内分子群の生化学反応経路の推定方法及び装置

発行国 日本国特許庁(JP)
公報種別 公開特許公報(A)
公開番号 特開2019-154299 (P2019-154299A)
公開日 令和元年9月19日(2019.9.19)
発明の名称または考案の名称 膜電位時系列を用いた細胞内分子群の生化学反応経路の推定方法及び装置
国際特許分類 C12Q   1/02        (2006.01)
C12M   1/34        (2006.01)
FI C12Q 1/02
C12M 1/34 D
請求項の数または発明の数 9
出願形態 OL
全頁数 25
出願番号 特願2018-044724 (P2018-044724)
出願日 平成30年3月12日(2018.3.12)
発明者または考案者 【氏名】作村 諭一
【氏名】山田 達也
【氏名】西山 誠
【氏名】ホン,キョンスウ
出願人 【識別番号】504143441
【氏名又は名称】国立大学法人 奈良先端科学技術大学院大学
個別代理人の代理人 【識別番号】110000822、【氏名又は名称】特許業務法人グローバル知財
審査請求 未請求
テーマコード 4B029
4B063
Fターム 4B029AA07
4B029BB11
4B029CC03
4B029FA12
4B063QA01
4B063QA05
4B063QA18
4B063QQ08
4B063QR42
4B063QS39
4B063QX05
要約 【課題】膜電位変化の時系列データといった間接情報を用いて、細胞が生きている状態下での細胞内に存在する分子群の生化学反応経路の推定方法及び装置を提供する。
【解決手段】誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測し、膜電位変化を生成する生化学反応経路の全候補を数理モデルで表し、数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、事前確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出し、事前確率分布で指定するパラメータの範囲で膜電位変化と同様な変化を示す数理モデルを抽出し、生化学反応経路を推定する。
【選択図】図1
特許請求の範囲 【請求項1】
観測対象の細胞内に存在する分子群の生化学反応経路の推定方法であって、
膜電位変化を生成する前記経路の数理モデルに対して、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとの少なくとも2種類のパラメータを設定し、誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データと同様な変化を示す前記数理モデルをベイズ推定により抽出し、前記経路を推定することを特徴とする生化学反応経路の推定方法。
【請求項2】
1)誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測する計測ステップと、
2)膜電位変化を生成する生化学反応経路の全候補を数理モデルで表すモデル化ステップと、
3)前記数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定するパラメータ設定ステップと、
4)前記膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、前記事前確率分布とから、前記細胞共通パラメータ及び前記細胞固有パラメータに関する事後確率分布を算出し、前記事前確率分布で指定するパラメータの範囲で前記膜電位変化と同様な変化を示す前記数理モデルを抽出し、前記経路を推定する推定ステップ、
を備えたことを特徴とする請求項1に記載の生化学反応経路の推定方法。
【請求項3】
前記細胞固有パラメータは、
セカンドメッセンジャーの細胞内基質内での拡散時定数、及び、
各イオンチャネルによる膜電位の変化量である、ことを特徴とする請求項1又は2に記載の生化学反応経路の推定方法。
【請求項4】
前記細胞共通パラメータは、
Hill係数、各分子の正規化反応速度、各イオンチャネル密度の反応速度、及び、各イオンチャネルの膜電位への寄与度の平均値である、ことを特徴とする請求項1~3の何れかに記載の生化学反応経路の推定方法。
【請求項5】
前記経路の全候補の最大数は、
前記分子群の内で経路を構成する分子間相互作用の経路数と、各分子からイオンチャネルへ向かう径路数との合計数をべき数とし、
分子間相互作用の態様数を基数として算出される、ことを特徴とする請求項1~4の何れかに記載の生化学反応経路の推定方法。
【請求項6】
分子間相互作用の態様は、
相互作用無し、活性作用有り、及び、抑制作用有り、の3態様であることを特徴とする請求項5に記載の生化学反応経路の推定方法。
【請求項7】
前記パラメータ設定ステップにおける事前確率分布は、
膜電位変化へのイオンチャネルの作用を抑制するか否かの制御パラメータと、
分子群の何れかの分子の作用を抑制するか否かの制御パラメータとの少なくとも何れかが、
更に含まれることを特徴とする請求項1~6の何れかに記載の生化学反応経路の推定方法。
【請求項8】
請求項1~7の何れかの生化学反応経路の推定方法の各ステップを、コンピュータに実行させるためのプログラム。
【請求項9】
観測対象の細胞内に存在する分子群の生化学反応経路の推定装置であって、
1)誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測する膜電位変化時系列データ計測部と、
2)膜電位変化を生成する前記経路の全候補を数理モデルで表し、前記数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、前記膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、前記事前確率分布とから、前記細胞共通パラメータ及び前記細胞固有パラメータに関する事後確率分布を算出し、前記事前確率分布で指定するパラメータの範囲で前記膜電位変化と同様な変化を示す前記数理モデルを抽出し、前記分子群の生化学反応経路を推定する生化学反応経路推定部、
を備えたことを特徴とする生化学反応経路の推定装置。
発明の詳細な説明 【技術分野】
【0001】
本発明は、細胞内の生化学反応経路の推定に関するものであり、特に、非直接的観測データである膜電位の時系列データを用いた生化学反応経路の推定に関するものである。
【背景技術】
【0002】
細胞の様々な機能は、細胞内の様々な分子群からなる生化学反応によって実現される。細胞の機能解明のためには、細胞内分子群の生化学反応経路を知ることが重要である。細胞内の生化学反応は試験管内の反応と大きく異なるため、細胞内という環境を維持した状態で反応経路を明らかにしなければならない。しかしながら、現在の生物学の観測技術では、細胞内分子群の生化学反応経路を直接観測することは非常に困難である。
【0003】
従来、細胞内に存在する分子を同定し、分子と分子が結合する割合を計測して、細胞内分子群の生化学反応経路を推定する方法として、主に次の2つの方法が知られている。1つ目の従来方法は、図36(1)に示すように、細胞1内の標的分子に抗体を結合させた複数の細胞をつぶし、電気泳動により分子間の結合割合を調べるものである(非特許文献1を参照)。この場合、細胞内の環境で、かつ細胞が生きている状態における分子間の親和性を計測することができないといった問題がある。また、細胞が生きている状態で時間変化する分子間相互反応を観測することも困難である。
2つ目の従来方法は、図36(2)に示すように、細胞1内の標的分子に蛍光タンパク(プローブ)を発現させ、細胞内の分子間相互反応を観測するものである(非特許文献2,3を参照)。この場合、分子に結合させる蛍光タンパクの種類が最大で3種程度であるため、1~2種類の分子間相互反応しか観測することができないといった問題がある。
【0004】
一方、神経細胞間の情報伝達は、軸索と樹状突起で構築されるシナプスを介して行われるが、神経系細胞の発達段階において、神経ネットワークを確立するために軸索誘導は不可欠である。軸索誘導は、神経細胞から軸策が誘導因子に従って伸長して標的細胞に向かう。神経突起の先端にある成長円錐は、細胞外の様々なシグナルに対して応答することにより、神経線維を正確な標的細胞まで導くことが知られている。成長円錐の誘導方向の決定の主な要因として誘導因子の存在がある。成長円錐は、細胞外のシグナルを細胞膜上の受容体により受取り、軸索が標的細胞へと誘導される。また、誘導因子のセカンドメッセンジャーの濃度によって誘導方向がスイッチングする現象が知られ、細胞内の分子群による分子ネットワーク(以下、分子システム)と誘導因子によるシグナルとの関連性など成長円錐の誘導メカニズムの解明が要望されている。西山らは、成長円錐の膜電位を人工的に変化させることにより、成長円錐の誘導方向を決定できることを示し(非特許文献4を参照)、成長円錐の膜電位が膜電位依存のイオンチャネルを制御して成長円錐の誘導方向を決定するといったモデルを提唱している。
【0005】
しかしながら、成長円錐の膜電位を制御する細胞内分子システムについては部分的にしか解明されていないのが実情である。神経突起先端の成長円錐のように、構造的に小さく、応力に弱い細胞では、実験が困難であり、限られた観測データに基づいて、細胞内分子システムの分子間相互作用のメカニズムを解明していかなければならない。
【0006】
本発明者の一人である山田は、西山らが示した成長円錐の生物学的知見を用いて、成長円錐の膜電位の時系列データを用いて、分子システムの同定を試みている(非特許文献5を参照)。非特許文献5では、誘導因子(Semaphorin 3A;Sema 3A)のセカンドメッセンジャー(細胞外のシグナルを細胞内へ伝達する際に変換される分子)であるcGMPの細胞内濃度を人工的に操作した条件下のシグナルによる膜電位変化の時系列データを用い、西山らが提唱する分子システムにおいて、未だ明らかでない分子間相互作用の解明を行っている。具体的には、西山らが計測した膜電位変化の時系列データに対し、生化学反応方程式を用いて数理モデル化する。そして、同一条件においても、細胞により膜電位変化の時系列データがばらつくことから、このデータのばらつきを数理モデルのモデルパラメータのばらつきとして扱い、ベイズ推定により、分子システムの同定を行っている。
【0007】
山田が行った分子システムの同定では、モデルパラメータとして、cGMPの拡散時定数、カリウムチャネルの膜電位寄与度、Hill係数、イオンチャネル濃度電位変換係数、定常速度係数、解離定数などを用い、これらのモデルパラメータを一括りとして全ての細胞に共通するモデルパラメータとして扱っていた。また、数理モデルの候補については、1~3種類のモデルパラメータを用いた限定されたモデル候補から尤もらしい分子システムを示す数理モデルを選択していた。このため、モデル候補に挙がっていない数理モデルや、細胞毎に存在する個性を考慮できず、細胞内分子システムの分子間相互作用のメカニズムを解明するには不十分であった。
【先行技術文献】
【0008】

【非特許文献1】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).
【非特許文献2】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).
【非特許文献3】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).
【非特許文献4】Makoto Nishiyama, et al., Membrane potential shifts caused by diffusible guidance signals direct growth-cone turning, Nature Neuroscience volume 11, pages 762-771 (2008).
【非特許文献5】山田達也,奈良先端科学技術大学院大学 情報科学研究科 情報生命科学専攻 修正論文,NAIST-IS-MT1051204(2012年2月).
【発明の概要】
【発明が解決しようとする課題】
【0009】
上述の如く、非特許文献4に開示された従来技術では、モデル候補に挙がっていない数理モデルや、細胞毎に存在する個性を考慮できず、細胞内分子システムの分子間相互作用のメカニズムを解明するには不十分であるという問題がある。
また、細胞内分子システムの生化学反応については部分的にしか解明されておらず、構造的に小さく、応力に弱い細胞では、限られた観測データに基づいて、分子間相互作用のメカニズムを解明していく必要がある。
【0010】
そこで、本発明は、膜電位変化の時系列データといった間接情報を用いて、細胞が生きている状態で発生している分子間相互反応の推定、すなわち、直接観測が困難な細胞内に存在する分子群の生化学反応経路の推定方法及び装置を提供することを目的とする。
【課題を解決するための手段】
【0011】
上記課題を達成すべく、本発明の生化学反応経路の推定方法は、観測対象の細胞内に存在する分子群の生化学反応経路の推定方法であって、膜電位変化を生成する生化学反応経路の数理モデルに対して、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとの少なくとも2種類のパラメータを設定する。そして、誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データと同様な変化を示す数理モデルをベイズ推定により抽出し、生化学反応経路を推定する。
細胞内の生化学反応経路の結果として表れる膜電位変化の計測は、従来から神経科学の技術として確立されている。本発明では、計測可能な膜電位を用い、膜電位変化の時系列データと同様な変化を示す数理モデルをベイズ推定により抽出して生化学反応経路を推定する。その際、数理モデルに対して、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとの少なくとも2種類のパラメータを設定することで、推定精度を向上させる。
生化学反応の伝搬をトリガーする化学的刺激を細胞に与えると、反応経路の出力として膜電位変化が起きる。この膜電位変化を生成する経路の候補を全て数理モデルで定量的に表現し、観測される膜電位変化と同様の変化を示す数理モデルを抽出する。個々の細胞が持つ個性により、観測される膜電位変化は様々な系列を持つため、抽出した数理モデルにおいて、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータの少なくとも2種類のパラメータを確率分布として表現し、ベイズ推定により最適な反応経路を推定する。
【0012】
本発明の生化学反応経路の推定方法では、具体的には、下記1)~4)のステップを備える。
1)誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測する計測ステップ。
2)膜電位変化を生成する生化学反応経路の全候補を数理モデルで表すモデル化ステップ。
3)上記2)の数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定するパラメータ設定ステップ。
4)上記1)の膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、上記3)の事前確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出し、事前確率分布で指定するパラメータの範囲で膜電位変化と同様な変化を示す数理モデルを抽出し、生化学反応経路を推定する推定ステップ。
【0013】
ここで、細胞固有パラメータは、具体的には、セカンドメッセンジャーの細胞内基質内での拡散時定数、及び、各イオンチャネルによる膜電位の変化量である。
また、細胞共通パラメータは、具体的には、Hill係数、各分子の正規化反応速度、各イオンチャネル密度の反応速度、及び、各イオンチャネルの膜電位への寄与度の平均値である。
【0014】
本発明の生化学反応経路の推定方法において、生化学反応経路の全候補の最大数は、分子群の内で経路を構成する分子間相互作用の経路数と、各分子からイオンチャネルへ向かう径路数との合計数をべき数とし、分子間相互作用の態様数を基数として算出される。また、分子間相互作用の態様は、具体的には、相互作用無し、活性作用有り、抑制作用有り、の3態様である。
例えば、細胞内分子システムに2つの分子が存在して相互作用がある場合は、分子群の内で経路を構成する分子間相互作用の経路数が2となる。また、細胞内に2つのイオンチャネルがあるとし、2つの分子からそれぞれ特定のイオンチャネルへ向かう場合には、各分子からイオンチャネルへ向かう径路数が2となる。また、分子間相互作用の態様数が、相互作用無し、活性作用有り、抑制作用有りの3態様の場合は、基数が3となる。その場合、生化学反応経路の全候補の最大数は、3=81となる。詳細については後述の実施例で説明する。
【0015】
本発明の生化学反応経路の推定方法において、パラメータ設定ステップにおける事前確率分布は、膜電位変化へのイオンチャネルの作用を抑制するか否かの制御パラメータと、分子群の何れかの分子の作用を抑制するか否かの制御パラメータとの少なくとも何れかが、更に含まれることでもよい。細胞外の環境に応じて、イオンチャネルの作用、分子の作用を抑制した場合の生化学反応経路を推定できる。
【0016】
本発明の生化学反応経路推定プログラムは、上記の生化学反応経路の推定方法の各ステップを、コンピュータに実行させるためのプログラムである。
【0017】
本発明の生化学反応経路の推定装置は、観測対象の細胞内に存在する分子群の生化学反応経路の推定装置であって、下記1)の膜電位変化時系列データ計測部と、下記2)の生化学反応経路推定部を備える。
1)誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測する膜電位変化時系列データ計測部。
2)膜電位変化を生成する生化学反応経路の全候補を数理モデルで表し、数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、パラメータの事前確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出し、事前確率分布で指定するパラメータの範囲で膜電位変化と同様な変化を示す数理モデルを抽出し、分子群の生化学反応経路を推定する生化学反応経路推定部。
【発明の効果】
【0018】
本発明によれば、観測困難な細胞内の生化学反応を直接観測することなく、計測容易な非直接的観測データである膜電位変化の時系列データといった間接情報を用いて、細胞が生きている状態で、細胞内に存在する分子群の生化学反応経路の推定できるといった効果がある。
【図面の簡単な説明】
【0019】
【図1】本発明の生化学反応経路の推定方法の概念図
【図2】本発明の生化学反応経路の推定方法の説明図(1)
【図3】本発明の生化学反応経路の推定方法の説明図(2)
【図4】細胞内の分子システムの情報伝達の説明図(1)
【図5】細胞内の分子システムの情報伝達の説明図(2)
【図6】細胞内の分子システムの情報伝達の説明図(3)
【図7】細胞内の分子システムの情報伝達の説明図(4)
【図8】数理モデル候補の説明図
【図9】本発明の生化学反応経路の推定方法のフロー図
【図10】他の細胞内の分子システムの情報伝達の説明図(1)
【図11】他の細胞内の分子システムの情報伝達の説明図(2)
【図12】実施例1の細胞内の分子システムの情報伝達の説明図(1)
【図13】実施例1の細胞内の分子システムの情報伝達の説明図(2)
【図14】実施例1の膜電位の計測についての説明図
【図15】膜電位時系列データの説明図
【図16】実施例1の細胞内の分子システムの数理モデル候補の説明図
【図17】cGMPの拡散シミュレーションの説明図
【図18】膜電位時系列データと数理モデルの尤度についての説明図
【図19】膜電位時系列データセットの説明図
【図20】数理モデル候補の尤度の推定結果の説明図(1)
【図21】数理モデル候補の尤度の推定結果の説明図(2)
【図22】実施例1の交差検証手順の説明図
【図23】平均適合度の評価の説明図
【図24】実施例1の細胞内の分子システムの情報伝達の説明図(3)
【図25】細胞共通パラメータの説明図
【図26】実施例2の交差検証手順の説明図
【図27】実施例2の膜電位の計測についての説明図
【図28】実施例2の数理モデル候補の尤度の推定結果の説明図(1)
【図29】実施例2の数理モデル候補の尤度の推定結果の説明図(2)
【図30】実施例3の予測可能性検証手順の説明図
【図31】実施例3の数理モデル候補の尤度の推定結果の説明図
【図32】実施例4の予測可能性検証手順の説明図
【図33】実施例4の数理モデル候補の尤度の推定結果の説明図
【図34】実施例5の定常膜電位変化の説明図(1)
【図35】実施例5の定常膜電位変化の説明図(2)
【図36】従来の生化学反応経路の推定方法の説明図

【発明を実施するための最良の形態】
【0020】
本発明の生化学反応経路の推定方法について、図1,2を参照して説明する。図1の概念図に示すように、本発明の生化学反応経路の推定方法では、観測対象の細胞1に対して、細胞膜2上にある受容体3と細胞外の誘導因子(図示せず)が反応し結合すると、二次的に情報伝達物質であるセカンドメッセンジャーが産生される。そして、セカンドメッセンジャーのシグナルにより、細胞内に存在する分子群に細胞周辺の情報を伝達していき、最終的に、細胞膜にある膜貫通タンパク質の一種であるイオンチャネル(7a,7b)を介して、情報を隣接細胞に伝える。イオンチャネル(7a,7b)は、細胞膜表面に多種類のイオンチャネル分子を備え、これらのイオンチャネル分子が細胞内外のイオンを出し入れする。
【0021】
本発明の生化学反応経路の推定方法は、上記のように観測対象の細胞内に存在する分子群の情報伝達の生化学反応経路の推定を行うものであり、観測データとしてイオンチャネル(7a,7b)の膜電位時系列データを計測し、存在する反応経路を推定する。
図2(1)に示すように、実験系において、生きている細胞の膜電位変化の時系列データを観測する。例えば、塩素(Cl)イオンチャネル(DCl)7aとナトリウム(Na)イオンチャネル(DNa)7bの膜電位時系列データを計測する。また、細胞内に存在する分子群の情報伝達経路に関して、候補の反応経路の数理モデルに対して、パラメータを設定し、観測した膜電位時系列データと誤差が少ない反応経路を抽出する。細胞により膜電位変化の時系列データがばらつくことから、このデータのばらつきを数理モデルのモデルパラメータのばらつきとして扱う。すなわち、本発明の生化学反応経路の推定方法は、図3に示すように、候補に挙げた複数の数理モデル(1,2,・・・,n)から、膜電位変化の時系列データと同様な変化を示す数理モデルiを、ベイズ推定により抽出し、生化学反応経路を推定する。
【0022】
ここで、細胞内に存在する分子群の情報伝達の生化学反応経路について、図4~7の細胞内の模式図の一例を参照して説明する。
図4に示すように、観測対象の細胞1に対して、誘導因子を注入(8a)し、細胞膜2上にある受容体3と誘導因子を結合させると、細胞内の分子群(図示せず)が活性化され、細胞膜上のイオンチャネル(7a,7b)の活性密度が調整され、その結果として、細胞外からのイオンの流入量が変化し膜電位変化を生じる。その膜電位変化の観測10して、細胞が生きている状態で、細胞内に存在する分子群の生化学反応経路を推定する。
【0023】
図5に示すように、誘導因子を注入(8a)すると、二次的に情報伝達物質であるセカンドメッセンジャー4が産生される(8b)。そして、セカンドメッセンジャー4のシグナルにより、例えば、細胞内に存在する第1/第2の分子(5,6)が活性化され(8c,8d)、イオンチャネル(7a,7b)の活性密度、すなわち活性状態の各イオンチャネル分子の濃度が調整される。
【0024】
図6に示すように、例えば、第1/第2の分子(5,6)からイオンチャネル(7a,7b)までに至る分子間の反応経路は未知な部分が多いとする。過去の実験などから、第1の分子5の活性が起因してイオンチャネル7aの活性密度が調整される反応経路(8e,9a)が存在すること、及び、第2の分子6の活性が起因してイオンチャネル7bの活性密度が調整される反応経路(8f,9b)が存在すること、これらが既知としても、それ以外の分子間相互作用については未知であるとする。このような場合に、本発明の推定方法は、未知の分子間相互作用の反応経路を推定できる。すなわち、図7に示すように、本発明の推定方法は、例えば、第1/第2の分子(5,6)からイオンチャネル(7a,7b)までの未知の分子間相互作用の反応経路(11a,11b,12,13)の有無を推定できる。
【0025】
本発明の推定方法では、第1/第2の分子(5,6)からイオンチャネル(7a,7b)までの未知の分子間相互作用の反応経路(11a,11b,12,13)について、全ての候補を数理モデルとして表現する。分子間相互作用として、相互作用無し(No interaction)、活性作用(Activation)、活性抑制作用(inhibition)の3通りがあるとする。第1の分子5から第2の分子6への反応経路11a、第2の分子6から第1の分子5への反応経路11b、第1の分子5からイオンチャネル7bへの反応経路12、及び、第2の分子6からイオンチャネル7aへの反応経路13の4通りがあり、それぞれの反応経路において、相互作用無し(No interaction)、活性作用(Activation)、活性抑制作用(inhibition)の3通りがある。そのため、反応経路(11a,11b,12,13)の全ての候補は、3=81通り存在することになる。
【0026】
図8に示すように、81通りの反応経路(11a,11b,12,13)をマトリックス状に整理する。図8において、右上に示す通り、マトリックス内に示す記号“-”,“→”,“⇒”は、それぞれ相互作用無し(No interaction)、活性作用(Activation)、活性抑制作用(inhibition)を意味する(矢印の向きは反応経路の向きに合せて示す)。
図8のマトリックスの上辺は、第1の分子5から第2の分子6への反応経路11aを示し、左辺は、第2の分子6から第1の分子5への反応経路11bを示し、右辺は、第1の分子5からイオンチャネル7bへの反応経路12を示し、下辺は、第2の分子6からイオンチャネル7aへの反応経路13を示している。例えば、Mは、反応経路11a,11b,12,13の全てが相互作用無しであることを意味し、Mは、反応経路11a,13の2つが活性抑制作用として働き,反応経路11b,12は相互作用無しであることを意味している。同様に、M73は、反応経路11b,12の2つが活性抑制作用として働き,反応経路11a,13は相互作用無しであることを意味し、M81は、反応経路11a,11b,12,13の全てが活性抑制作用として働くことを意味している。例えば、Mの下のM18(図示せず)では、反応経路11bが活性作用として働き、反応経路11a,13の2つが活性抑制作用として働き、反応経路12は相互作用無しであることを意味している。
【0027】
本発明の推定方法では、上述の未知の反応経路(11a,11b,12,13)の全ての候補を数理モデルとして表現し、誘導因子のセカンドメッセンジャーのシグナルによる膜電位変化の時系列データを計測結果と突き合わせ、81通りの候補の数理モデルに対する尤もらしさから、膜電位変化と同様な変化を示す数理モデルを抽出し、生化学反応経路を推定する。
本発明の推定方法では、数理モデルにおける、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定する。本発明の推定方法では、細胞内の分子群の相互反応は、個々の細胞に依存せず全ての細胞で共通しているプロセスと、個々の細胞の個性に依存しているプロセスが存在することを前提とする。細胞で共通しているプロセスを数理モデルに表現するものが、個々の細胞に依存しない細胞共通パラメータであり、細胞の個性に依存しているプロセスを数理モデルに表現するものが、個々の細胞に依存する細胞固有パラメータである。
図7に細胞1内には、細胞で共通しているプロセスと、細胞の個性に依存しているプロセスとが共存している概念を符号1a,1bで示している。
【0028】
本発明の推定方法では、数理モデルにおいて、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、その事前確率分布と、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出して、81通りの候補の数理モデルから、膜電位変化と同様な変化を示す数理モデルを抽出する。生化学反応経路を推定する。
【0029】
本発明の推定方法のフローを図9に示す。観測対象の細胞に誘導因子を注入し(S01)、注入した誘導因子のセカンドメッセンジャーの濃度変化によって誘導方向が変わり、細胞内の分子群の相互作用によって、膜電位変化が生じる(S02)。公知の計測方法を用いて、膜電位時系列データを計測する(計測ステップ;S03)。そして、膜電位変化を生成する生化学反応経路の全候補を数理モデルで表現する(モデル化ステップ;S04)。表現した数理モデルにおいて、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定する(パラメータ設定ステップ;S05)。なお、モデル化ステップ(S04)、パラメータ設定ステップ(S05)は、計測ステップ(S03)の以前、又は、併行して行っても構わない。
そして、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布と、個細胞共通パラメータと細胞固有パラメータとに関する事前確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出し、計測した膜電位変化の時系列データを用いて、膜電位変化と同様な変化を示す数理モデルを抽出し、反応経路を推定する(推定ステップ;S06)。
【0030】
上述の説明では、誘導因子のセカンドメッセンジャーのシグナルによって生じる反応経路、すなわち、第1/第2の分子(5,6)からイオンチャネル(7a,7b)までの未知の分子間相互作用の反応経路(11a,11b,12,13)について、数理モデルが81通りの存在することを示した。
本発明の推定方法は、様々なケースに応用できる。例えば、図10に示すように、誘導因子のセカンドメッセンジャー4のシグナルにより、3つの分子(15~17)が活性化され(8c,8d,8e)、3つのイオンチャネル(7a,7b,7c)の活性状態の各イオンチャネル分子の濃度が調整され、膜電位変化が生じるとする。上述と同様に、分子間相互作用として、相互作用無し(No interaction)、活性作用(Activation)、活性抑制作用(inhibition)の3通りがあるとする。
【0031】
図11に示すように、3つの分子(15~17)の相互反応経路が6つ、3つの分子(15~17)それぞれから3つのイオンチャネル(7a,7b,7c)への反応経路が9つ存在する。そのため、反応経路の全ての候補は、315通り存在することになる。但し、既に何れかの反応経路が確認できており、既知の反応経路(図7の反応経路9a,9bのように既知の経路)が存在するのであれば、候補数を減らせることができる。また、イオンチャネルは、透過するイオンの選択性により、ナトリウム(Na)イオンチャネル、塩素(Cl)イオンチャネル、カルシウム(Ca)イオンチャネル、カリウム(K)イオンチャネルが存在するが、誘導因子のセカンドメッセンジャーのシグナルによって膜電位変化が生じるイオンチャネルに絞りこんで反応経路の候補を決定できる。
【0032】
以下、本発明の実施形態の一例を、図面を参照しながら詳細に説明していく。なお、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更及び変形が可能である。
【実施例1】
【0033】
上述したように、西山らは、成長円錐の膜電位を人工的に変化させることにより、成長円錐の誘導方向を決定できることを示し、成長円錐の膜電位が膜電位依存のイオンチャネルを制御して成長円錐の誘導方向を決定するといった細胞内分子システムのモデルを提唱している。
本実施例では、西山らが提唱する分子システムにおいて、未知の反応経路の推定を行った結果を示す。
【0034】
先ず、誘導因子Semaphorin 3A(以下、Sema 3A)のセカンドメッセンジャーであるcGMPの細胞内濃度を人工的に操作した。実際には、cGMPと同等の効果を持つアナログである8-Br-cGMPを用いた。図12に示すように、誘導因子のSema 3Aは、成長円錐の細胞膜2上の受容体3に受容されると、cGMPへとシグナル変換される。cGMPはリン酸基を有するエネルギー運搬体であり、PKG(protein kinase G)やCNGC(cyclic nucleotide gated ion channel)の分子と結合する(8c,8d)。PKGやCNGCの分子が直接的または間接的に塩素イオンチャネル(ClC)7a及びナトリウムイオンチャネル(NaC)7bへ作用し、膜電位変化を生じさせる。ここで、過去の西山らの実験によって、PKGはナトリウムイオンチャネル(NaC)7bを介して膜電位の脱分極(静止膜電位から膜電位を上げる働き)を引き起こし、また、CNGCは塩素イオンチャネル(ClC)7aを介して膜電位の過分極(静止膜電位から膜電位を下げる働き)を引き起こすことが既知である。しかし、それ以外の分子間相互作用については判明しておらず、図13に示すような反応経路(9a,9b)の存在は既知であるが、反応経路(11a,11b,12,13)が存在するか否かは未知である。
本実施例では、図13に示す反応経路(11a,11b,12,13)の存在について、膜電位時系列データを用いて推定する。
【0035】
図14に示すように、cGMPのアナログである8-Br-cGMPを、マイクロピペット30を用いて観測対象の成長円錐細胞31に注入し、その後の膜電位変化を計測した。実験は生物にとって通常に近い条件下で実施した。8-Br-cGMPの濃度を変えて、幾つかの濃度で実験を行った。注入後10分程度の膜電位時系列データを計測した。なお、成長円錐の膜電位の計測は、マイクロピペット30を固定し、Whole Cell Recordingによって計測した。図15に、8-Br-cGMPの濃度を変えて膜電位変化を計測した様子を示す。図15のグラフにおいて、横軸は時間(秒)、縦軸は成長円錐膜電位(mV)を示す。8-Br-cGMPの濃度を変えると、膜電位時系列データは異なる。
【0036】
図16は、PKGとCNGCの2つの分子からClCとNaCの2つのイオンチャネルまでの未知の分子間相互作用の反応経路(11a,11b,12,13)について、全ての候補をマトリックスに表現したものである。分子間相互作用として、相互作用無し(No interaction)、活性作用(Activation)、活性抑制作用(inhibition)の3通りがある。CNGCからPKGへの反応経路、PKGからCNGCへの反応経路、CNGCからNaCへの反応経路、及び、PKGからClCへの反応経路の4通りがあり、それぞれの反応経路において、相互作用無し、活性作用、活性抑制作用の3通りがあり、反応経路の全ての候補は、3=81通り存在する。
【0037】
図16のマトリックスの上辺は、CNGCからPKGへの反応経路11aを示し、左辺は、PKGからCNGCへの反応経路11bを示し、右辺は、CNGCからNaCへの反応経路12を示し、下辺は、PKGからClCへの反応経路13を示している。Mは、反応経路11a,11b,12,13の全てが相互作用無しであることを意味し、Mは、反応経路11a,13の2つが活性抑制作用として働き,反応経路11b,12は相互作用無しであることを意味し、M73は、反応経路11b,12の2つが活性抑制作用として働き、反応経路11a,13は相互作用無しであることを意味し、M81は、反応経路11a,11b,12,13の全てが活性抑制作用として働くことを意味している。
【0038】
未知の反応経路(11a,11b,12,13)の全ての候補を数理モデルとして表現し、膜電位変化の時系列データを計測結果と突き合わせ、81通りの数理モデルに対する尤もらしさから、膜電位変化と同様な変化を示す数理モデルを抽出する方法について、以下に詳細に説明する。
まず、数理モデルにおいて、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定する。具体的な細胞共通パラメータと細胞固有パラメータを、それぞれ下記表1,2に示す。
【0039】
【表1】
JP2019154299A_000003t.gif

【0040】
【表2】
JP2019154299A_000004t.gif

【0041】
細胞内の分子群の相互反応において、個々の細胞に依存せず全ての細胞で共通しているプロセスと、個々の細胞の個性に依存しているプロセスが存在するとし、細胞で共通しているプロセスを細胞共通パラメータ(上記表1の10個のパラメータ)で表現し、細胞の個性に依存しているプロセスを細胞固有パラメータ(上記表2の4個のパラメータ)で表現する。
【0042】
そして、ベイズ推定を用いて、81通りの候補の数理モデルから、膜電位変化と同様な変化を示す数理モデルを抽出する。具体的には、個々の細胞に依存しない細胞共通パラメータと、個々の細胞に依存する細胞固有パラメータとに関する事前確率分布を設定し、その事前確率分布と、膜電位変化の時系列データの数理モデルに対する尤もらしさを表す確率分布とから、細胞共通パラメータ及び細胞固有パラメータに関する事後確率分布を算出する。以下、具体的な数式を示しながら説明する。
【0043】
実験では、アフリカツメガエルの脊髄ニューロンの成長円錐細胞を用い、成長円錐細胞の中心付近にcGMPのアナログ(8-Br-cGMP)の溶液で満たされた記録ピペットを取り付け、8-Br-cGMPを成長円錐全体に拡散させた。実験における8-Br-cGMPの刺激濃度は数μMであり、成長円錐の内在的なcGMP(数10nM)との間に大きな濃度差があるため、ピペットから成長円錐細胞内への8-Br-cGMPの流入には拡散過程を考慮した。成長円錐Sにおける8-Br-cGMP濃度の時系列は、時間tの単純指数関数として近似され、下記式(1)によって拡散過程として表すことができる。
ここで、Smaxはピペット中の8-Br-cGMP濃度であり、これが入力濃度刺激の最大値になる。τsは拡散時定数である。 拡散時定数τsは、成長円錐細胞の体積に依存し、上記表2に示す細胞固有パラメータである。
【0044】
【数1】
JP2019154299A_000005t.gif

【0045】
拡散時定数τsの変化を表現するために、図17に示す3つの第1~第3区画(1,2,3)からなるモデル成長円錐における数値計算から推定された確率分布としてモデル化した。内側の2つの第1,2区画(1,2)は1μmの厚さを有し、外側の第3区画(3)はガウスノイズε~N(0,1)によって修正された1~2μmの可変厚さを有する。第n区画の8-Br-cGMP濃度S(μM)はフィックの法則に従うべきであり、下記式(2)として表すことができる。
【0046】
【数2】
JP2019154299A_000006t.gif

【0047】
ここで、D、An-1,n、dn-1,n、Vは、8-Br-cGMP(1μm/s)の拡散係数、n-1番目とn番目の区画中点間の距離、および第n区画の容積にそれぞれ対応するものである。上記式(2)の境界条件S=Smaxを用いて数値的に計算し、第2区画(3)の時間経過が上記数式1の指数関数で近似した。種々の成長円錐サイズにおけるモンテカルロシミュレーションによって、拡散時定数τsは、平均が40秒である非負のガウス分布とした。
【0048】
上述の如く、反応経路の数理モデルでは、8-Br-cGMPの結合時のCNGCやPKGの活性化の反応経路の存在を考慮した。以下、X,Yとして表されるCNGCおよびPKGの下流因子(DF)の活性レベルが、8-Br-cGMPの濃度レベルが増加するにつれて増加し、DFの活性レベルは、正常な状態(または飽和)に達するまでの時間で観測した。上記数式SのHill型関数として、下記式(3),(4)に示すように、XとYを表現する準定常近似を適用した。ここで、Hill係数n,mは、下流プロセスに対して効果的なトリガーを行う。 X及びYは、[0,1]の範囲で正規化され、実際の活動レベルはその最大値を掛けて表現される。
【0049】
【数3】
JP2019154299A_000007t.gif

【0050】
上記式(3),(4)における解離定数KおよびKは、細胞内の8-Br-cGMP濃度とDF活性との間の反応の順方向速度(KxbおよびKyb)に逆方向速度(KxfおよびKyf)の比によって定義される。 順方向速度および逆方向速度に直線近似を適用すると、有効解離定数は、下記式(5),(6)で表すことができる。ここで、分母と分子はそれぞれ順反応速度と逆反応速度を表す。 数式中の8つのパラメータの値は、細胞共通パラメータである。
【0051】
【数4】
JP2019154299A_000008t.gif

【0052】
膜電位(MP)変化は、イオンチャネル分子の密度に依存する。このため、数理モデルは、CNGCおよびPKGによるClCおよびNaCのイオンチャネル分子の密度の調節を考慮している。イオンチャネル分子の密度のダイナミクスは、それぞれDCl、DNaとして表されるイオンチャネル密度の比の変化によって定義される。DCl、DNaの調整はXやYに依存するため、細胞共通パラメータは、数理モデルにおける信号処理の速度(微小時間スケール)を制限することになる。 したがって、生化学反応のダイナミクスは、常微分方程式(ODE:Ordinary Differential Equation)によって記述され、下記数式で表すことができる。ここで、添え字付きのkは、正規化された変数X,Yの定数または関数のいずれかであり、有効な順方向反応速度または逆方向反応速度である。
【0053】
【数5】
JP2019154299A_000009t.gif

【0054】
線形近似によって、生化学反応のダイナミクスは、下記式(9)~(12)として表すことができる。下記式において、反応速度k**0は一定であり(XまたはYとは無関係)、他の反応速度はXまたはYによって与えられる最大速度である。下記式(9)~(12)の全てのパラメータは、細胞共通パラメータであり、上記式(5)および(6)のパラメータを用いて、数理モデルを表す。
【0055】
【数6】
JP2019154299A_000010t.gif

【0056】
上記式(9)~(12)において、イオンチャネルからCNGCおよびPKGへのフィードバックの反応経路を仮定しなかった。いくつかの状態を伴うフィードバック反応経路は、発散、振動、双安定性などのシステム状態の不安定さをもたらすが、そのような不安定な応答は、実験では膜電位時系列データでは観察されなかった。
【0057】
イオンチャネル密度(DClおよびDNa)が膜電位(MP)を調整し、イオンチャネル密度からMPへの信号変換については、HodgkinおよびHuxleyの式によって十分に確立されていることから、膜電位は下記式(13)の線形表現として表すことができる。
【0058】
【数7】
JP2019154299A_000011t.gif

【0059】
ここで、Cはイオンチャネル容量と膜電位に依存する膜容量、Vは膜電位、VNa、V、VClは反転電位、gNa、g、gClはチャネルコンダクタンスである。静止状態でのカリウムコンダクタンスgは、支配的な透過性(g>>gCl、gNa)を有する。上記式(13)をゼロに設定することによって、定常状態および静止状態での膜電位(MP)は、下記式(14)で表すことができる。
【0060】
【数8】
JP2019154299A_000012t.gif

【0061】
Na、g、gCl(上付きバー)は、定常状態におけるチャネルコンダクタンスである。数理モデルでは、DClおよびDNaは8-Br-cGMPの下流因子によってイオンチャネルの分子密度が調節される。 静止電位における膜電位Vの小さな変化のため、8-Br-cGMPのシグナルによる膜電位変化は、以下の直線近似式に従うことになる。
【0062】
【数9】
JP2019154299A_000013t.gif

【0063】
ClおよびDNaは、それぞれ8-Br-cGMPのシグナルによって誘導されたClCおよびNaCの有効濃度の正規化された変化であり、それらの変化が最大に達した場合には1つになる。また、パラメータV、AClおよびANaは個々の細胞特性に依存し、その値は各細胞に特異的であり、上記表2に示す細胞固有パラメータである。
【0064】
成長円錐体積およびイオンチャネル密度のパラメータが膜電位変化に影響を及ぼすことから、ベイズ推定を用いる。数理モデルのパラメータは、細胞共通パラメータと細胞固有パラメータ と実験環境として強制的に導入された実験条件パラメータの3つに分類される。
生化学反応速度、Hill係数およびAClおよびANaの平均(ACl0およびANa0)を含む細胞共通パラメータセット(θ)は、すべての細胞に共通であるため、全ての膜電位時系列データセットから推定した。一方、個々の細胞の特性(例えば、成長円錐体積、チャネル密度)に大きく依存する細胞固有パラメータセット(φ)= {V、ACl、ANa、τ}を各膜電位時系列データについて推定した。実験条件パラメータセット(c)= {ηCl、ηNa、Smax}は、薬理学的条件および適用された8-Br-cGMP濃度などの実験的に誘導されたパラメータである。
【0065】
ベイズ推定のアプローチにより、事前分布(上記表1に示すように、殆どガウス分布の非負部分)によって与えられた制約下で、計測した膜電位時系列データからモデルパラメータを推定した。モデルパラメータの事後分布を下記式(16)で表した。
【0066】
【数10】
JP2019154299A_000014t.gif

【0067】
ここで、実験条件パラメータセットcは、見やすさのために省略されている。エビデンスE(V,M)は、膜電位時系列データセットVおよび数理モデルMによって決定されるモデル評価基準である。事前分布のパラメータp(φ、θ|M)は、上述の表1,2に示したものである。与えられたデータセットに対する適応レベルを表す尤度は、ガウス分布の積によって与えられる。
【0068】
【数11】
JP2019154299A_000015t.gif

【0069】
ΔV(t)は、観測された膜電位とV(t)との間の誤差を示す。ここで、iは膜電位時系列データのインデックスであり、I=16、T=300~900である。i番目の膜電位時系列データのσiはデータ前処理中に計算された観測ノイズの大きさである。
なお、膜電位の観測データは、細胞の自発的興奮等の影響により、観測ノイズ以外の予期しないトレンドが生じるため、これらのノイズをデータの前処理によりデータ補正したものを推定に用いる時系列データとしている。具体的には、明らかにはずれ値であるとわかる区間については、範囲を指定し、線形補間によりデータを内挿した。
【0070】
次に、推定結果について説明する。
数理モデル間で、計測した膜電位時系列データの尤度を比較し、上述の81通りの候補の数理モデルのそれぞれのモデル適合度を調べ、モデル適合度を評価してモデル抽出を行った。事前確率分布によって制約された細胞共通パラメータ(θ)と細胞固有パラメータ(φ)を用いて数理モデルの膜電位時系列データを計算した。上述の表1,2に示したように、パラメータの事前確率分布のほとんどにゼロ平均のガウス分布の非負部分を使用した。
【0071】
そして、単一の実験的膜電位時系列データと数理モデルの膜電位時系列データとの間の誤差を尤度関数として比較した。実験的膜電位時系列データに細胞依存性ノイズが存在するため、各ノイズレベルで誤差を標準化した。サンプリングされた膜電位時系列データと同じ膜電位時系列データの平滑化された時間平均との間の誤差からノイズレベル(σi)を推定した。
このプロセスは、すべての膜電位時系列データの全データ点について反復され、すべてのデータセットのモデル適合度として全尤度を得るために全ての尤度が乗算された。パラメータ値に依存するモデル適合度は、実験的膜電位時系列データと数理モデルの膜電位時系列データとの間の類似性のみを提供し、パラメータの妥当性を無視する。したがって、数理モデルがデータセットにあてはまるときにパラメータが不合理な値を取る可能性があるため、モデルの適合度だけでは、もっともらしい数理モデルを選択するには不十分である。
【0072】
上述の如く、ベイズ推定を使用して、データセットのオーバーフィットからパラメータを回避し、パラメータの妥当性を組み込むことによって、数理モデルの妥当性を計算した。数理モデルの妥当性は、モデル適合度としての尤度とパラメータの妥当性としての事前確率との積として表される。
ここで、尤度として表されるモデル適合度は、図18(4)に示す尤度(likelihood)グラフに示すように定義される。定数パラメータを無視して上記式(18)の自然対数をとることによって得られる。エビデンスの対数を数理モデルの妥当性の指標とし、エビデンスの実際の計算は、マルコフ連鎖モンテカルロ法(MCMC)による事後分布からのサンプリングにより実行した。
【0073】
実験的な膜電位時系列データセットの数は16であり、その内訳は、10μM 8-Br-cGMPによる時系列データ数n=7、膜電位変化へのClCの寄与をブロッキングするDNDSを実験条件として加えた時系列データ数n=5、膜電位変化へのNaCの寄与をブロッキングするSTXを実験条件として加えた時系列データ数n=4である。図19に示すように、全て(16セット)の実験的な膜電位時系列データセットから得られる細胞共通パラメータをθallとして表す。また、細胞に依存する細胞固有パラメータ(φi)は、モデル適応度を計算するための全ての膜電位時系列データセットからi番目の膜電位時系列データからのみ推定されたものを表す。10μM 8-Br-cGMPから誘発される膜電位時系列データセットから細胞内の分子システム内のシグナル伝達ネットワークを同定することを目指して、検証、特異性、および予測可能性のための全データセットから推定された細胞共通パラメータ(θall)を利用したテストを行った。
【0074】
図20に示すように、各々の数理モデルの実験的膜電位時系列データセットから計算されたエビデンスは、より小さなエビデンス(グレーの薄い領域)を有するグループaと、より大きなエビデンスを有するグループb(グレーの濃い領域)を有していた。
【0075】
図21に示すように、各グループ内の数理モデル間のエビデンスの違いは、a,bの2つのグループ間のものよりもずっと小さく、グループbの数理モデルに一般的に存在する反応経路が重要な役割を果たす可能性があることを示唆している。 注目すべきことは、グループb内に存在する最も単純な分子相互作用を表すMは、膜電位変化中にPKGを介してClCを抑制するcGMPのシグナル伝達に関与することを明らかにしていることである。
【0076】
細胞共通パラメータがどの程度安定しているかを評価するために、交差検証(CV;Cross-validation)を実行することによって、推定された細胞共通パラメータの特異性を確認した。図22に示すように、最初に16個のデータセットを2つのサブセットに分割する。1つのサブセットは、15のデータセットであり、細胞共通パラメータ(θiMAP、i番目の膜電位時系列データを除く)を訓練する。もう1つのサブセットは、膜電位時系列データセットを1つ(抜き出したデータセット)だけ有し、細胞固有パラメータ(φi; i番目の膜電位時系列データセット)を訓練する。抜き出した膜電位時系列データセットは、訓練されていないデータとしての可能性を計算するために使用される。これは、細胞共通パラメータに組み込まれていないためである。全データセットサンプル(i=1,2,・・・,16)について16回訓練と計算手順を繰り返し、モデル適合度を得るために16の尤度を乗算した。
【0077】
細胞共通パラメータの推定値が安定している場合には、Mの細胞共通パラメータは交差検証によって大きく変更されるべきではなく、訓練されていないデータセットに対するモデル適応度のマトリックス表現は訓練データセットのそれと同様でなければならない。
図23(1)と(2)を比較すると、データポイントあたりの対数尤度で表されるモデル適合度の平均値(θiMAPとの平均適合度)は、グループaでは小さく、グループbでは大きく見える。
【0078】
交差検証において、各々の訓練ごとに細胞共通パラメータ(θiMAP; i=1,2,・・・,16)を推定した。図25に示すように、Mの細胞共通パラメータ(θallMAP)がθiMAPの±1の標準偏差(s.d.)内に入っている。これらの結果は、本推定方法が、交差検証のデータセットの組み合わせによって著しく影響を受けることなく、細胞共通パラメータが適度に安定していることを示している。
【0079】
以上述べたように、誘導因子Sema 3AのセカンドメッセンジャーであるcGMPのアナログである8-Br-cGMPを用いて、成長円錐の細胞に刺激を与えることにより、塩素イオンチャネルとナトリウムイオンチャネルの動態により変化する膜電位時系列データを計測し、計測した時系列データ用いて、本発明の推定方法を適用し、未知の生化学反応経路を推定できた。
【実施例2】
【0080】
図26は、交差検証の手順を示している。図26に示す検証手順は、上述の図22に示す検証手順と異なり、エビデンス評価に全く用いていないデータを用いて細胞固有パラメータφを推定した。最初に16個のデータセットを用いて、細胞共通パラメータ(θallMAP)を訓練し、その後、エビデンス評価に全く用いていないn個のデータセットを用いて、細胞固有パラメータ(φi; i番目の膜電位時系列データセット)を訓練し、候補の数理モデルを評価した。
【0081】
本実施例では、下記表3の試薬を用いて、特定の分子の働きを抑制した条件下で膜電位時系列データを計測した。図27に、膜電位時系列データ計測の実験手順を示す。図27に示すように、時点Aにおいてチャネルブロッカー(STX/DNDS)やPKG抑制薬(KT5823)を導入し、十分な時間(30分以上)経過後の時点Bにおいて、cGMPのアナログである8-Br-cGMPを導入すると同時に膜電位の観測を行った。
【0082】
【表3】
JP2019154299A_000016t.gif

【0083】
図28(1)は、10μM 8-Br-cGMPのデータセットを用い、チャネルブロッカー(STX/DNDS)を導入した条件下で、チャネルブロッカーの効果(ηCl、ηNa)の値(存在/非存在)を入れ替えたデータセットを用いた検証結果、図28(2)は、5μM 8-Br-cGMPのデータセットを用い、刺激を弱め、(1)と同様に、チャネルブロッカー(STX/DNDS)の導入条件下で、チャネルブロッカーの効果(ηCl、ηNa)の値を入れ替えたデータセットを用いた検証結果を示している。図28に示すように、共にcGMPによって駆動される細胞内分子システムが機能していることが確認できた。
【0084】
一方、図29(1)は、10μM 8-Br-cGMPのデータセットを用い、PKG抑制薬の導入条件下でのデータセットを用いた検証結果、図29(2)は、Netrin(脊髄交連神経細胞の軸索を誘引する分泌タンパク質として同定された軸索ガイダンス分子)の環境下でのデータセットを用いた検証結果を示している。共に、cGMPによって駆動される本来の分子システムとは異なるシステムが機能していることが確認できた。
【実施例3】
【0085】
図30は、予測可能性検証の手順を示している。10μM 8-Br-cGMPの16個のデータセットを2つのサブセットに分割し、1つのサブセット(15のデータセット)を用いて、細胞共通パラメータ(θiMAP、i番目の膜電位時系列データを除く)を訓練する。もう1つのサブセット(i番目の膜電位時系列データ)において、時系列の過渡期(0~250秒)の膜電位時系列データを用いて細胞固有パラメータφを推定し、推定に用いていない残りの膜電位時系列データ(250秒~)に対する予測誤差を尤度により評価した。すなわち、図31(1)に示すグラフにおいて、グレイ領域が時系列の過渡期(0~250秒)の領域であり、この領域内の膜電位時系列データを用いてパラメータφを推定し、残りの領域の膜電位時系列データを尤度評価に用いた。図31(1)のグラフにおいて、黒線は実験データを表し、グレイの実線および破線はそれぞれ、M,Mによる予測した膜電位時系列データを表す。
図31(2)は、各候補の数理モデルに対し評価した尤度を示している。図31(2)から、cGMPによって駆動される細胞内分子システムが機能していると言える。
【実施例4】
【0086】
図32は、予測可能性検証のその他の手順を示している。10μM 8-Br-cGMPの16個のデータセットを用いて、細胞共通パラメータ(θallMAP)を訓練する。そして、エビデンス評価に全く用いていない11個のデータセットであって、上述の実施例と同様に、時系列の過渡期(0~250秒)の膜電位時系列データを用いて細胞固有パラメータφを推定する。エビデンス評価に全く用いていない11個のデータセットは、5μM 8-Br-cGMPにより刺激を弱めたデータである。推定に用いていない残りの膜電位時系列データ(250秒~)に対する予測誤差を尤度により評価した。すなわち、図33(1)に示すグラフにおいて、グレイ領域が時系列の過渡期(0~250秒)の領域であり、この領域内の膜電位時系列データを用いてパラメータφを推定し、残りの領域の膜電位時系列データを尤度評価に用いた。図33(1)のグラフにおいて、黒線は実験データを表し、グレイの実線および破線はそれぞれ、M,Mによる予測した膜電位時系列データを表す。
図33(2)は、各候補の数理モデルに対し評価した尤度を示している。図33(2)に示すように、cGMPによって駆動される細胞内分子システムが機能していることが確認できた。
【実施例5】
【0087】
膜電位は、定常状態では、ほぼ一定の値を維持する。これは、8-Br-cGMP刺激のレベルに双方向の依存性を有するからである。実施例1における図16のマトリックスのカラムMにおいて強調されたClCのPKG抑制を含む細胞内分子システムが、cGMP依存性の双方向の膜電位変化の誘発に必要か否かについて確認した。異なる3通りの濃度のcGMPと、10μMの8-Br-cGMP(n=16)から推定されるθallMAPを用いてMを計算した。
モデルシミュレーションは、(1)低濃度(0.5μM)の8-Br-cGMPによって刺激されて生じる過分極、(2)適度な8-Br-cGMPによって刺激され徐々に回復するシャープな過分極、(3)高濃度(20μM)の8-Br-cGMPによって刺激され、脱分極に変換されたシャープな過分極の3通り実施した。
【0088】
次に、モデルの定常状態の膜電位(図34の長時間経過後の円点を参照)と実験データを比較した。図34は、数理モデルMによる各cGMP濃度に対する膜電位時系列データの一例を示している。図35は、数理モデルMによるcGMP濃度依存性の膜電位変化と実験データの比較を示している。図中、●印は、図34と対応している。シミュレートされたcGMP濃度依存性の膜電位変化は、双方向の膜電位変化を示す(図35の点線のグラフを参照)。
モデルシミュレーションは双方向の膜電位変化を示したが、脱分極の大きさははるかに大きく、脱分極の発生は実験的に観測されたものよりはるかに低い8-Br-cGMP濃度で現れた。この差は、実験での8-Br-cGMPの浴適用(bath)と比較して、記録ピペットを通した8-Br-cGMPによる直接刺激との差から生じるものである。
8-Br-cGMPの膜浸透の単純なHill様モデルを導入したとき、8-Br-cGMPは、実験的に観測されたものにはるかに近い双方向の膜電位変化を再現することができた(図35の黒四角を参照)。
【産業上の利用可能性】
【0089】
本発明は、神経細胞などの細胞内分子群の生化学反応経路の推定に有用である。
【符号の説明】
【0090】
1 細胞
2 細胞膜
3 受容体
4 セカンドメッセンジャー
5 第1の分子
6 第2の分子
7a~7c イオンチャネル
8a~8f,9a,9b,11~13 反応経路
10 膜電位変化の観測
30 マイクロピペット
31 成長円錐細胞
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16
【図18】
17
【図19】
18
【図20】
19
【図21】
20
【図22】
21
【図23】
22
【図24】
23
【図25】
24
【図26】
25
【図27】
26
【図28】
27
【図29】
28
【図30】
29
【図31】
30
【図32】
31
【図33】
32
【図34】
33
【図35】
34
【図36】
35