Quick start guide for molecular circadian clock model


Gonze-Goldbeter model[1]

例として,アカパンカビが創りだす mRNA およびタンパク質発現に関する概日振動子モデルを考えます.このモデルは,Goldbeter 等によって提案され,Gonze によって詳細な解析が行われました.このモデルに周期的外力(ここでは光による明暗サイクルを想定しています)が印加された系を取り上げましょう.この時,ある条件下では,この概日振動子は周期的な外力へ同調(引き込み)します.しかし,周期外力の周期,その振幅を変化させると,周期解の Neimark-Sacker 分岐が発生します.この発生パラメータ集合を (omega,vsmax)-パラメータ平面上で求めてみましょう.この概日振動子モデルのモデル方程式 は

GG_eq1.png

として記述されます.また vs(t) は

GG_eq2.png

によって記述される周期的な正弦波関数を仮定します.ここで GG_eq3.png と置き換えておきます. 方程式の詳細は原著論文[1]または[2],[3]を参照してください.

次の3つのプログラムを作成し,順番に用いることによって計算を行います.

  1. PP
  2. FIX
  3. BF

環境設定

  • まずインストールを参考に BUNKI をインストールし,環境設定を行って下さい.
  • 補足.1
    • Symbolic Math Toolbox (MATLAB 追加パッケージ)を利用出来ない方は,デモ用として以下のサンプルを用意してあります.このモデルの分岐解析に必要な解析ツール,それら設定ファイルなどが 入った Project ディレクトリをひとまとめにしたファイルです.適当な場所に展開して下さい.
    • 注意:BUNKI コマンドで展開したディレクトリを読込むと SE が起動します.ここで system create を行うと,Symbolic Math Toolbox が無い場合はエラーを返されます.Symbolic Math Toolbox が無い場合は,プログラムの作成 の項は飛ばして,PP に進んでください.
  • 補足.2
    • プログラム BUNKI を展開したディレクトリ内にサンプル集を集めた sample ディレクトリーがあります.
      /BUNKI/Sample/Gonze-Goldbeter
      このディレクトリ下に置かれた種々のファイルは,以下の説明で使用される設定ファイルです.ファイルを読み込むことで,各フィールドに書き込む手間が省けます.用意された設定ファイルは以下の通り:
      • fileGG.sed : モデル方程式が記述されたファイル:SE 用の設定ファイル
      • filePd_point.ppd : 安定周期解の例:PP 用の設定ファイル
      • filePd_point.fixd : Neimark-Sacker 分岐を検出:FIX 用の設定ファイル
      • filePd_point_NS.bfd : 周期解の Neimark-Sacker 分岐点用設定ファイル:BF 用の設定ファイル

プログラムの作成

これからの作業はサブ・ディレクトリで行うことにしましょう. ディレクトリの名前はなんでもよいです.まず MATLAB command window 内で以下を実行します.

>> BUNKI

実行後,Select BUNKI Project Directory という window が開きます.

project.png

ここで解析を行う Project を選択します.この段階では,まだ何も Project を生成していないので, Name 欄に識別用の Project 名を命名して OK ボタンを押します.ここでは, この方程式の解析を中心的に進めた Gonze, Goldbeter にちなんで, Project 名を GG として話を進めます*1

project2.png

新規に Project を生成する場合,ディレクトリ(or フォルダ)を生成しても良いかを尋ねられます. Yes を選択して下さい. SE(system editor) が起動します.

startSE.png

