準備

インストールについては別記事にありますのでそちらを参考にして下さい.

Windows 7かMac OS X上でE-Cell4を使ってみる

使い方

E-Cell4の基本要素

E-Cell4を使う上で理解しなければならない3つの要素としてModel, World, Simulatorがあり, 計算における異なる概念をあらわしている.

Modelは, その名の通り扱うモデル, すなわち問題を記述したものである.
Worldは, 状態をあらわしたものである. 例えば初期状態や計算途中のある時点における状態などである.
Simulatorは, 計算解法 (ソルバ, アルゴリズム)をあらわす.

これら3つの要素は独立しているが, WorldとSimulatorはセットで扱うことが多い (理由については後述する).

E-Cell4におけるModelの基本

まず, Modelからみていこう. Modelは扱う問題を記述したものだと書いたが, 大きく2つの要素によって構成されている.

Species, すなわち分子種は, モデルに登場する実体をあらわす.
ReactionRule, すなわち反応則は, Species間の相互作用をあらわす.

E-Cell4では主に化学反応を扱うので, それを例に言えば:

Speciesは, 基質や酵素などどんな種類の分子が存在し, それらがどのような性質 (分子の大きさなど)を持つかをあらわす.
ReactionRuleは, 結合やかい離などそれら分子間にどのような反応があるかをあらわす.

ということである.

Species

実際にE-Cell4を使ってみる. E-Cell4を使ってSpeciesを作るのは非常に簡単で:

from ecell4 import *

A = Species("A")
B = Species("B")

とすれば良い. それぞれ, AとBという名前の分子種をあらわす. 1個の特定の分子ではなく, そういう種類の分子があることを意味している.

注意点として, Speciesの名前 (シリアルSerialと呼ぶ)にはいくつかのルールがあり, 特殊記号 '()', '.', '_'や数字, 空白の使用には注意が必要である¹.

1行目はE-Cell4をPython上で簡単に使うための方法で, 几帳面な人であれば

import ecell4

A = ecell4.Species("A")
B = ecell4.Species("B")

としても良い. このページは簡単な使い方なので今後は最初の方式を取ったものとして1行目は省略する.

さて, Speciesには任意の属性を持たせることで, どのような性質を持つのかを指定できる.

A = Species("A")
A.set_attribute("radius", "0.005")
A.set_attribute("D", "1")
A.set_attribute("location", "cytoplasm")

1つ目に与える引数が属性名で, 2つ目がその値で, どちらも文字列でなければならない. 属性のうち特に分子径, 拡散係数, 属する場所は頻繁に用いられるので簡単に

A = Species("A", "0.005", "1", "cytoplasm")  # XXX: serial, radius, D, location

でまとめて与えることができる. 逆に与えられたSpeciesの属性を知りたいときは以下のようにする.

# X = Species("A", "0.005", "1")

print(X.serial())  # will return 'A'
print(X.get_attribute("D"))  # will return '1'

ReactionRule

Speciesができたところで今度は反応を作ってみる.

反応はおおざっぱに言って, 何かの分子種が何かの分子種になることで, 前者をReactants, 後者をProductsと呼ぶ. あとはその反応の速さ, 速度定数, を決めれば良い. 速度定数は数値で与える.

rr = ReactionRule()
rr.add_reactant(Species("A"))
rr.add_product(Species("B"))
rr.set_k(1.0)

これで, 分子種Aが分子種Bに変る反応を作れた. ここでは反応を定義しているので, Speciesに属性値はなくても良い. 結合もこれと同じに記述できる.

rr = ReactionRule()
rr.add_reactant(Species("A"))
rr.add_reactant(Species("B"))
rr.add_product(Species("C"))
rr.set_k(1.0)

AとBの結合によってCができる. より読み手に伝わりやすい書き方として以下の便利関数も用意されている.

rr1 = create_unimolecular_reaction_rule(Species("A"), Species("B"), 1.0)
rr2 = create_binding_reaction_rule(Species("A"), Species("B"), Species("C"), 1.0)
rr3 = create_binding_reaction_rule(Species("C"), Species("A"), Species("B"), 1.5)

逆に与えられたReactionRuleがどのようなものかを簡単に知りたい場合はas_string関数で見られる.

print(rr3.as_string())  # will return 'C+A>B|1.5'

