8. 「1. E-Cell4を用いたシミュレーションの概要」についてより詳しく

一度1. E-Cell4を用いたシミュレーションの概要を読んだ後は、World と Simulator を使うのは難しくありません。 volume と {『C』 : 60} は World に相当し、ソルバーは以下の Simulator に相当します。

%matplotlib inline
from ecell4 import *

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

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

ここでは run_simulation がその内部で行っていることを分解して示します。 run_simulation はデフォルトでODEシミュレータを使用するので、ODEWorld を段階的に作成してみましょう。

8.1. ODEWorld の作成

World はこのようにして作ることができます。

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

Real3 は座標ベクトルです。

この例では、 ODEWorld コンストラクタの最初の引数はキューブです。

run_simulation 引数のように、ode.ODEWorld 引数にはvolumeを使用できないことに注意してください。

今度はシミュレーション用のキューブボックスを作成した後、分子をキューブに投げ込みましょう。

w = ode.ODEWorld(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
print(w.t(), w.num_molecules(Species('C')))  # must return (0.0, 60)

(0.0, 60)

分子を追加するには add_molecules、分子を削除するにはremove_molecules、分子の数を知るには num_molecules を使います。

各メソッドの第一引数は、あなたが知りたい Species です。

t メソッドで現在時刻を取得できます。

しかし、ODEソルバーの分子数は実数であり、これらの _molecules 関数では整数に対してのみ動作します。

ODEで実数を扱うときは、 set_value と get_value を使います。

8.2. Real3 の使い方

Simulator の詳細の前に、Real3 について詳しく説明します。

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

<ecell4.core.Real3 object at 0x7f05b0e6f630>
(1.0, 2.0, 3.0)

Real3 オブジェクトの内容を直接 print することはできません。

Real3 をPythonタプルまたはリストに一度変換する必要があります。

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))  # must print 1.73205080757
print(dot_product(pos1, pos3))  # must print 9.0

1.73205080757
9.0

dot_product のような基本的な関数も使えます。

もちろん、Real3 をnumpyの配列に変換することもできます。

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

[ 1.  2.  3.]

Integer3 は整数のトリプレットになります。

g = Integer3(1, 2, 3)
print(tuple(g))

(1, 2, 3)

もちろん、単純な算術演算を Integer3 に適用することもできます。

print(tuple(Integer3(1, 2, 3) + Integer3(4, 5, 6)))  # => (5, 7, 9)
print(tuple(Integer3(4, 5, 6) - Integer3(1, 2, 3)))  # => (3, 3, 3)
print(tuple(Integer3(1, 2, 3) * 2))  # => (2, 4, 6)
print(dot_product(Integer3(1, 2, 3), Integer3(4, 5, 6)))  # => 32
print(length(Integer3(1, 2, 3)))  # => 3.74165738677

(5, 7, 9)
(3, 3, 3)
(2, 4, 6)
32
3.74165738677

8.3. ODESimulator の作成と実行

下記のように Model と World を使って Simulatorを作成することができます。

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()

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

run メソッドを呼び出すと、シミュレーションが実行されます。

この例では、シミュレーションは10秒間実行されます。

また World の状態は下記のようにしてチェックできます。

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

(10.0, 30)

Species C の数が60から30に減少していることがわかります。

World はある時点での状態を表しているため、World でシミュレーションの遷移を見ることはできません。

時系列の結果を得るには、 Observer を使います。

w = ode.ODEWorld(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
sim = ode.ODESimulator(m, w)

obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)
print(obs.data())  # must return [[0.0, 0.0, 60.0], ..., [10.0, 29.994446899691276, 30.005553100308752]]