さて,いよいよ解析プログラムの作成です*2

  • システムの記述
    • まず解析する方程式がどういった種類の系であるかを選択します.GG model は,3変数, 非自律系方程式 ですので,Dimension は 3, type は Nonautonomous をチェックします.

      startSE2.png

    • Equations 欄に モデル方程式を記述します.その際,以下のルールに従って記述して下さい.
      • 状態変数:各状態変数は x[i], i = 1,2,...,n として下さい.すなわち,M->x[1]fc.png->x[2]fn.png->x[3]と読みかえて式を記述して下さい.
      • パラメータ: パラメータであることを識別するために, $ 記号でパラメータを挟む形で記述して下さい.例えば,$para$,$Vm$,$alpha$,etc,.$の記号で挟まれたものはパラメータとして認識されます.そのパラメータ名の命名は自由です.
      • 正弦波外力: ここで利用できる周期的外力は, 単一周波数をもつ正弦波に限ります. すなわち cos(t) または, sin(t) です.各方程式は,角周波数 omega で正規化されます*3

        startSE3.png

  • 文法チェック&プログラム生成
    • 以下が入力終了後の様子です.方程式の記述に矛盾が無いかをチェックするために, check ボタンを押して下さい.文法チェックが行われます.

      startSE4.png

    • 文法チェックを行った後,エラーが無ければ,いよいよプログラムを生成します.どのプログラムを生成するかを右下の check ボタンで制御します.分岐曲線を求める為には,PPFIXBF 全てにチェックを入れて下さい.Create ボタンを押して下さい.プログラム生成を開始します.

      startSE5.png

      (a): システム生成開始時に status bar が表示される.

      startSE6.png

      (b): システム生成完了後,log window に正常生成のメッセージが出力される.



PP

ここでは,パラメータと状態の初期値を読み込んで位相面図を描き,安定な周期解をつかまえます.

まず SE から PP(phase portrait) プログラムを起動しましょう. Program -> PP を選択します.

pp.png

描画 window とコントロールパネルが起動します.初回起動時には,描画 window は, (x[1],x[2])平面となっています*4.ここでは,描画 window の変更は必要ないでしょう.

pp2.png

次にコントロールパネルです.まず初期状態点を入力します*5.対応する状態(x[1],x[2],x[3])にそれぞれスタートさせたい初期点を記述します.ここでは原点(0,0,0)からスタートさせてみましょう!

描画の前にパラメータ値を記述することを忘れないで下さい.初回起動時には全てのパラメータに1 が入力されます.それでは,各パラメータを以下のように設定します.

parametervalue
KI1
Kd0.13
Km0.5
k10.5
k20.6
ks0.5
vd1.4
vm0.505
vsmax3
vsmin1.6
omega0.261


それでは,Startボタンを押してシミュレーションを開始しましょう.原点からスタートした解軌道は,ある周期解に収束していく様子が描かれたでしょうか?青の実線解軌道を,赤点は,時刻 2pi*周期毎にサンプリングされた点,すなわちPoincare写像点を表しています.

pp3.png

パラメータはリアルタイムに変化することができます.取り敢えずパラメータ omega を変えてみましょう. omega は周期的外力の周期に関係します.現在この外力の周期は約 24 [h](T = 2π/omega) に設定されています.ちょっとデフォルトのパラメータ変化が0.1と大きいので,この点をまず調整しましょう.

メニュー: Setting -> Parameter を選択して下さい.

pp4.png

すると,パラメータのキー設定*6画面が開きます.ここで最後に表示される omega の増減サイズを 0.1 から 0.001 へ変更しましょう.変更が完了したらOKです.

pp5.png

さてパラメータ omega を増加することで何が起こるでしょうか? omega を大きくしていくと安定だった周期解はやがて不安定化します.例えば, 次の例は, omega = 0.317 になったときの周期解の変化の様子を示しています.

pp7.png

もう少し見やすくしましょう. メインパネル Plot 内にある Orbit のチェックを外します. すると表示されていた解軌道が表示されなくなり,位相平面上には点列だけが示される様になります.十分時間が経過すると,これらの点列はある閉曲線を描きます.これが不変閉曲線 (Invariant Colosed Curve:ICC) と呼ばれるものです.この様な安定性変化の背景には,周期解の Neimark-Sacker 分岐*7の発生が関与しています.この場合は,omega の増加によって,安定な周期解が不安定化し,安定な ICC が観測できる Neimark-Sacker 分岐(Supercritical NS 分岐)が発生したことになります.

