海津一成
In [1]:
%matplotlib inline
from ecell4 import *
from ecell4_base.core import *
from ecell4_base import *
In [2]:
import seaborn
seaborn.set(font_scale=1.5)
In [3]:
import matplotlib as mpl
# mpl.rc("figure", figsize=(8.0, 5.5))
mpl.rc("figure", figsize=(6, 4))
Model: いわゆるモデル。分子の特徴量（拡散定数、大きさなど）や反応を記述したもの。具体的にはNetworkModelを使う。World: ある時刻における状態をあらわすもの。分子の数や位置、空間の表現法など。使うWorldはアルゴリズムによって異なる. 例えば、ode.Worldみたいな感じ。ode.Simulatorみたいな感じ。Modelは扱う問題を記述したもので、大きく2つの要素によって構成されている。
Species: すなわち分子種。モデルに登場する実体をあらわす。ReactionRule: すなわち反応則。分子種（Species）間の相互作用をあらわす。E-Cell4では主に化学反応を扱うので、それを例に言えば
Species: 基質や酵素などどんな種類の分子が存在し、それらがどのような性質 (分子の大きさなど)を持つかをあらわす。ReactionRule: 結合やかい離などそれら分子間にどのような反応があるかをあらわす。
In [4]:
A = Species('A')
B = Species('B')
とすれば良い。
それぞれ、AとBという名前の分子種をあらわす。1個の特定の分子ではなく、そういう種類の分子があることを意味している。
Speciesには任意の属性を持たせることで、どのような性質を持つのかを指定できる。
In [5]:
A = Species("A")
A.set_attribute("radius", 0.005)
A.set_attribute("D", 1)
A.set_attribute("location", "cytoplasm")
1つ目に与える引数が属性名で、2つ目がその値で、属性名は文字列で値は文字列か数値、真偽値でなければならない。
属性のうち特に分子径、拡散係数、属する場所は頻繁に用いられるので簡単に
In [6]:
A = Species("A", 0.005, 1, "cytoplasm") # XXX: serial, radius, D, location
でまとめて与えることができる。
与えたSpeciesの属性を知りたいときは
In [7]:
print(A.serial())
print(A.get_attribute("radius").magnitude)
print(A.get_attribute("D").magnitude)
print(A.get_attribute("location"))
とすれば良い。数値属性は単位を伴ったQuantityとして返されるのでmagnitudeによって値だけ取り出せる。
In [8]:
del A, B
Speciesができたところで今度は反応を作ってみる。
反応はおおざっぱに言って、何かの分子種が何かの分子種になること。
前者をReactants、後者をProductsと呼ぶ。あとはその反応の速さ（速度定数）を数値で与えれば良い。
In [9]:
rr = ReactionRule()
rr.add_reactant(Species("A"))
rr.add_product(Species("B"))
rr.set_k(1.0)
これで分子種Aが分子種Bに変る反応を作れた。ここでは反応を定義しているので、Speciesに属性値はなくても良い。
確認してみる。与えられたReactionRuleがどのようなものかを簡単に知りたい場合はas_string関数で見られる。
In [10]:
print(rr.as_string())
結合もこれと同じに記述できる。AとBの結合によってCができる。
In [11]:
rr = ReactionRule()
rr.add_reactant(Species("A"))
rr.add_reactant(Species("B"))
rr.add_product(Species("C"))
rr.set_k(1.0)
print(rr.as_string())
より読み手に伝わりやすい書き方として以下の便利関数も用意されている。
In [12]:
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_unbinding_reaction_rule(
Species("C"), Species("A"), Species("B"), 1.0)
以上でModelに含まれる要素を作れるようになった。
In [13]:
sp1 = Species("A", 0.005, 1.0)
sp2 = Species("B", 0.005, 1.0)
sp3 = Species("C", 0.010, 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)
あとはこれらをModelに追加していく。
Speciesの追加はadd_species_attributeで、ReactionRuleの追加はadd_reaction_ruleで行う。
In [14]:
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_attributesとreaction_rulesを使う。list_speciesはまた別物なので注意！
In [15]:
print(m.species_attributes())
print(m.reaction_rules())
ちなみに、特にSpeciesの属性を必要としない計算(常微分方程式やGillespie法)ではadd_species_attributeは省略しても良い。
WorldとSimulatorの説明がまだだが、早くシミュレーションを走らせたい人のためにここで一度計算を行なってみる。
In [16]:
# XXX: 'm' is a NetworkModel, which is described in the section above.
run_simulation(10, {'C': 60}, volume=1.0, model=m)
run_simulationは分子数を記録して時系列をプロットしてくれる。
Cを60分子とした。volumeは計算する空間の容積。modelはモデル。先程のNetworkModelを与えた。この味気ないプロットに満足できない人のためには確率論的計算手法というものがある。
In [17]:
run_simulation(10, {'C': 60}, volume=1.0, model=m, solver='gillespie')
全く同じモデルを異る手法で計算することができた。これは問題と計算手法が分離されているからこそなせる技である。
E-Cell4を使ってSpeciesとReactionRuleを作り、NetworkModelに登録して簡単に計算することができるようになった。
しかし、モデルの記述はこのように簡単な結合とかい離だけであっても煩雑でモデルの読み手に伝わりづらい。
そこで、E-Cell4ではこのモデルを記述するための独自の構文をいくつか勝手に用意している。
まずはReactionRuleについて見てみると、
In [18]:
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つだけ与える。また、可逆反応については特別な記法があり、それを使えばより簡単になる。
In [19]:
with reaction_rules():
A + B == C | (0.01, 0.3)
run_simulation(10, {'C': 60})
可逆反応の記号は==で、セパレータの後にはtupleを用いて2つの速度定数を与える。
run_simulation関数は引数modelを指定しない場合は自分でget_model()を呼んでモデルとするので、このように省略可能な場合がある。
In [20]:
with reaction_rules():
A > | 1.0 # XXX: will throw SyntaxError
> A | 1.0 # XXX: will throw SyntaxError
これらを表現するために特殊記号~を用いる。
ReactionRuleを記述する際にSpeciesの前に~記号を置くと反応上のその分子種の量論係数が0になる。簡単に言えば、単に無いものとみなされる。
In [21]:
with reaction_rules():
A > ~A | 1.0 # XXX: create_degradation_reaction_rule
~A > A | 1.0 # XXX: create_synthesis_reaction_rule
In [22]:
m = get_model()
Modelが理解できればWorldとSimulatorはさほど難しくはない。
さきほどのrun_simulation関数で与えたvolumeと{'C': 60}などがWorldで、solverがSimulatorに相当する。
実際にrun_simulationがしている作業を自分で書けるようにする。
run_simulationは標準で常微分方程式ソルバodeを使う。
In [23]:
w = ode.World(Real3(1, 1, 1))
Real3は3つの実数値をまとめた3次元ベクトル。ones()としてもReal3(1, 1, 1)が得られる。Worldは1つ目の引数として直方体の三辺の長さを受け取る。この例では一辺が1の立方体とした。
In [24]:
print(w.t())
print(w.edge_lengths())
print(w.volume())
t関数で得られる。edge_lengths関数。volume関数。箱ができたので次は実際に分子を加えてみる.
In [25]:
w.add_molecules(Species('C'), 60)
print(w.num_molecules(Species('C')))
add_molecules。remove_molecules。num_molecules。どれも1つ目の引数は問題とする分子種Species。
ただし、odeソルバでは分子数は実数値として扱われるが、これらの関数では整数値しか扱うことができない。
実数を扱いたい場合はset_valueとget_valueを使う。
In [26]:
with reaction_rules():
A + B == C | (0.01, 0.3)
m = get_model()
In [27]:
sim = ode.Simulator(w, m)
あとはrun関数を呼べば計算が行われる。
In [28]:
sim.run(10)
Worldの状態がどう変わったかをみてみよう。計算が進み、Cの数が減っている。
In [29]:
print(w.t(), w.num_molecules(Species('C')))
Worldはあくまでも現在もしくはある時点での状態を表すもので、これだけでは途中の経過を追うことができない。
時系列を得るためにはObserverを使う。
In [30]:
obs = FixedIntervalNumberObserver(1, ('A', 'B', 'C'))
sim.run(10, obs)
print(obs.data())
Observerには様々な種類があるが、時系列を得るにはFixedIntervalNumberObserverを使う。
決まった時間幅で分子数を記録する。1つ目の引数が時間幅で、2つ目の引数が分子種。その結果はdata関数で引き出せる。
viz.plot_number_observer関数で簡単に時系列をプロットすることもできる（単にshowとしても良い）。
In [31]:
w = ode.World(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
sim = ode.Simulator(w, m)
obs = FixedIntervalNumberObserver(0.1, ('A', 'B', 'C'))
sim.run(10, obs)
viz.plot_number_observer(obs)
In [32]:
# ode.World -> gillespie.World
w = gillespie.World(ones())
w.add_molecules(Species('C'), 60)
# ode.Simulator -> gillespie.Simulator
sim = gillespie.Simulator(w, m)
obs = FixedIntervalNumberObserver(0.1, ('A', 'B', 'C'))
sim.run(10, obs)
show(obs)
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を同時に表示させて比較する。