代謝経路のモデル化について -

目次

• 代謝モデル化ソフトウェア　Scrumpy について
• An auxin-driven polarized transport model for phyllotaxis.　　　Jonsson H, et al.　　　Proc Natl Acad Sci U S A. 2006 Jan 31;103(5):1633-8 の、Fig.3 に出てくるグラフのようなものを表現する小さなモデル(セルが5個しかない)
• ローレンツ方程式に基づく？モデル

代謝モデル化ソフトウェア　Scrumpy について

植物の代謝経路全体のモデル化に成功した論文が発表された。

A_Genome-scale_metabolic_model_of_Arabidopsis_and_some_of_its_properties. 　　　Poolman MG, Miguet L, Sweetlove LJ, Fell DA.　　　Plant Physiol. 2009 Nov;151(3):1570-81. Epub 2009 Sep 15.

Poolman 博士は この他にも代謝のモデル化に関して論文を多数発表している。自らScrumpy という Matabolic modeling package を開発している。

今では、生物の細胞をモデル化し、数値計算することが簡単にでき、Scilab のように使えるソフトウェアが公開されている。　その一つである、CellDesignerを使えるようにする手順

今となっては Scrumpy はインストールしにくく（新しいPCに入れたら一部の機能が使えなくなり、グラフが出なかった。うまく直せていない）機能もたいしたことはないが、一応記録として置いておく。

ウインドウが出てきたら

 MyModel = ScrumPy.Model()

と打つ。コピーしておいてウインドウから CTRL-Y でペーストできる。 ファイルダイアログが出てくるので、適当なモデルファイル (.spy) を指定する。 例として添付されている、Simple.spy がわかりやすい。 "Scrumpy" メニューが出現する。まず Compile をする。エラーがなければなにもおきない。次に Simulate すると、四角のウインドウが一つ、steps, step size を設定するウインドウが一つ出てくる。 さらに solve をする。「Get steady state」とメッセージボックスが出る。　出ないときは隠れているので、探し出して　OK を押さないと、反応がなくなってしまう。 つぎに Inspector をすると、パネルが出てくる つぎに dynamic monitor をする。グラフ画面が出てくる。右に代謝産物の名前がリストされる。ダブルクリックしてセレクトする。 セレクトしたら、steps, step size を設定するウインドウの 「OK」を押す。するとグラフが描かれる。同時にコンソール画面に「integrate to yy : ended at yy」と表示される。 Simple.spy では、「Get steady state」と返事が返ってくる。そのことと一致して、ステップを進めていってもそれぞれの反応の速度(流速 flux)、化合物の量は一定に近い値をとり変化しない。 エラーがででうまくグラフが描かれないときは、モデルの読み込みからやり直す。

メニューと同じことを、コマンドで行うことができる。

  定常状態を求める
  MyModel = Scrumpy.Model() で、現在読み込まれているモデルを指定した。ので、
  
  MyModel.FindSS()　と打つ。メニューから実行するのと同じことが起きる。
  実行例
  >>> MyModel = ScrumPy.Model()
  >>> MyModel.FindSS()
  >>> MyModel["mGH"]
  2.7237045658365977
  >>> MyModel["x_RS"]=1
  >>> MyModel.FindSS()
  >>> MyModel["mGH"]
  3.5650850653913384

.py ファイルを作り、計算手順を指定することができる。

 # MyModel = ScrumPy.Model()　あらかじめモデルを読み込んでおく　その後で、このファイルを IDLE から実行する

 i = 0.5
 while i <= 1:
    MyModel["x_RS"] = i
    MyModel.FindSS()
    print "RS = ", i,"  GS = ", MyModel["mGS"]
    print "    Total = ", MyModel["mGS"]+MyModel["mGS2"]
    i = i + 0.1

このソフトウェアは代謝のモデル化のために開発されたものだが、Flux (流量)を用いてモデル化できる物事なら、どんなことでも分析させ平衡状態を求めることができる。または時間的な状態変化（時間発展）を調べることができる。「バスタブモデル」というモデル化の方法に相当するらしい。

現象から方程式を創り出す--立式で鍛える論理的思考法(4)砂漠化問題のバスタブモデル [in Japanese]　　　o 西浦 廉政 o 寺本 敬　　　 数理科学 48(5), 70-77,　2010-05 　　　サイエンス社

バスタブには流入口と流出口がある。それぞれの流量 (flux) がパラメーターの変化によって調節される。複数のバスタブが組み合わさることで複雑なモデルが構成できる。バスタブを平面的、立体的に並べることで二次元、三次元の構造を表現できる。拡散を表現したければ、「となりのバスタブへ、濃度差に応じて流れる」というようになる。

この「バスタブモデル」というのは、要するに微分方程式によるモデル化に相当するが、生物学的な物事を表現するのに考えやすく、非常に適している。西浦、寺本両博士も感染症のモデル、砂漠化問題のモデルについて解説している。さらに優れたソフトウェアが開発されることで、微分方程式の解き方がわからなくても計算機が解いてくれるようになった。

バスタブで生物学の現象をモデル化した論文の例：

