!pip install ecell

Collecting ecell
  Downloading ecell-4.1.2-cp36-cp36m-manylinux1_x86_64.whl (45.5MB)
    100% |################################| 45.5MB 29kB/s  eta 0:00:01
Installing collected packages: ecell
Successfully installed ecell-4.1.2

5. シミュレーションを記録して視覚化する方法

ここでは、シミュレーション結果のログを取る方法と、それを視覚化する方法について説明します。

%matplotlib inline
import math
from ecell4 import *

5.1. Observer を用いたシミュレーションのログの取り方

E-Cell4は、 Observerというロギングのためのクラスを提供しています。 ObserverクラスはSimulatorの run関数を呼び出すときに与えられます。

def create_simulator(f=gillespie.GillespieFactory()):
    m = NetworkModel()
    A, B, C = Species('A', '0.005', '1'), Species('B', '0.005', '1'), Species('C', '0.005', '1')
    m.add_reaction_rule(create_binding_reaction_rule(A, B, C, 0.01))
    m.add_reaction_rule(create_unbinding_reaction_rule(C, A, B, 0.3))
    w = f.create_world()
    w.bind_to(m)
    w.add_molecules(C, 60)
    sim = f.create_simulator(w)
    sim.initialize()
    return sim

最も一般的な Observerの1つはFixedIntervalNumberObserverです。これは与えられた時間間隔で分子の数を記録します。 FixedIntervalNumberObserverは、時間間隔と、ログの対象とするのSpeciesのシリアルのリストを必要とします。

obs1 = FixedIntervalNumberObserver(0.1, ['A', 'B', 'C'])
sim = create_simulator()
sim.run(1.0, obs1)

FixedIntervalNumberObserverのdata関数はログに記録されたデータを返します。

print(obs1.data())

[[0.0, 0.0, 0.0, 60.0], [0.1, 1.0, 1.0, 59.0], [0.2, 3.0, 3.0, 57.0], [0.30000000000000004, 4.0, 4.0, 56.0], [0.4, 6.0, 6.0, 54.0], [0.5, 8.0, 8.0, 52.0], [0.6000000000000001, 9.0, 9.0, 51.0], [0.7000000000000001, 10.0, 10.0, 50.0], [0.8, 10.0, 10.0, 50.0], [0.9, 10.0, 10.0, 50.0], [1.0, 11.0, 11.0, 49.0]]

targets()はコンストラクタの引数として指定した Speciesのリストを返します。

print([sp.serial() for sp in obs1.targets()])

['A', 'B', 'C']

NumberObserverは、反応が起こったときのすべてのステップの後に分子の数を記録します。 このオブザーバーはすべての反応を記録するのに便利ですが、 odeでは利用できません。

obs1 = NumberObserver(['A', 'B', 'C'])
sim = create_simulator()
sim.run(1.0, obs1)
print(obs1.data())

[[0.0, 0.0, 0.0, 60.0], [0.043181315145277815, 1.0, 1.0, 59.0], [0.0999354928781079, 2.0, 2.0, 58.0], [0.10051642633328928, 3.0, 3.0, 57.0], [0.1576666300428589, 4.0, 4.0, 56.0], [0.22160562023979907, 5.0, 5.0, 55.0], [0.24838212056833223, 6.0, 6.0, 54.0], [0.38083606360319644, 7.0, 7.0, 53.0], [0.44115742167963257, 8.0, 8.0, 52.0], [0.478654279832333, 9.0, 9.0, 51.0], [0.524776672724245, 10.0, 10.0, 50.0], [0.5831861079171226, 11.0, 11.0, 49.0], [0.6016715021229188, 12.0, 12.0, 48.0], [0.6458184477531915, 13.0, 13.0, 47.0], [0.7164454853598512, 14.0, 14.0, 46.0], [0.8389493519306314, 15.0, 15.0, 45.0], [0.9250934945924608, 16.0, 16.0, 44.0], [0.9314408219193694, 17.0, 17.0, 43.0], [1.0, 17.0, 17.0, 43.0]]

TimingNumberObserverでは、ログの時刻をそのコンストラクタの引数として与えることができます。