pp8.png

さて,元の設定に戻します. omega を減少させます.パラメータの減少ボタンを押しても良いですし,パラメータ値を直接元の数値に戻してやることも可能です.やはり安定な周期解に解軌道が収束する様子が観察できるでしょう.

pp9.png

次に,Plot 欄にある states flow にチェックを入れてみます.こうすることで,数値積分を行っている現在の状態変数値リストをリストボックス内でチェックできます.リストボックス内には,2π周期毎にサンプリングされた状態変数値が出力されます.この例では,GG 系は Poincare 写像の安定固定点へ収束していきますので,リストボックス内には,同じ数値が流れていくのを確認できます.

pp10.png

また,この時,Clear ボタンを押すことで,描画 window 内に表示されていた過渡応答分の軌道が Clear されて定常状態に達した時の様子を見ることができます.Orbit のチェックを外して,解軌道の描画も停止させましょう.十分時間が経過した後,安定固定点へ収束します.そのため,描画 window 上では安定固定点に相当する一点だけが表示されます.

pp11.png

それでは,この周期解(Poincare 写像の固定点)の分岐を調べてみましょう.パラメータを変化させた時の周期解に対する固有値(特性乗数)を評価します.この目的には FIX ツールが利用できます.FIX を動かす為の初期値を PP から与えます.ここで重要な事は,定常状態に至った時点での状態変数値,パラメータ値を FIX に与えなければならないということです.

定常状態へ至ったかどうかは,上述した states flow スイッチを入れて,状態変数値に変化がなくなったことを確認するか,Staus ボタンを数度クリックしてみて,リストボックスに表示される情報に変化が無いことを確認する必要があります.

pp12.png

変化がなくなった事を確認した後,メニュ−バーの Tools から, Export current status を選択します.この操作により,定常状態に至った時点での状態変数値,パラメータ値,その他の設定値が保持され,次のプログラムに引き継ぐことが可能になります(FIX の項目を参照).

pp13.png

(a): export 選択.

pp14.png

(b): 設定値の export 完了.

FIX

ここでは,パラメータを変化させて周期解を追跡し,そのときの固有値を表示させることが目的です.

FIX 動作の設定値はユーザ側でそれぞれ調整する必要があります.ここでは簡単のために,そのまま動作させてみましょう.FIX プログラムは,PP プログラムで保持された値 (Export Current status) をデータを読み込んでニュートン法への初期値として実行するのが通常の方法です*8

まず PP コントロールパネルからFIX プログラムを起動しましょう. Program -> FIX を選択します.

fix.png

初回起動時には,各フィールドにはデフォルト値が入力された状態で起動します.

fix2.png

PPで求められた周期解の固定点に関する情報を初期設定値としてインポートします.メインパネル Tools -> Import initial point を選択します.

fix3.png

選択後,PP を利用して得たパラメータ値,固定点の座標値などが各項目に反映されます.

fix4.png

PP でパラメータ omega を大きくすることで, 系のダイナミクスに変化が起きることが判ったので,変化させるパラメータを omega とし、その増分ステップサイズを設定します.

fix5.png

Start ボタンを押して計算を開始します.計算開始と同時に,別ウィンドウが開きます.これは,Gauss 平面上のどの位置に固有値が位置しているかを示す図です.固定点の Neimark-Sacker 分岐の発生条件は,複素共役な固有値が"単位円を横切るとき"です.パラメータの変化とともに,ウィンドウ上の1対の固有値が"単位円"へと近づいていく様子を観察できます.

パラメータを連続的に変化させていき,分岐点が検出されるとプログラムは停止します*9.上の例では,複素共役な固有値が絶対値1となったため Neimark-Sacker 分岐の発生を検出した様子を示しています.

