Once you read through 1. Brief Tour of E-Cell4 Simulations, it is NOT difficult to use World and Simulator.
volume and {'C': 60} is equivalent of the World and solver is the Simulator below.
In [1]:
%matplotlib inline
from ecell4 import *
with reaction_rules():
A + B == C | (0.01, 0.3)
run_simulation(10.0, {'C': 60}, volume=1.0)
Here we give you a breakdown for run_simulation.
run_simulation use ODE simulator by default, so we create ode.World step by step.
In [2]:
from ecell4_base.core import *
from ecell4_base import *
In [3]:
w = ode.World(Real3(1, 1, 1))
Real3 is a coordinate vector.
In this example, the first argument for ode.World constructor is a cube.
Note that you can NOT use volume for ode.World argument, like run_simulation argument.
Now you created a cube box for simulation, next let's throw molecules into the cube.
In [4]:
w = ode.World(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
print(w.t(), w.num_molecules(Species('C'))) # must return (0.0, 60)
Use add_molecules to add molecules, remove_molecules to remove molecules, num_molecules to know the number of molecules.
First argument for each method is the Species you want to know.
You can get current time by t method.
However the number of molecules in ODE solver is real number, in these _molecules functions work only for integer number.
When you handle real numbers in ODE, use set_value and get_value.
In [5]:
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)
You can not print the contents in Real3 object directly.
You need to convert Real3 to Python tuple or list once.
In [6]:
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
You can use basic function like dot_product.
Of course, you can convert Real3 to numpy array too.
In [7]:
import numpy
a = numpy.asarray(tuple(Real3(1, 2, 3)))
print(a) # must print [ 1. 2. 3.]
Integer3 represents a triplet of integers.
In [8]:
g = Integer3(1, 2, 3)
print(tuple(g))
Of course, you can also apply simple arithmetics to Integer3.
In [9]:
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
In [10]:
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.Simulator(w, m)
sim.run(10.0)
then call run method, the simulation will run.
In this example the simulation runs for 10 seconds.
You can check the state of the World like this.
In [11]:
print(w.t(), w.num_molecules(Species('C'))) # must return (10.0, 30)
You can see that the number of the Species C decreases from 60 to 30.
World describes the state at a timepoint, so you can NOT tack the transition during the simulation with the World.
To obtain the time-series result, use Observer.
In [12]:
w = ode.World(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
sim = ode.Simulator(w, m)
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]]
There are several types of Observers for E-Cell4.
FixedIntervalNumberObserver is the simplest Observer to obtain the time-series result.
As its name suggests, this Observer records the number of molecules for each time-step.
The 1st argument is the time-step, the 2nd argument is the molecule types.
You can check the result with data method, but there is a shortcut for this.
In [13]:
show(obs)
This plots the time-series result easily.
We explained the internal of run_simulation function.
When you change the World after creating the Simulator, you need to indicate it to Simulator.
So do NOT forget to call sim.initialize() after that.
In [14]:
from ecell4 import *
with reaction_rules():
A + B == C | (0.01, 0.3)
m = get_model()
# ode.World -> gillespie.World
w = gillespie.World(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
# ode.Simulator -> gillespie.Simulator
sim = gillespie.Simulator(w, m)
obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)
show(obs)
World and Simulator never change the Model itself, so you can switch several Simulators for one Model.