obs1 = TimingNumberObserver([0.0, 0.1, 0.2, 0.5, 1.0], ['A', 'B', 'C'])
sim = create_simulator()
sim.run(1.0, obs1)
print(obs1.data())

[[0.0, 0.0, 0.0, 60.0], [0.1, 1.0, 1.0, 59.0], [0.2, 3.0, 3.0, 57.0], [0.5, 8.0, 8.0, 52.0], [1.0, 13.0, 13.0, 47.0]]

run関数は複数のオブザーバを一度に受け付けます。

obs1 = NumberObserver(['C'])
obs2 = FixedIntervalNumberObserver(0.1, ['A', 'B'])
sim = create_simulator()
sim.run(1.0, [obs1, obs2])
print(obs1.data())
print(obs2.data())

[[0.0, 60.0], [0.05163820372508536, 59.0], [0.10222920269451034, 58.0], [0.14249526069798346, 57.0], [0.1709401154327075, 56.0], [0.1786522132631355, 55.0], [0.32345579486892295, 54.0], [0.3320212006041937, 53.0], [0.37295796722048336, 52.0], [0.4297993097403142, 51.0], [0.442315663109651, 50.0], [0.456452361786781, 49.0], [0.759999628617218, 48.0], [0.7601731295722447, 47.0], [0.847168492783278, 46.0], [0.856013633150855, 47.0], [0.8614538661854448, 46.0], [0.870506834401381, 45.0], [0.978711403259659, 44.0], [1.0, 44.0]]
[[0.0, 0.0, 0.0], [0.1, 1.0, 1.0], [0.2, 5.0, 5.0], [0.30000000000000004, 5.0, 5.0], [0.4, 8.0, 8.0], [0.5, 11.0, 11.0], [0.6000000000000001, 11.0, 11.0], [0.7000000000000001, 11.0, 11.0], [0.8, 13.0, 13.0], [0.9, 15.0, 15.0], [1.0, 16.0, 16.0]]

FixedIntervalHDF5ObservedrはWorldのデータ全体を一定の間隔で出力ファイルに記録します。 2番目の引数は、出力ファイル名の接頭辞です。 filename()は、次に保存するようにスケジュールされたファイルの名前を返します。 ％02dのような多くのフォーマット文字列では、ファイル名にステップカウントを使用することができます。 書式文字列を使用しないと、最新のデータがファイルに上書きされます。

obs1 = FixedIntervalHDF5Observer(0.2, 'test%02d.h5')
print(obs1.filename())
sim = create_simulator()
sim.run(1.0, obs1)
print(obs1.filename())

test00.h5
test06.h5

w = load_world('test05.h5')
print(w.t(), w.num_molecules(Species('C')))

1.0 44

FixedIntervalCSVObserverの使用方法は、FixedIntervalHDF5Observerの使用方法とほぼ同じです。 FixedIntervalCSVObserverは半径(r)とSpeciesのシリアル番号(sid)を持つパーティクルの位置(x、y、z)をCSVファイルに保存します。

obs1 = FixedIntervalCSVObserver(0.2, "test%02d.csv")
print(obs1.filename())
sim = create_simulator()
sim.run(1.0, obs1)
print(obs1.filename())

test00.csv
test06.csv

下記は出力CSVファイルの最初の10行です。

print(''.join(open("test05.csv").readlines()[: 10]))

x,y,z,r,sid
0.17508278810419142,0.12197625380940735,0.53627523058094084,0,0
0.24078186159022152,0.62765785818919539,0.7254820850212127,0,0
0.40445487247779965,0.16720660566352308,0.92957371729426086,0,0
0.24083335651084781,0.59580284473486245,0.51183571317233145,0,0
0.90290508698672056,0.4137285016477108,0.39729581261053681,0,0
0.51044316007755697,0.26319809840060771,0.40705726365558803,0,0
0.84177179471589625,0.085374288959428668,0.64473599148914218,0,0
0.54177056765183806,0.098819603677839041,0.17240164172835648,0,0
0.28669063560664654,0.71829434810206294,0.91102162329480052,0,0