NetworkModel

以上でModelに含まれる要素を作れるようになったので, これらをModelに登録していく.

sp1 = Species("A", "0.005", "1")
sp2 = Species("B", "0.005", "1")
sp3 = Species("C", "0.01", "0.5")
rr1 = create_binding_reaction_rule(Species("A"), Species("B"), Species("C"), 0.01)
rr2 = create_unbinding_reaction_rule(Species("C"), Species("A"), Species("B"), 0.3)

m = NetworkModel()
m.add_species_attribute(sp1)
m.add_species_attribute(sp2)
m.add_species_attribute(sp3)
m.add_reaction_rule(rr1)
m.add_reaction_rule(rr2)

Speciesの追加はadd_species_attributeで, ReactionRuleの追加はadd_reaction_ruleで行う. これで簡単な結合とかい離のモデルを作ることができた.

モデルの中身を取り出すときにはspecies_attributesとreaction_rulesを使う².

print(m.species_attributes())
print(m.reaction_rules())

# will return like:
# [<ecell4.core.Species object at 0x7f36443bfa98>, <ecell4.core.Species object at 0x7f36443bfab0>, <ecell4.core.Species object at 0x7f36443bfac8>]
# [<ecell4.core.ReactionRule object at 0x7f36443bfa98>, <ecell4.core.ReactionRule object at 0x7f36443bfab0>]

特にSpeciesの属性を必要としない計算 (常微分方程式やGillespie法)ではadd_species_attributeは省略しても良い.

E-Cell4を使ってModelを簡単に計算してみる

WorldとSimulatorの説明がまだだが, 早くシミュレーションを走らせたい人のためにここで一度計算を行なってみる.

E-Cell4には簡単なデモのための便利関数としてrun_simulationがある. IPythonなどmatplotlibが利用可能な環境で以下のようにする.

# XXX: 'm' is a NetworkModel, which is described in the section above.

run_simulation(10, {'C': 60}, volume=1.0, model=m)

run_simulationは時刻tで与えた時点で分子数を記録して時系列をプロットしてくれる. 今回は10秒間, 0.1秒刻みとした³. 2つ目の引数は初期値で, 分子種Cが60分子存在している状態から始める. volumeは計算する空間の容積で, modelに今回作ったNetworkModelを与えている.

計算がうまくいっていれば上のようなプロットが見られるはずだ. オメデトウ!

この味気ないプロットに満足できない人のためには確率論的計算手法というものがある⁴.

run_simulation(10, {'C': 60}, volume=1.0, model=m, solver='gillespie')

全く同じモデルを異る手法で計算することができた. これは問題と計算手法が分離されているからこそなせる技である.

Modelを記述するE-Cell4の特殊記法

以上, E-Cell4を使ってSpeciesとReactionRuleを作り, NetworkModelに登録して簡単に計算することができるようになった. しかし, モデルの記述はこのように簡単な結合とかい離だけであっても煩雑でモデルの読み手に伝わりづらい.

そこで, E-Cell4ではこのモデルを記述するための独自の構文をいくつか勝手に用意している. 今回はReactionRuleについてだけ見ていこう.

with reaction_rules():
    A + B > C | 0.01  # equivalent to create_binding_reaction_rule
    C > A + B | 0.3   # equivalent to create_unbinding_reaction_rule

m = get_model()

特殊記法を使う時にはwithステートメントを利用する. この下でインデントを続けている限りは特殊構文が利用できる. 'reaction_rules'のあとに'()'を忘れずに.

構文については見ての通りで, セパレータ'|'の後に反応速度定数として数値を1つだけ与える. 構文自体はあくまでもPythonであるため, 不要な改行は許されない.

また, 可逆反応については特別な記法があり, それを使えばより簡単になる.

equilibrium1.py

import numpy
from ecell4 import *

with reaction_rules():
    A + B == C | (0.01, 0.3)

run_simulation(10, {'C': 60}, volume=1.0)

\frac{\mathrm{d[A]}}{\mathrm{d}t}=\frac{\mathrm{d[B]}}{\mathrm{d}t}=-0.01\mathrm{[A][B]}+0.3\mathrm{[C]}\\
\frac{\mathrm{d[C]}}{\mathrm{d}t}=+0.01\mathrm{[A][B]}-0.3\mathrm{[C]}