A thermodynamic model of transcriptome formation.　　　Konishi T.　　　Nucleic Acids Res. 2005 Nov 24;33(20):6587-92. Print 2005.　　　

小西博士による論文の Figure 1 では、バスタブ（mRNA の量を表現する）を用いた図式化がされている。 さらに流入速度（合成）、流出速度（分解）を数式で表現し、それによって予想される物事について議論されている。

細胞内の合成や分解は質量作用の法則に基づいて表現できる。それらには反応定数が必ず付随する。アレニウスの式によって、温度と活性化エネルギーから反応定数が決定できる。

因子（この場合はプロモーターや mRNA）が多種類ある場合、活性化エネルギーはそれぞれの因子ごとに異なる。統計的に見ると、それらは正規分布する（サイコロを12回振って、出た目を足し合わせて得られる値をたくさん集めると正規分布に近くなることが知られている。それと同じパターン）。

アレニウスの式の指数関数の肩に活性化エネルギーが乗っている。そのために、反応定数の方は対数正規分布になる。結果的に転写産物全体は対数正規分布になりやすい（というようなことが書いてあるような気がする）。 転写というのは非常に複雑な、いくつもの化学反応が関わっている過程だが小西博士の成果に基づいてシンプルな形で表現できる。

「複雑系にひそむ規則性 : 対数正規分布を軸にして(交流)」　　http://ci.nii.ac.jp/naid/110008723070　　　という解説が読める。この解説によると、「ある変数の成長（時間変化または変化の系列）を生み出す微視的なプロセスが乗算過程で近似されるときに、対数正規分布が出現する」というようなことが書いてある。

あるステップ（X 番目）での変数の値　＝　[その変数の初期値] * A1 * A2 * A3 * ... * AX

A1 〜 AX は、それぞれのステップでの、変化を表す定数

両辺の対数を取る。ln(ABCDE) = ln(A) + ln(B) + ln(C) + ln(D) + ln(E) という具合だから、

ln(あるステップ（X 番目）での変数の値) = ln([その変数の初期値]) + ln(A1) + .... + ln(AX)

右側の値が全部無相関で、個数が十分多ければ、左側の値が正規分布になる。

細胞内には非常に多くの種類のmRNAが存在する。それをモデル化しようとすれば、バスタブを複数並べればよいことになる。それぞれのバスタブについて、パラメーターが決められる。モデルを複雑に、実際に近くしようとすれば、「バスタブ同士がつながっている」 「翻訳産物のバスタブを追加する」ようにする。

西浦、寺本両博士の解説では、74ページの図2のCに、モデルの模式図が書いてある。こういう図は生物学の論文によく見られる。今まではこういう図を書いてそれでわかったような気になって終わりになっていた。 しかし、こういう図を75ページにあるようにモデル式に書き直し、さらにソフトウェアの助けを借りて分析ができるかもしれない。 西浦、寺本両博士は、非常に素人にもわかりやすく解説を書かれている。すばらしいことである。素人でもうまくモデルを作り、ソフトウェアの力で分析することによって、そのモデルがあり得ない挙動を示していないかどうか、チェックできるかもしれない。いままで見落とされていたことがわかるようになるかもしれない。 新しい重要な問題が見いだされるかもしれない。なにか有用なことが予想できるかもしれない。他人の論文にケチをつけるのに有効かもしれない。しかしソフトウェアに頼り切るのでは大間違いをする事もあり得る、限界が出てくるので数学の勉強もしなければならない。

微分方程式を使ったモデルの場合、「解の安定性」が重要な因子になるそうである。それを調べるには「ヤコビ行列の固有値を求める方法」が使われるそうである。何か行列があれば、それを分解してみたり固有値、固有ベクトルを求めてみたりすることは定石である。 様々な分野で固有値に何らかの意味がでてくることが知られている。「分岐図」「定常状態の数」も重要な因子だそうである。西浦、寺本両博士の解説でも安定性や分岐図について解説されている。

代謝（複数の化学反応をまとめた系）のモデル化では、化学量論行列　stoichiometric 行列が重要な因子になる。化学反応系の揺らぎについて研究されている。 http://www.sat.t.u-tokyo.ac.jp/~ryotat/nc.pdf　　　遺伝子発現の揺らぎのモデル予測: 構造的アプローチ　　　冨岡亮太　博士らの研究

微分方程式を用いたモデルは様々な分野で用いられる。経済学を例にして解説した資料　http://www.fides.dti.ne.jp/q-gmz/Math/Math_5.pdf　　経済学のモデルでは均衡点が安定かどうかが重要になる。ヤコビ行列が出てくる。

生物の細胞内には様々なフィードバックの仕組みが存在する。細胞内の転写、翻訳、生合成には時間がかかるので、「時間遅れのあるフィードバック」という仕組みが多くなる。そういう時間遅れが、システムの性質を決めるのに重要だという研究成果が既に発表されている。　　　http://www.science-plaza.or.jp/about/sankangaku/forum10/bessyo.pdf　　　別所博士のすばらしい成果　　　Scrumpy だと、「時間遅れ」を表現するのが難しい。これも次々と優れたソフトウェアが開発されているので問題なくなっている。

