In [1]:
%matplotlib inline
from ecell4 import *
from ecell4_base.core import *
In [2]:
sp1 = Species("A")
print(sp1.serial())
There are some naming conventions for the name of Species
.
This naming convention requires careful use of special symbols (e.g. parenthesis ()
, dot .
, underbar _
), numbers and blank.
Species
has a set of APIs for handling its attributes:
In [3]:
sp1.set_attribute("radius", 0.005)
sp1.set_attribute("D", 1)
sp1.set_attribute("location", "cytoplasm")
print(sp1.has_attribute("radius"))
print(sp1.get_attribute("radius").magnitude)
print(sp1.get_attribute("radius").units)
print(sp1.has_attribute("location"))
print(sp1.get_attribute("location"))
sp1.remove_attribute("radius")
print(sp1.has_attribute("radius"))
The arguments in set_attribute
is the name of attribute and its value. The name of an attribute is given as a string, and its value is a string, integer, float, or boolean. get_attribute
returns the value set. For an integer or float attribute, get_attribute
returns a quantity object, which is a pair of value (magnitude
) and unit (units
).
There is a shortcut to set the attributes above at once because radius
, D
(a diffusion coefficient) and location
are frequently used.
In [4]:
sp1 = Species("A", 0.005, 1, "cytoplasm") # serial, radius, D, location
The equality between Species
is just evaluated based on their serial
:
In [5]:
print(Species("A") == Species("B"), Species("A") == Species("A"))
A Species consists of one or more UnitSpecies:
In [6]:
sp1 = Species()
usp1 = UnitSpecies("C")
print(usp1.serial())
sp1.add_unit(usp1)
sp1.add_unit(UnitSpecies("A"))
sp1.add_unit(UnitSpecies("B"))
print(sp1.serial(), len(sp1.units()))
A Species
can be reproduced from a serial. In a serial
, all UnitSpecies
are joined with the separator, dot .
. The order of UnitSpecies
affects the Species
comparison.
In [7]:
sp1 = Species("C.A.B")
print(sp1.serial())
print(Species("A.B.C") == Species("C.A.B"))
print(Species("A.B.C") == Species("A.B.C"))
UnitSpecies
can have sites. Sites consists of a name
, state
and bond
, and are sorted automatically in UnitSpecies
. name
must be unique in a UnitSpecies
. All the value have to be string. Do not include parenthesis, dot and blank, and not start from numbers except for bond
.
In [8]:
usp1 = UnitSpecies("A")
usp1.add_site("us", "u", "")
usp1.add_site("ps", "p", "_")
usp1.add_site("bs", "", "_")
print(usp1.serial())
UnitSpecies
can be also reproduced from its serial. Please be careful with the order of sites where a site with a state must be placed after sites with no state specification:
In [9]:
usp1 = UnitSpecies()
usp1.deserialize("A(bs^_, us=u, ps=p^_)")
print(usp1.serial())
Of course, a site of UnitSpecies
is available even in Species
' serial.
In [10]:
sp1 = Species("A(bs^1, ps=u).A(bs, ps=p^1)")
print(sp1.serial())
print(len(sp1.units()))
The information (UnitSpecies
and its site
) is used for rule-based modeling. The way of rule-based modeling in E-Cell4 will be explained in 7. Introduction of Rule-based Modeling local ipynb readthedocs.
In [11]:
rr1 = ReactionRule()
rr1.add_reactant(Species("A"))
rr1.add_reactant(Species("B"))
rr1.add_product(Species("C"))
rr1.set_k(1.0)
Here is a binding reaction from A
and B
to C
. In this reaction definition, you don't need to set attributes to Species
. The above series of operations can be written in one line using create_binding_reaction_rule(Species("A"), Species("B"), Species("C"), 1.0)
.
You can use as_string
function to check ReactionRule
:
In [12]:
rr1 = create_binding_reaction_rule(Species("A"), Species("B"), Species("C"), 1.0)
print(rr1.as_string())
You can also provide components to the constructor:
In [13]:
rr1 = ReactionRule([Species("A"), Species("B")], [Species("C")], 1.0)
print(rr1.as_string())
Basically ReactionRule
represents a mass action reaction with the rate k
. ode
solver also supports rate laws thought it's under development yet. ode.ODERatelaw
is explained in 6. How to Solve ODEs with Rate Law Functions local ipynb readthedocs.
In [14]:
sp1 = Species("A", 0.005, 1)
sp2 = Species("B", 0.005, 1)
sp3 = Species("C", 0.01, 0.5)
In [15]:
rr1 = create_binding_reaction_rule(Species("A"), Species("B"), Species("C"), 0.01)
rr2 = create_unbinding_reaction_rule(Species("C"), Species("A"), Species("B"), 0.3)
You can put the Species
and ReactionRule
with add_species_attribute
and add_reaction_rule
.
In [16]:
m1 = NetworkModel()
m1.add_species_attribute(sp1)
m1.add_species_attribute(sp2)
m1.add_species_attribute(sp3)
m1.add_reaction_rule(rr1)
m1.add_reaction_rule(rr2)
Now we have a simple model with the binding and unbinding reactions. You can use species_attributes
and reaction_rules
to check the Model
.
In [17]:
print([sp.serial() for sp in m1.species_attributes()])
print([rr.as_string() for rr in m1.reaction_rules()])
The Species
attributes are required for the spatial Model
, but not required for the nonspatial Model
(i.e. gillespie
or ode
). The attribute pushed first has higher priority than one pushed later. You can also attribute a Species
based on the attributes in a Model
.
In [18]:
sp1 = Species("A")
print(sp1.has_attribute("radius"))
sp2 = m1.apply_species_attributes(sp1)
print(sp2.has_attribute("radius"))
print(sp2.get_attribute("radius").magnitude)
For your information, all functions related to Species
, ReactionRule
and NetworkModel
above are also available on C++ in the same way.
Finally, you can solve this model with run_simulation
as explained in 1. Brief Tour of E-Cell4 Simulations local ipynb readthedocs :
In [19]:
run_simulation(10.0, model=m1, y0={'C': 60})
As shown in 1. Brief Tour of E-Cell4 Simulations local ipynb readthedocs, E-Cell4 also provides the easier way to build a model using with
statement:
In [20]:
with species_attributes():
A | B | {'radius': 0.005, 'D': 1}
C | {'radius': 0.01, 'D': 0.5}
with reaction_rules():
A + B == C | (0.01, 0.3)
m1 = get_model()
For reversible reactions, <>
is also available instead of ==
on Python 2, but deprecated on Python3. In the with
statement, undeclared variables are automaticaly assumed to be a Species
. Any Python variables, functions and statement are available even in the with
block.
In [21]:
from math import log
ka, kd, kf = 0.01, 0.3, 0.1
tau = 10.0
with reaction_rules():
E0 + S == ES | (ka, kd)
if tau > 0:
ES > E1 + P | kf
E1 > E0 | log(2) / tau
else:
ES > E0 + P | kf
m1 = get_model()
del ka, kd, kf, tau
Meanwhile, once some variable is declared even outside the block, you can not use its name as a Species
as follows:
In [22]:
A = 10
try:
with reaction_rules():
A + B == C | (0.01, 0.3)
except Exception as e:
print(repr(e))
del A
where A + B == C
exactly means 10 + B == C
.
In the absence of left or right hand side of ReactionRule
like synthesis or degradation, you may want to describe like:
with reaction_rules():
A > | 1.0 # XXX: will raise SyntaxError
> A | 1.0 # XXX: will raise SyntaxError
However, this is not accepted because of SyntaxError
on Python. To describe this case, a special operator, tilde ~
, is available. ~
sets a stoichiometric coefficient of the following Species
as zero, which means the Species
is just ignored in the ReactionRule
.
In [23]:
with reaction_rules():
~A > A | 2.0 # equivalent to `create_synthesis_reaction_rule`
A > ~A | 1.0 # equivalent to `create_degradation_reaction_rule`
m1 = get_model()
print([rr.as_string() for rr in m1.reaction_rules()])
The following Species
' name is not necessarily needed to be the same as others. The model above describes $[A]'=2-[A]$:
In [24]:
from math import exp
run_simulation(10.0, model=m1, opt_args=['-', lambda t: 2.0 * (1 - exp(-t)), '--'])
A chain of reactions can be described in one line. To split a line into two or more physical lines, wrap lines in a parenthesis:
In [25]:
with reaction_rules():
(E + S == ES | (0.5, 1.0)
> E + P | 1.5)
m1 = get_model()
print([rr.as_string() for rr in m1.reaction_rules()])
The method uses global variables in ecell4.util.decorator
(e.g. REACTION_RULES
) to cache objects created in the with
statement:
In [26]:
import ecell4.util.decorator
with reaction_rules():
A + B == C | (0.01, 0.3)
print(ecell4.util.decorator.REACTION_RULES) #XXX: Only for debugging
get_model()
print(ecell4.util.decorator.REACTION_RULES) #XXX: Only for debugging
Python decorator functions are also available. Decorator functions improve the modularity of the Model
.
In [27]:
@species_attributes
def attrgen1(radius, D):
A | B | {'radius': radius, 'D': D}
C | {'radius': radius * 2, 'D': D * 0.5}
@reaction_rules
def rrgen1(kon, koff):
A + B == C | (kon, koff)
attrs1 = attrgen1(0.005, 1)
rrs1 = rrgen1(0.01, 0.3)
print(attrs1)
print(rrs1)
In contrast to the with
statement, do not add parentheses after the decorator here. The functions decorated by reaction_rules
and species_attributes
return a list of ReactionRule
s and Species
respectively. The list can be registered to Model
at once by using add_reaction_rules
and add_species_attributes
.
In [28]:
m1 = NetworkModel()
m1.add_species_attributes(attrs1)
m1.add_reaction_rules(rrs1)
print(m1.num_reaction_rules())
This method is modular and reusable relative to the way using with
statement.