In [1]:
!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:011 day, 23:59:59 68% |##################### | 31.0MB 52.8MB/s eta 0:00:01 91% |############################# | 41.9MB 51.9MB/s eta 0:00:01 98% |############################### | 44.8MB 59.9MB/s eta 0:00:01
Installing collected packages: ecell
Successfully installed ecell-4.1.2
4. シミュレーションの走らせ方¶
セクション2と3では、モデルを構築し、初期状態を設定する方法について説明しました。
それでは、シミュレーションを実行してみましょう。
「Spatiocyte.SpatiocyteSimulator」、「egfrd.EGFRDSimulator」、「bd.BDSimulator」、「meso.MesoscopicSimulator」、「gillespie.GillespieSimulator」、および「ode.ODESimulator」の6つの「シミュレータ」クラスが存在します。
それぞれの
Simulator
クラスはWorld
の対応する型だけを受け入れますが、それらのすべてが同じ
Model
を受け入れます。
In [2]:
import ecell4
4.1. シミュレーターのセットアップの方法¶
アルゴリズム固有の引数を持つ初期化(いわゆるコンストラクタ関数)を除いて、すべてのシミュレータは同じAPIを持っています。
In [3]:
from ecell4.core import *
from ecell4.gillespie import GillespieWorld, GillespieSimulator
from ecell4.ode import ODEWorld, ODESimulator
from ecell4.spatiocyte import SpatiocyteWorld, SpatiocyteSimulator
from ecell4.bd import BDWorld, BDSimulator
from ecell4.meso import MesoscopicWorld, MesoscopicSimulator
from ecell4.egfrd import EGFRDWorld, EGFRDSimulator
Simulator
を構築する前に、Simulator
のタイプに対応するModel
とWorld
を用意してください。
In [4]:
from ecell4 import species_attributes, reaction_rules, get_model
with species_attributes():
A | B | C | {'D': '1', 'radius': '0.005'}
with reaction_rules():
A + B == C | (0.01, 0.3)
m = get_model()
In [5]:
w1 = GillespieWorld()
w2 = ODEWorld()
w3 = SpatiocyteWorld()
w4 = BDWorld()
w5 = MesoscopicWorld()
w6 = EGFRDWorld()
シミュレータは、構築時に下記の順序でモデルとワールドの両方を必要とします。
In [5]:
sim1 = GillespieSimulator(m, w1)
sim2 = ODESimulator(m, w2)
sim3 = SpatiocyteSimulator(m, w3)
sim4 = BDSimulator(m, w4)
sim5 = MesoscopicSimulator(m, w5)
sim6 = EGFRDSimulator(m, w6)
Model
をWorld
にバインドした場合は、シミュレータを作成するために要るのはWorld
だけです。
In [6]:
w1.bind_to(m)
w2.bind_to(m)
w3.bind_to(m)
w4.bind_to(m)
w5.bind_to(m)
w6.bind_to(m)
In [7]:
sim1 = GillespieSimulator(w1)
sim2 = ODESimulator(w2)
sim3 = SpatiocyteSimulator(w3)
sim4 = BDSimulator(w4)
sim5 = MesoscopicSimulator(w5)
sim6 = EGFRDSimulator(w6)
もちろん、Simulator
にバインドされたModel
とWorld
は、Simulator
から以下のように引き出すことができます。
In [8]:
print(sim1.model(), sim1.world())
print(sim2.model(), sim2.world())
print(sim3.model(), sim3.world())
print(sim4.model(), sim4.world())
print(sim5.model(), sim5.world())
print(sim6.model(), sim6.world())
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.gillespie.GillespieWorld object at 0x7f6b329d1a08>)
(<ecell4.ode.ODENetworkModel object at 0x7f6b329d1af8>, <ecell4.ode.ODEWorld object at 0x7f6b329d1ae0>)
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.spatiocyte.SpatiocyteWorld object at 0x7f6b329d1a08>)
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.bd.BDWorld object at 0x7f6b329d1ae0>)
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.meso.MesoscopicWorld object at 0x7f6b329d1ac8>)
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.egfrd.EGFRDWorld object at 0x7f6b329d1ae0>)
World
を自分で更新した後は、シミュレーションを実行する前にSimulator
の内部状態を初期化する必要があります。
In [9]:
w1.add_molecules(Species('C'), 60)
w2.add_molecules(Species('C'), 60)
w3.add_molecules(Species('C'), 60)
w4.add_molecules(Species('C'), 60)
w5.add_molecules(Species('C'), 60)
w6.add_molecules(Species('C'), 60)
In [10]:
sim1.initialize()
sim2.initialize()
sim3.initialize()
sim4.initialize()
sim5.initialize()
sim6.initialize()
一定のステップ間隔を有するアルゴリズムの場合は上記に加えdt
も必要です。
In [11]:
sim2.set_dt(1e-6) # ODESimulator. This is optional
sim4.set_dt(1e-6) # BDSimulator
4.2. シミュレーションの実行¶
シミュレーションを実行するために、Simulatorはstep
とrun
の2つのAPIを提供します。
step()
は、シミュレータが期待する時間(next_time()
)の間シミュレーションを進めます。
In [12]:
print(sim1.t(), sim1.next_time(), sim1.dt())
print(sim2.t(), sim2.next_time(), sim2.dt()) # => (0.0, 1e-6, 1e-6)
print(sim3.t(), sim3.next_time(), sim3.dt())
print(sim4.t(), sim4.next_time(), sim4.dt()) # => (0.0, 1e-6, 1e-6)
print(sim5.t(), sim5.next_time(), sim5.dt())
print(sim6.t(), sim6.next_time(), sim6.dt()) # => (0.0, 0.0, 0.0)
(0.0, 0.032171880483933143, 0.032171880483933143)
(0.0, 1e-06, 1e-06)
(0.0, 1.6666666666666667e-05, 1.6666666666666667e-05)
(0.0, 1e-06, 1e-06)
(0.0, 0.029999537310460653, 0.029999537310460653)
(0.0, 0.0, 0.0)
In [13]:
sim1.step()
sim2.step()
sim3.step()
sim4.step()
sim5.step()
sim6.step()
In [14]:
print(sim1.t())
print(sim2.t()) # => 1e-6
print(sim3.t())
print(sim4.t()) # => 1e-6
print(sim5.t())
print(sim6.t()) # => 0.0
0.0321718804839
1e-06
1.66666666667e-05
1e-06
0.0299995373105
0.0
last_reactions()
は、最後のステップで発生するReactionRule
とReactionInfo
のペアのリストを返します。
各アルゴリズムにはReactionInfoの独自の実装があります。
詳細については、help(module.ReactionInfo)
を参照してください。
In [15]:
print(sim1.last_reactions())
# print(sim2.last_reactions())
print(sim3.last_reactions())
print(sim4.last_reactions())
print(sim5.last_reactions())
print(sim6.last_reactions())
[(<ecell4.core.ReactionRule object at 0x7f6b329d1ae0>, <ecell4.gillespie.ReactionInfo object at 0x7f6b329d1af8>)]
[]
[]
[(<ecell4.core.ReactionRule object at 0x7f6b329d1ae0>, <ecell4.meso.ReactionInfo object at 0x7f6b329d1af8>)]
[]
step(upto)
は、next_time
がupto
より小さい場合はnext_time
のシミュレーションを進め、それ以外の場合はupto
まで実行します。
またstep(upto)
は時間が上限(upto
)に達していないかどうかを返します。
In [16]:
print(sim1.step(1.0), sim1.t())
print(sim2.step(1.0), sim2.t())
print(sim3.step(1.0), sim3.t())
print(sim4.step(1.0), sim4.t())
print(sim5.step(1.0), sim5.t())
print(sim6.step(1.0), sim6.t())
(True, 0.11785771174833615)
(True, 2e-06)
(True, 3.3333333333333335e-05)
(True, 2e-06)
(True, 0.0709074540101994)
(True, 0.0)
upto
の時間ちょうどまでシミュレーションをは走らせるには,
それがTrue
を返す間 step(upto)
を呼びます。
In [17]:
while sim1.step(1.0): pass
while sim2.step(0.001): pass
while sim3.step(0.001): pass
while sim4.step(0.001): pass
while sim5.step(1.0): pass
while sim6.step(0.001): pass
In [18]:
print(sim1.t()) # => 1.0
print(sim2.t()) # => 0.001
print(sim3.t()) # => 0.001
print(sim4.t()) # => 0.001
print(sim5.t()) # => 1.0
print(sim6.t()) # => 0.001
1.0
0.001
0.001
0.001
1.0
0.001
これはちょうど run
の動作と同じものです。
run(tau)
はt()+tau
の時間までシミュレーションを進めます。
In [19]:
sim1.run(1.0)
sim2.run(0.001)
sim3.run(0.001)
sim4.run(0.001)
sim5.run(1.0)
sim6.run(0.001)
In [20]:
print(sim1.t()) # => 2.0
print(sim2.t()) # => 0.002
print(sim3.t()) # => 0.002
print(sim4.t()) # => 0.002
print(sim5.t()) # => 2.0
print(sim6.t()) # => 0.02
2.0
0.002
0.002
0.002
2.0
0.002
num_steps
は、シミュレーション中のステップ数を返します。
In [21]:
print(sim1.num_steps())
print(sim2.num_steps())
print(sim3.num_steps())
print(sim4.num_steps())
print(sim5.num_steps())
print(sim6.num_steps())
37
2001
120
2001
34
581
4.3. アルゴリズムをカプセル化してファクトリクラスに¶
Model
の移植性とWorld
と
Simulator
の一貫したAPIのおかげで、アルゴリズムに共通するスクリプトを書くのは非常に簡単です。
しかし、アルゴリズムを切り替えるときには、コード内のクラスの名前を1つずつ書き直す必要があります。
この問題を回避するために、E-Cell4はアルゴリズムごとに
Factory
クラスも提供しています。 Factory
はWorld
と
Simulator
を構造化に必要な引数でカプセル化します。
Factory
クラスを使うことで、あなたのスクリプトはアルゴリズムの変更に対してポータブルかつロバストになります。
In [22]:
from ecell4.gillespie import GillespieFactory
from ecell4.ode import ODEFactory
from ecell4.spatiocyte import SpatiocyteFactory
from ecell4.bd import BDFactory
from ecell4.meso import MesoscopicFactory
from ecell4.egfrd import EGFRDFactory
Factory
はcreate_world
と
create_simulator
の2つの関数を提供します。
In [23]:
def singlerun(f, m):
w = f.create_world(Real3(1, 1, 1))
w.bind_to(m)
w.add_molecules(Species('C'), 60)
sim = f.create_simulator(w)
sim.run(0.01)
print(sim.t(), w.num_molecules(Species('C')))
上記のsinglerun
はアルゴリズムから独立しています。
したがって、単にFactory
を切り替えるだけで、結果を簡単に比較することができます。
In [24]:
singlerun(GillespieFactory(), m)
singlerun(ODEFactory(), m)
singlerun(SpatiocyteFactory(), m)
singlerun(BDFactory(bd_dt_factor=1), m)
singlerun(MesoscopicFactory(), m)
singlerun(EGFRDFactory(), m)
(0.01, 60)
(0.01, 59)
(0.01, 60)
(0.01, 60)
(0.01, 60)
(0.01, 60)
World
やSimulator
を初期化するためにいくつかのパラメータを指定する必要があるとき、
run_simulation
はsolver
の代わりに
Factory
を受け取ります。
In [25]:
from ecell4.util import run_simulation
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=GillespieFactory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=ODEFactory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=SpatiocyteFactory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=BDFactory(bd_dt_factor=1))[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=MesoscopicFactory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=EGFRDFactory())[-1])
[0.01, 0.0, 0.0, 60.0]
[0.01, 0.17972919304001073, 0.17972919304001067, 59.82027080696036]
[0.01, 0.0, 0.0, 60.0]
[0.01, 0.0, 0.0, 60.0]
[0.01, 0.0, 0.0, 60.0]
[0.01, 0.0, 0.0, 60.0]