In [2]:
%matplotlib inline
import numpy
from ecell4 import *
with reaction_rules():
A + B == C | (0.01, 0.3)
y = run_simulation(numpy.linspace(0, 10, 100), {'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.
You can create world like this.
In [3]:
w = ode.ODEWorld(Real3(1, 1, 1))
Real3 is a coordinate vector. In this example, the first argument for ODEWorld constructor is a cube. Note that you can NOT use volume for ode.ODEWorld argument, like run_simulation argument.
Now you created a cube box for simulation, next let's throw in molecules into it.
In [4]:
w = ode.ODEWorld(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
print(w.t(), w.num_molecules(Species('C'))) # will 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. If you use real number in ODE, use set_value and get_value.
Before the detail of Simulator, we explaing more about Real3.
In [5]:
pos = Real3(1, 2, 3)
print(pos) # will print <ecell4.core.Real3 object at 0x7f44e118b9c0>
print(tuple(pos)) # will print (1.0, 2.0, 3.0)
You can not print 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)) # will print 1.73205080757
print(dot_product(pos1, pos3)) # will print 9.0
You can use basic function like dot_product. Of course you can convert Real3 to numpy array too.
In [7]:
a = numpy.asarray(tuple(Real3(1, 2, 3)))
print(a) # will print [ 1. 2. 3.]
You can create a Simulator with Model and World like
In [8]:
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)
then call run method, the simulation will run. In this example the simulation runs for 10seconds.
You can check the state of the World like this.
In [9]:
print(w.t(), w.num_molecules(Species('C'))) # will return (10.0, 30)
You can see that the number of the Species C decreases from 60 to 30.
World describes the state of a timepoint, so you can NOT see the transition of the simulation with the World. To obtain the time-series result, use Observer.
In [10]:
obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)
print(obs.data()) # will return [[0.0, 0.0, 60.0], ..., [10.0, 29.994446899698026, 30.005553100301967]]
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 [11]:
viz.plot_number_observer(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
In [12]:
sim.initialize()
In [13]:
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 and Simulator do NOT change the Model, so you can switch several simulators for 1 model.