可逆反応の記号は'=='で, セパレータの後にはtuple⁵を用いて2つの速度定数を与える⁶. 'run_simulation'関数は引数'model'を指定しない場合は自分で'get_model()'を呼んでモデルとするので, このように省略可能な場合がある.

発現や分解の記述

ここで特殊な例について説明を加えておく. 発現や分解の反応では, 反応の左辺か右辺がなにもないことになるが, ecell4の特殊記法では残念ながらこれらをそのまま書くことはできない. 例えば

with reaction_rules():
    A > | 1.0  # XXX: will throw SyntaxError
    > A | 1.0  # XXX: will throw SyntaxError

と書くと, 'SyntaxError: invalid syntax'というエラーが表示されるはずだ.

これらを表現するために特殊記号'~'を用いる. ReactionRuleを記述する際にSpeciesの前に'~'記号を置くことで反応上のその分子種の量論係数が0になる. 簡単に言えば, 単に無いものと見なされると考えれば良い⁷.

with reaction_rules():
    A > ~A | 1.0  # XXX: create_degradation_reaction_rule
    ~A > A | 1.0  # XXX: create_synthesis_reaction_rule

\frac{\mathrm{d[A]}}{\mathrm{d}t}=1.0-1.0\mathrm{[A]}

常微分方程式ソルバでWorldとSimulator

Modelが理解できれば, WorldとSimulatorはさほど難しくはない. さきほどの'run_simulation'関数で与えた'volume'と"{'C': 60}"などがWorldで, 'solver'がSimulatorに相当する情報だと思えばよろしい. 以下で実際に'run_simulation'がしている作業を自分で書けるようにする.

'run_simulation'は標準で'ode'と呼ばれるソルバを使う. これは常微分方程式ソルバと呼ばれるものに相当する⁸. そこでこのodeソルバを例にとって以下では説明する.

ODEWorldの準備

まずはWorldを作る.

w = ode.ODEWorld(Real3(1, 1, 1))

'Real3'は3つの実数値をまとめたもので3次元ベクトルに相当する. Worldは1つ目の引数として直方体の三辺の長さを受け取る. この例では計算を行う空間を一辺が1の立方体とした. 注意として, 'run_simulation'とは異なり, 容積のみを与えることはできない.

箱ができたので次は実際に分子を加えてみる.

w.add_molecules(Species('C'), 60)
print(w.t(), w.num_molecules(Species('C')))  # will return (0.0, 60)

分子を加えるときはadd_molecules, 減らすときはremove_molecules, その数を知りたければnum_moleculesを使えばよい. どれも1つ目の引数は問題とする分子種 Speciesである. 現在時刻は't'関数によって得られる.

ただし, odeソルバでは分子数は実数値として扱われるが, これらの関数では整数値しか扱うことができない. 実数を扱いたい場合はset_valueとget_valueを使う.

さて, これでWorldの準備はできた.

Real3の使い方

ここでSimulatorにうつる前に少し, 三次元ベクトルReal3の使い方について説明しておく.

pos = Real3(1, 2, 3)
print(pos)  # will print <ecell4.core.Real3 object at 0x7f44e118b9c0>
print(tuple(pos))  # will print (1.0, 2.0, 3.0)

Real3はそのままprintしても中身は見えないので, 一度Pythonのtupleかlistに変換してから出力しなければならない.

pos1 = Real3(1, 1, 1)
x, y, z = pos[0], pos[1], pos[2]
pos2 = pos1 + pos1
pos3 = pos1 * 3
pos4 = pos1 / 5
print(length(pos1))  # will print 1.73205080757
print(dot_product(pos1, pos3))  # will print 9.0

など基本的な関数は利用できるが, もちろん

a = numpy.asarray(tuple(Real3(1, 2, 3)))
print(a)  # will print [ 1.  2.  3.]

のようにしてnumpyのarrayにするのも良いだろう.

ODESimulatorの作成と実行

SimulatorはModelとWorldを与えることで作れる.

sim = ode.ODESimulator(m, w)
sim.run(10.0)

あとは'run'関数を呼べば計算が行われる. 上の例では10秒分の計算を実行する. Worldの状態がどう変わったかをみてみよう. 計算が進み, Cの数が減っていることがわかる.