粒子シミュレーションのために、E-Cell4は、分子の軌道を追跡するためのオブザーバ（Observer）を提供し、これはFixedIntervalTrajectoryObserverと名付けられています。 ParticleIDが特に指定されていないときは、すべての軌跡を記録します。 シミュレーション中にパーティクルIDが反応のために失われると、FixedIntervalTrajectoryObserverはパーティクルをそれ以上トレースするのを単に止めます。

sim = create_simulator(spatiocyte.SpatiocyteFactory(0.005))
obs1 = FixedIntervalTrajectoryObserver(0.01)
sim.run(0.1, obs1)

print([tuple(pos) for pos in obs1.data()[0]])

[(0.46540305112880387, 0.739008344562721, 0.48), (0.2122891110412088, 0.9266471820493494, 0.665), (0.49806291436591293, 0.8689121551303868, 0.685), (0.7838367176906171, 0.8313843876330611, 0.64), (0.9961258287318259, 0.7967433714816836, 0.5), (0.7920016834998943, 0.9122134253196088, 0.63), (0.7348469228349536, 0.7534421012924616, 0.515), (0.9716309313039941, 0.851591647054698, 0.525), (1.1349302474895393, 0.9122134253196088, 0.63), (0.9961258287318259, 0.9439676901250381, 0.725), (0.718516991216399, 0.8833459118601275, 0.93)]

ほとんどの場合、Worldは各平面の周期的境界を想定しています。 境界条件に起因するエッジでのパーティクルの大きなジャンプを避けるために、 FixedIntervalTrajectoryObserverは位置のシフトを維持しようとします。 従って、Observerに格納されたパーティクルの位置は、Worldに対して与えられた立方体に必ずしも限定されません。 境界条件の拡散を正確に追跡するには、ロギングのステップ間隔を十分小さくする必要があります。 もちろん、このオプションは無効にすることができます。 help(FixedIntervalTrajectoryObserver)を参照してください。

5.2. ログデータの可視化

このセクションでは、Observerによって記録されたログデータの可視化ツールについて説明します。

まず、時系列データのためにviz.plot_number_observerはNumberObserverとFixedIntervalNumberObserverおよびTimingNumberObserverによって提供されるデータをプロットします。

viz.plot_number_observerの詳しい使用法については、help(viz.plot_number_observer)を参照してください。

obs1 = NumberObserver(['C'])
obs2 = FixedIntervalNumberObserver(0.1, ['A', 'B'])
sim = create_simulator()
sim.run(10.0, [obs1, obs2])

viz.plot_number_observer(obs1, obs2)

プロットのスタイルを設定したり、プロットに任意の関数を追加することもできます。

viz.plot_number_observer(obs1, '-', obs2, ':', lambda t: 60 * math.exp(-0.3 * t), '--')

位相平面でのプロットは、x軸とy軸を指定することによっても利用できます。

viz.plot_number_observer(obs2, 'o', x='A', y='B')

空間シミュレーションでは、Worldの状態を可視化するために、viz.plot_worldが利用可能です。 この関数は、インタラクティブな方法で3次元ボリューム内のパーティクルの点をプロットします。 描画領域の右ボタンをクリックすると、画像を保存できます。

sim = create_simulator(spatiocyte.SpatiocyteFactory(0.005))
# viz.plot_world(sim.world())
viz.plot_world(sim.world(), interactive=False)

FixedIntervalHDF5Observerとして与えられた一連のHDF5ファイルからムービーを作ることもできます。

注： viz.plot_movieは、オプションinteractive=Falseの時に追加ライブラリ ffmpegを必要とします。

sim = create_simulator(spatiocyte.SpatiocyteFactory(0.005))
obs1 = FixedIntervalHDF5Observer(0.02, 'test%02d.h5')
sim.run(1.0, obs1)
viz.plot_movie(obs1)

最後に、 FixedIntervalTrajectoryObserverに対応して、viz.plot_trajectoryはパーティクル軌道の可視化を提供します。

sim = create_simulator(spatiocyte.SpatiocyteFactory(0.005))
obs1 = FixedIntervalTrajectoryObserver(1e-3)
sim.run(1, obs1)
# viz.plot_trajectory(obs1)
viz.plot_trajectory(obs1, interactive=False)