fix6.png

停止した時のパラメータ値,状態変数値は分岐曲線を追跡するプログラム BF の初期点として利用できます.停止した際のパラメータ値などを BF 用に Export します.ここで Export BF point を選択すると,状態・パラメータ値,その他必要な設定値が保持されます.この保持された値を次のプログラムに引き継ぎます(BF の項目を参照).

fix7.png

(a): Export 選択.

fix8.png

(b): リストボックス内に Export した情報が表示される.

BF

ここでは,パラメータを変化させて分岐曲線を求めます.

BF 動作の設定値はユーザ側でそれぞれ調整する必要があります.ここでは簡単のために,そのまま動作させてみましょう.BF プログラムは,FIX プログラムで検出された分岐点での状態変数やパラメータ値を Newton 法の初期値として読み込んで実行するのが通常の方法です*10

まず FIX メインパネルから BF プログラムを起動しましょう. Program -> BF を選択します.

bf.png

初回起動時には,各フィールドにはデフォルト値が入力された状態で起動します.

bf2.png

FIX で求めた分岐点に関する各種の情報を初期値として Import します.メインパネル Tools -> Import BF point を選択します.

bf3.png

選択後,分岐の種類, 固有値の偏角情報などの付属情報が各項目に入力されます.(omega,vsmax) 平面上の分岐曲線を追跡することを考え, X パラメータとして omega を, Y パラメータとして vsmax にチェックを入れます.また step size,end などを設定します.計算結果をデータファイルとして残す場合は, Output に出力ファイル名を書いておきます.

bf4.png

ここまで設定が終われば,後は Start ボタンを押して分岐曲線の追跡を開始します.そのとき別ウィンドウとして,Eigenvalues window と Bifurcation diagram window が起動します.Eigenvalues window は,分岐曲線を計算している際,Neimark-Sacker 分岐に対応する固有値だけではなく,それ以外の固有値の配置変化を視覚化する目的で表示されます.例えば,余次元2の分岐点が発生したことなどを視覚的に理解するために使用します.Bifurcation diagram window は,Neimark-Sacker 分岐曲線を追跡していく様子をリアルタイムに確認するための window です.

bf5.png

(a):固有値配置図とリストボックスの様子.

bf6.png

(b):(omega,vsmax) 平面の分岐図.リアルタイムに計算結果を見ることができる.


参考文献


  1. J.C. Leloup ,D. Gonze, A. Goldbeter, Limit cycle models for circadian rhythms based on transcriptional regulation in Drosophila and Neurospora, J. Biol. Rhythms. 14, pp.443-448, 1999.
  2. D. Gonze and A. Goldbeter, Entrainment versus chaos in a model for a circadian oscillator driven by light-dark cycles, J. Stat. Phys. 101, pp.649-663, 2000.
  3. G. Kurosawa and A. Goldbeter, Amplitude of circadian oscillations entrained by 24-h light-dark cycles, J. Theor. Biol., 242, pp.478-488, 2006.