Scilab/scicos というソフトは、生物学向けに作られたわけではないが、一般的な微分方程式を用いたモデルを表現するために使われている。解説書も出ている。

今では、生物の細胞をモデル化し数値計算することが簡単にでき、Scilab のように使える優れたソフトウェアが公開されている。　その一つである、CellDesignerを使えるようにする手順

ためになる論文

Computer simulation of the dynamic behavior of the glutathione-ascorbate redox cycle in chloroplasts.　　　Valero E, González-Sánchez MI, Maciá H, García-Carmona F.　　　Plant Physiol. 2009 Apr;149(4):1958-69. Epub 2009 Feb 25. PMID: 19244456

試しに Scrunpyのモデルを作った。たいしたものはできていない。

An auxin-driven polarized transport model for phyllotaxis.　　　Jonsson H, et al.　　　Proc Natl Acad Sci U S A. 2006 Jan 31;103(5):1633-8 の、Fig.3 に出てくるグラフのようなものを表現する小さなモデル(セルが5個しかない)

 #
 # Auxin
 #
 A1_to_A2: 					# Cell 1 to 2, 2 to 1
 A1 <> A2					# Reversible
 D*(A1 - A2) + T*(A1*P1_2 - A2*P2_1) # Reaction rate equation = dA2/dt

 A1_to_A5:
 A1 <> A5
 D*(A1 - A5) + T*(A1*P1_5 - A5*P5_1)

 A2_to_A3:
 A2 <> A3
 D*(A2 - A3) + T*(A2*P2_3 - A3*P3_2)

 A3_to_A4:
 A3 <> A4
 D*(A3 - A4) + T*(A3*P3_4 - A4*P4_3)
 
 A4_to_A5:
 A4 <> A5
 D*(A4 - A5) + T*(A4*P4_5 - A5*P5_4)
 
 #
 # PIN1
 #
 PIN1_to_2:
 P1_0 <> P1_2 		# P1_0 means internal PIN1 in cell 1
 (k1 * A2 * P1_0 - k2 * P1_2)

 PIN1_to_5:
 P1_0 <> P1_5 		# P1_0 means internal PIN1 in cell 1
 (k1 * A5 * P1_0 - k2 * P1_5)

 PIN2_to_3:
 P2_0 <> P2_3 		# P1_0 means internal PIN1 in cell 1
 (k1 * A3 * P2_0 - k2 * P2_3)

 PIN2_to_1:
 P2_0 <> P2_1 		# P1_0 means internal PIN1 in cell 1
 (k1 * A1 * P2_0 - k2 * P2_1)

 PIN3_to_4:
 P3_0 <> P3_4 		# P1_0 means internal PIN1 in cell 1
 (k1 * A4 * P3_0 - k2 * P3_4)

 PIN3_to_2:
 P3_0 <> P3_2 		# P1_0 means internal PIN1 in cell 1
 (k1 * A2 * P3_0 - k2 * P3_2)

 PIN4_to_5:
 P4_0 <> P4_5 		# P1_0 means internal PIN1 in cell 1
 (k1 * A5 * P4_0 - k2 * P4_5)

 PIN4_to_3:
 P4_0 <> P4_3 		# P1_0 means internal PIN1 in cell 1
 (k1 * A3 * P4_0 - k2 * P4_3)

 PIN5_to_1:
 P5_0 <> P5_1 		# P1_0 means internal PIN1 in cell 1
 (k1 * A1 * P5_0 - k2 * P5_1)

 PIN5_to_4:
 P5_0 <> P5_4 		# P1_0 means internal PIN1 in cell 1
 (k1 * A4 * P5_0 - k2 * P5_4)

# # Initialize #

 A1 = 1
 A2 = 0.98
 A3 = 1
 A4 = 1
 A5 = 1

 P1_2 = 0.5
 P1_5 = 0.5
 P2_3 = 0.5
 P2_1 = 0.5
 P3_4 = 0.5
 P3_2 = 0.5
 P4_5 = 0.5
 P4_3 = 0.5
 P5_1 = 0.5
 P5_4 = 0.5

 P1_0 = 10
 P2_0 = 10
 P3_0 = 10
 P4_0 = 10
 P5_0 = 10

 k1 = 1
 k2 = 0.01

 D = 1
 T = 1

ローレンツ方程式に基づく？モデル

 「step 0.01, 100 step 進める」で、dynamic monitor をすると、確かに値が複雑な変化をする。ローレンツ方程式もヤコビ行列で分析できるらしい。

 #
 # Lorenz equation ???
 #
 re1:
    x <> x_1
    D * (x - y)
 #
 re2:
    y <> y_1
    y + x * z - R * x
 #
 re3:
    z <> z_1
    B * z - x * y
 #
 #
 # Initialize
 #
 x = 1
 x_1 = 1
 y = 1
 y_1 = 1
 z = 1
 z_1 = 1

 D = 10
 B = 8/3
 R = 28

vim: set ts=8 sts=2 sw=2 et ft=a111_modified_flexwiki textwidth=0 lsp=12: