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]