[[0.0, 0.0, 60.0], [0.1, 1.7722206098711988, 58.227779390128795], [0.2, 3.4860124889661757, 56.51398751103382], [0.30000000000000004, 5.1376332715496495, 54.862366728450354], [0.4, 6.724090809612153, 53.27590919038786], [0.5, 8.243129756755453, 51.75687024324456], [0.6000000000000001, 9.69320376680592, 50.3067962331941], [0.7000000000000001, 11.073435590968808, 48.92656440903121], [0.8, 12.383567691608423, 47.616432308391595], [0.9, 13.62390591657045, 46.37609408342957], [1.0, 14.795258681171735, 45.20474131882829], [1.1, 15.898873899780316, 44.10112610021971], [1.2000000000000002, 16.936375633755194, 43.06362436624483], [1.3, 17.90970211303103, 42.090297886969], [1.4000000000000001, 18.821046466966358, 41.17895353303366], [1.5, 19.67280118133671, 40.32719881866331], [1.6, 20.46750699887015, 39.53249300112988], [1.7000000000000002, 21.207806714875233, 38.792193285124796], [1.8, 21.896404094818397, 38.10359590518163], [1.9000000000000001, 22.536027939969333, 37.463972060030684], [2.0, 23.129401186565122, 36.8705988134349], [2.1, 23.679214801318086, 36.320785198681946], [2.2, 24.188106157761755, 35.81189384223828], [2.3000000000000003, 24.658641522840725, 35.34135847715931], [2.4000000000000004, 25.093302252968826, 34.9066977470312], [2.5, 25.49447428958285, 34.50552571041718], [2.6, 25.86444054777012, 34.13555945222991], [2.7, 26.205375806607968, 33.79462419339207], [2.8000000000000003, 26.51934373464207, 33.48065626535796], [2.9000000000000004, 26.80829570797557, 33.19170429202446], [3.0, 27.07407111744664, 32.92592888255339], [3.1, 27.318398885233666, 32.68160111476636], [3.2, 27.542899949227504, 32.457100050772524], [3.3000000000000003, 27.749090501388586, 32.25090949861144], [3.4000000000000004, 27.938385796523626, 32.0616142034764], [3.5, 28.112104375214862, 31.88789562478517], [3.6, 28.271472567998387, 31.728527432001645], [3.7, 28.417629170003387, 31.582370829996645], [3.8000000000000003, 28.551630196373424, 31.44836980362661], [3.9000000000000004, 28.674453642480778, 31.32554635751925], [4.0, 28.787004191633095, 31.212995808366934], [4.1000000000000005, 28.890117821767866, 31.109882178232162], [4.2, 28.984566277309003, 31.015433722691025], [4.3, 29.071061377146354, 30.928938622853675], [4.4, 29.150259141926327, 30.8497408580737], [4.5, 29.22276372627192, 30.777236273728107], [4.6000000000000005, 29.289131148912805, 30.710868851087223], [4.7, 29.349872817459094, 30.650127182540935], [4.800000000000001, 29.405458847152726, 30.594541152847302], [4.9, 29.456321176884504, 30.543678823115524], [5.0, 29.50285648638057, 30.49714351361946], [5.1000000000000005, 29.545428921493947, 30.45457107850608], [5.2, 29.58437263404594, 30.415627365954087], [5.300000000000001, 29.61999414522251, 30.38000585477752], [5.4, 29.652574540327944, 30.347425459672085], [5.5, 29.682371504475235, 30.317628495524794], [5.6000000000000005, 29.70962120751236, 30.29037879248767], [5.7, 29.734540047958276, 30.265459952041752], [5.800000000000001, 29.757326264061593, 30.242673735938435], [5.9, 29.778161421010864, 30.221838578989164], [6.0, 29.797211782455758, 30.20278821754427], [6.1000000000000005, 29.814629574218237, 30.18537042578179], [6.2, 29.8305541480645, 30.169445851935528], [6.300000000000001, 29.84511305246321, 30.15488694753682], [6.4, 29.858423017246, 30.141576982754028], [6.5, 29.87059085870597, 30.129409141294058], [6.6000000000000005, 29.881714310971407, 30.11828568902862], [6.7, 29.89188278945169, 30.10811721054834], [6.800000000000001, 29.90117809152673, 30.098821908473298], [6.9, 29.90967503947212, 30.09032496052791], [7.0, 29.917442070090726, 30.082557929909303], [7.1000000000000005, 29.92454177537786, 30.075458224622167], [7.2, 29.931031398107503, 30.068968601892525], [7.300000000000001, 29.936963286002193, 30.063036713997835], [7.4, 29.942385307810326, 30.057614692189702], [7.5, 29.947341234455514, 30.052658765544514], [7.6000000000000005, 29.951871088073982, 30.048128911926046], [7.7, 29.956011461640173, 30.043988538359855], [7.800000000000001, 29.959795811548137, 30.04020418845189], [7.9, 29.963254725453346, 30.036745274546682], [8.0, 29.966416167388665, 30.033583832611363], [8.1, 29.969305702114955, 30.030694297885073], [8.200000000000001, 29.971946700366555, 30.028053299633473], [8.3, 29.97436052665158, 30.02563947334845], [8.4, 29.9765667110567, 30.023433288943327], [8.5, 29.978583106421084, 30.021416893578944], [8.6, 29.98042603207785, 30.019573967922177], [8.700000000000001, 29.982110405343864, 30.017889594656165], [8.8, 29.983649861761148, 30.01635013823888], [8.9, 29.985056865065513, 30.014943134934516], [9.0, 29.98634280775227, 30.01365719224776], [9.1, 29.9875181030258, 30.01248189697423], [9.200000000000001, 29.988592268884133, 30.011407731115895], [9.3, 29.98957400499266, 30.010425995007367], [9.4, 29.990471262976794, 30.009528737023235], [9.5, 29.99129131068186, 30.00870868931817], [9.600000000000001, 29.992040790927835, 30.007959209072194], [9.700000000000001, 29.99272577521901, 30.007274224781018], [9.8, 29.993351812847337, 30.00664818715269], [9.9, 29.993923975779538, 30.00607602422049], [10.0, 29.994446899691276, 30.005553100308752]]

E-Cell4にはいくつかのタイプの Observer があります。

FixedIntervalNumberObserver は、時系列の結果を得るための最もシンプルな Observer です。

その名前が示唆するように、この Observer は各時間ステップの分子の数を記録します。

第1引数は時間ステップ、第2引数は分子種です。

その結果は data メソッドで確認できますが、これにはショートカットがあります。

viz.plot_number_observer(obs)

上記は時系列の結果を簡単にプロットします。

ここでは run_simulation 関数の内部について説明しました。

Simulator を作成した後に World を変更した場合は Simulator にそれを示す必要があります。

またそのあとで sim.initialize() を呼び出すことを忘れないでください。

8.4. ソルバーの差し替えについて

run_simulation の内部を上記で示したので、ソルバーを確率論的手法に替えることは難しくないでしょう。

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 自体を決して変更しないので、一つの Model に対して 複数の Simulator を差し替えることができます。