*1 Project 名は自由に命名出来ます.ユーザーが識別しやすい名前にすれば良いです.
*2 サンプルディレクトリの下に fileGG.sed という SE 用の設定ファイルがあります.このファイルを読み込む事で,あとはプロクラムの生成をするだけの状態にすることが可能です.すなわち,これから行う諸々の設定を省略することができます.ファイルの読み込み方法は,SE設定ファイルを読み込むには?を参照して下さい.
*3 自動的にパラメータ omega が付加されます
*4 この設定は,コントロールパネル Setting から変更することが出来ます
*5 サンプルディレクトリの下に filePd_point.ppd という PP 用の設定ファイルがあります.このファイルは以下の説明にある周期解をみるためのパラメータ設定や初期条件が保存された PP 設定ファイルです.あとは読み込んだ後にスタートボタンを押すだけでシミュレーションが可能です.ファイルの読み込み方法は,PPの設定ファイルを読込むには?を参照して下さい.なお,filechoas.ppd という PP 用の設定ファイルがあります.このファイルは,カオス振動が発生する場合の初期設定ファイルです.興味のある場合,この設定ファイルを PP で読込んでシミュレーションを楽しんで下さい.
*6 詳しくは キー入力によってパラメータを変更したい を参照して下さい.
*7 固定点の Hopf 分岐と呼ばれたりします.
*8 一方,パラメータ,状態変数値が記述さたファイルを読み込んで実行する方法もあります.サンプルディレクトリの下に filePd_point.fixd という FIX 用の設定ファイルがあります.このファイルは以下で説明される,周期解のパラメータ設定や初期条件が保存された FIX 用の設定ファイルです.読み込んだ後にスタートボタンを押すだけで数値計算を開始することが可能です.ファイルの読み込み方法は,FIX の設定ファイルを読込むには?を参照して下さい.
*9 この動作の詳細は Setting から option パネルを動作させてコントロールすることが可能です.詳細はマニュアルまたは,個別対策マニュアルを参考にして下さい.
*10 一方,パラメータ,状態変数値が予め記述さたファイルを読み込んで実行する方法もあります.サンプルディレクトリの下に filePd_point_NS.bfd という BF 用の設定ファイルがあります.このファイルは以下で説明される固定点の Neimark-Sacker 分岐曲線を追跡するためのパラメータ設定や初期条件が保存された BF 設定ファイルです.読み込んだ後にスタートボタンを押すだけで分岐曲線の追跡が可能です.ファイルの読み込み方法は,BFの設定ファイルを読込むには?を参照して下さい.

Attach file: fileeq_img05.gif 219 download [Information] fileeq_img04.gif 225 download [Information] fileeq_img03.gif 227 download [Information] fileeq_img02.gif 232 download [Information] fileeq_img01.gif 225 download [Information] fileGG-sample.zip 294 download [Information] filePd_point_NS.bfd 367 download [Information] filePd_point.fixd 398 download [Information] filePd_point.ppd 416 download [Information] filepp12.png 392 download [Information] filestartSE4.png 400 download [Information] fileGG.sed 437 download [Information] filepp9.png 408 download [Information] fileproject.png 395 download [Information] filefix8.png 378 download [Information] filefix4.png 392 download [Information] filepp5.png 382 download [Information] filebf3.png 384 download [Information] filefix.png 379 download [Information] fileGG_eq1.png 364 download [Information] filepp8.png 428 download [Information] filefn.png 359 download [Information] filepp4.png 391 download [Information] filepp13.png 378 download [Information] fileproject2.png 371 download [Information] filebf6.png 403 download [Information] filestartSE5.png 390 download [Information] filebf2.png 379 download [Information] filechaos.ppd 374 download [Information] filepp6.png 241 download [Information] filepp2.png 380 download [Information] fileGG_eq2.png 391 download [Information] filefix5.png 387 download [Information] filestartSE.png 396 download [Information] filefix3.png 382 download [Information] filestartSE3.png 395 download [Information] filepp11.png 378 download [Information] filepp14.png 396 download [Information] filestartSE6.png 381 download [Information] filefix2.png 372 download [Information] filebf4.png 379 download [Information] filefix7.png 389 download [Information] filepp10.png 380 download [Information] filestartSE2.png 379 download [Information] filebf5.png 368 download [Information] filefix6.png 380 download [Information] fileGG_eq3.png 368 download [Information] filepp.png 389 download [Information] filebf.png 385 download [Information] filefc.png 400 download [Information] filepp3.png 398 download [Information] filepp7.png 393 download [Information]

Front page   Edit Freeze Diff Backup Upload Copy Rename Reload   New List of pages Search Recent changes   Help   RSS of recent changes
Last-modified: 2009-07-23 (Thu) 20:18:53 (3376d)