print(w.t(), w.num_molecules(Species('C')))  # will return (10.0, 30)

Worldはあくまでも現在もしくはある時点での状態を表すものであるため, これだけでは途中の経過を追うことができない.

時系列を得るためにはObserverを使う.

obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)
print(obs.data())  # will return [[0.0, 0.0, 60.0], ..., [10.0, 29.994446899698026, 30.005553100301967]]

Observerには様々な種類があるが, 時系列を得るための簡単なものとしてFixedIntervalNumberObserverがある. その名の通り, 決まった時間刻みで分子の数を記録しておいてくれる. 1つ目の引数に時間幅を与え, 2つ目の引数で分子種を指定する⁹. その結果は'data'関数で引き出すことができるが, もっと簡単に

viz.plot_number_observer(obs)

とすることで時系列をプロットすることもできる.

以上で'run_simulation'関数内で行われている処理について説明が済んだ.

最後に, Simulatorを作った後にWorldに手を加えた場合はそのことをSimulatorに報せる必要があるため, 'sim.initialize()'を呼ぶことを忘れてはいけない.

ソルバを切り替える

より詳しくは後述するが, 'run_simulation'で試した通り, 確率論的手法に切り替えることはさほど難しくない.

equilibrium2.py

from ecell4 import *

with reaction_rules():
    A + B == C | (0.01, 0.3)

m = get_model()

# ode.ODEWorld -> gillespie.GillespieWorld
w = gillespie.GillespieWorld(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)

# ode.ODESimulator -> gillespie.GillespieSimulator
sim = gillespie.GillespieSimulator(m, w)
obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)

viz.plot_number_observer(obs)

WorldやSimulatorはModel自体を変更することはないため, 1つのモデルで様々な手法や条件を変えて計算を行うことも可能だ.

練習問題

Modelに関する練習: モデルといえば振動モデルが華だが, 今回は質量作用則以外の反応式は使えない. 質量作用則だけを用いて振動モデルを作るためにはいくつの反応が必要だろうか? できるだけ少ない反応数で作る方法はないか?
run_simulationに関する練習: gillespieソルバは確率論的手法と呼ばれているが, 分子の数が多くなれば, 決定論的な結果に近づくと考えられる. 実際に容積や分子数を増やしてみて, その結果がどうなるか確認する.
WorldとSimulatorに関する練習: 'run_simulation'関数では連続した計算しかできない. WorldとSimulatorを使って, 計算を5秒間行った後, 分子の数を変えて再び5秒の計算を行ってみる. また, その結果をObserverを用いてプロットしてみよう. ただし, Worldの状態を変えた場合はSimulatorを'initialize()'しなければならないので注意.
Observerに関する練習: 'viz.plot_number_observer'関数は複数のNumberObserverを受け取って重ねてプロットすることができる. そこで, 1つのモデルを作成し, odeとgillespieのそれぞれでシミュレーションを行い, 結果を記録した2つのObserverを同時に表示させて比較する.

詳しくは後述するが, 空白を含んだり, 数字や'_'からはじまってはいけない. ’A_B’などはよろしい. ↩
似た関数にlist_speciesがあるが, species_attributesとは違いモデルに含まれるSpecies全てをリスト化して返す. この違いについては後述する. ↩
数値だけ与えると自動的に100分割して記録される. 自分の好きな時刻に記録したければリストで与えるか, numpy.linspaceなどを利用すると良い. http://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html ↩
Gillespie D.T., "Exact stochastic simulation of coupled chemical reactions", J. Phys. Chem., 81(25), 2340–2361 (1977) ↩
http://docs.python.jp/2/tutorial/datastructures.html#tut-tuples ↩
必ず'()'か'[]'を使うこと. RuntimeErrorとなってしまう. ↩
この場合は, '~'の後の分子種は無視されるので'A'でなくて何でも良いのだが, 存在する分子種を書いた方が都合が良い. ↩
rosenbrock4. http://headmyshoulder.github.io/odeint-v2/doc/boost_numeric_odeint/odeint_in_detail/steppers.html#boost_numeric_odeint.odeint_in_detail.steppers.stepper_overview ↩
分子種はSpeciesではなく文字列で与えることに注意. ↩

E-Cell4の簡単な使い方 初級編

準備

使い方