Here, we explain the basics of World
classes. In E-Cell4, six types of World classes are supported now: spatiocyte.SpatiocyteWorld
, egfrd.EGFRDWorld
, bd.BDWorld
, meso.MesoscopicWorld
, gillespie.GillespieWorld
, and ode.ODEWorld
.
In the most of softwares, the initial condition is supposed to be a part of Model
. But, in E-Cell4, the initial condition must be set up as World
separately from Model
. World
stores an information about the state at a time point, such as a current time, the number of molecules, coordinate of molecules, structures, and random number generator. Meanwhile, Model
contains the type of interactions between molecules and the common properties of molecules.
In [1]:
from ecell4_base.core import *
In [2]:
from ecell4_base import *
World
classes accept different sets of arguments in the constructor, which determine the parameters specific to the algorithm. However, at least, all World
classes can be instantiated only with their size, named edge_lengths
. The type of edge_lengths
is Real3
, which represents a triplet of Real
s. In E-Cell4, all 3-dimensional positions are treated as a Real3
. See also 8. More about 1. Brief Tour of E-Cell4 Simulations.
In [3]:
edge_lengths = Real3(1, 2, 3)
w1 = gillespie.World(edge_lengths)
w2 = ode.World(edge_lengths)
w3 = spatiocyte.World(edge_lengths)
w4 = bd.World(edge_lengths)
w5 = meso.World(edge_lengths)
w6 = egfrd.World(edge_lengths)
World
has getter methods for the size and volume.
In [4]:
print(tuple(w1.edge_lengths()), w1.volume())
print(tuple(w2.edge_lengths()), w2.volume())
print(tuple(w3.edge_lengths()), w3.volume())
print(tuple(w4.edge_lengths()), w4.volume())
print(tuple(w5.edge_lengths()), w5.volume())
print(tuple(w6.edge_lengths()), w6.volume())
spatiocyte.World
(w3
) would have a bit larger volume to fit regular hexagonal close-packed (HCP) lattice.
Next, let's add molecules into the World. Here, you must give Species
attributed with radius and D to tell the shape of molecules. In the example below 0.0025 corresponds to radius and 1 to D. Positions of the molecules are randomly determined in the World
if needed. 10 in add_molecules function represents the number of molecules to be added.
In [5]:
sp1 = Species("A", 0.0025, 1)
w1.add_molecules(sp1, 10)
w2.add_molecules(sp1, 10)
w3.add_molecules(sp1, 10)
w4.add_molecules(sp1, 10)
w5.add_molecules(sp1, 10)
w6.add_molecules(sp1, 10)
After model is bound to world, you do not need to rewrite the radius and D once set in Species (unless you want to change it).
In [6]:
m = NetworkModel()
m.add_species_attribute(Species("A", 0.0025, 1))
m.add_species_attribute(Species("B", 0.0025, 1))
w1.bind_to(m)
w2.bind_to(m)
w3.bind_to(m)
w4.bind_to(m)
w5.bind_to(m)
w6.bind_to(m)
w1.add_molecules(Species("B"), 20)
w2.add_molecules(Species("B"), 20)
w3.add_molecules(Species("B"), 20)
w4.add_molecules(Species("B"), 20)
w5.add_molecules(Species("B"), 20)
w6.add_molecules(Species("B"), 20)
Similarly, remove_molecules
and num_molecules_exact
are also available.
In [7]:
w1.remove_molecules(Species("B"), 5)
w2.remove_molecules(Species("B"), 5)
w3.remove_molecules(Species("B"), 5)
w4.remove_molecules(Species("B"), 5)
w5.remove_molecules(Species("B"), 5)
w6.remove_molecules(Species("B"), 5)
In [8]:
print(w1.num_molecules_exact(Species("A")), w2.num_molecules_exact(Species("A")), w3.num_molecules_exact(Species("A")), w4.num_molecules_exact(Species("A")), w5.num_molecules_exact(Species("A")), w6.num_molecules_exact(Species("A")))
print(w1.num_molecules_exact(Species("B")), w2.num_molecules_exact(Species("B")), w3.num_molecules_exact(Species("B")), w4.num_molecules_exact(Species("B")), w5.num_molecules_exact(Species("B")), w6.num_molecules_exact(Species("B")))
Unlike num_molecules_exact
, num_molecules
returns the numbers that match a given Species
in rule-based fashion. When all Species
in the World
has no molecular interaction, num_molecules
is same with num_molecules_exact
.
In [9]:
print(w1.num_molecules(Species("A")), w2.num_molecules(Species("A")), w3.num_molecules(Species("A")), w4.num_molecules(Species("A")), w5.num_molecules(Species("A")), w6.num_molecules(Species("A")))
print(w1.num_molecules(Species("B")), w2.num_molecules(Species("B")), w3.num_molecules(Species("B")), w4.num_molecules(Species("B")), w5.num_molecules(Species("B")), w6.num_molecules(Species("B")))
World
holds its simulation time.
In [10]:
print(w1.t(), w2.t(), w3.t(), w4.t(), w5.t(), w6.t())
w1.set_t(1.0)
w2.set_t(1.0)
w3.set_t(1.0)
w4.set_t(1.0)
w5.set_t(1.0)
w6.set_t(1.0)
print(w1.t(), w2.t(), w3.t(), w4.t(), w5.t(), w6.t())
Finally, you can save
and load
the state of a World
into/from a HDF5 file.
In [11]:
w1.save("gillespie.h5")
w2.save("ode.h5")
w3.save("spatiocyte.h5")
w4.save("bd.h5")
w5.save("meso.h5")
w6.save("egfrd.h5")
del w1, w2, w3, w4, w5, w6
In [12]:
w1 = gillespie.World()
w2 = ode.World()
w3 = spatiocyte.World()
w4 = bd.World()
w5 = meso.World()
w6 = egfrd.World()
print(w1.t(), tuple(w1.edge_lengths()), w1.volume(), w1.num_molecules(Species("A")), w1.num_molecules(Species("B")))
print(w2.t(), tuple(w2.edge_lengths()), w2.volume(), w2.num_molecules(Species("A")), w2.num_molecules(Species("B")))
print(w3.t(), tuple(w3.edge_lengths()), w3.volume(), w3.num_molecules(Species("A")), w3.num_molecules(Species("B")))
print(w4.t(), tuple(w4.edge_lengths()), w4.volume(), w4.num_molecules(Species("A")), w4.num_molecules(Species("B")))
print(w5.t(), tuple(w5.edge_lengths()), w5.volume(), w5.num_molecules(Species("A")), w5.num_molecules(Species("B")))
print(w6.t(), tuple(w6.edge_lengths()), w6.volume(), w6.num_molecules(Species("A")), w6.num_molecules(Species("B")))
In [13]:
w1.load("gillespie.h5")
w2.load("ode.h5")
w3.load("spatiocyte.h5")
w4.load("bd.h5")
w5.load("meso.h5")
w6.load("egfrd.h5")
print(w1.t(), tuple(w1.edge_lengths()), w1.volume(), w1.num_molecules(Species("A")), w1.num_molecules(Species("B")))
print(w2.t(), tuple(w2.edge_lengths()), w2.volume(), w2.num_molecules(Species("A")), w2.num_molecules(Species("B")))
print(w3.t(), tuple(w3.edge_lengths()), w3.volume(), w3.num_molecules(Species("A")), w3.num_molecules(Species("B")))
print(w4.t(), tuple(w4.edge_lengths()), w4.volume(), w4.num_molecules(Species("A")), w4.num_molecules(Species("B")))
print(w5.t(), tuple(w5.edge_lengths()), w5.volume(), w5.num_molecules(Species("A")), w5.num_molecules(Species("B")))
print(w6.t(), tuple(w6.edge_lengths()), w6.volume(), w6.num_molecules(Species("A")), w6.num_molecules(Species("B")))
del w1, w2, w3, w4, w5, w6
All the World
classes also accept a HDF5 file path as an unique argument of the constructor.
In [14]:
print(gillespie.World("gillespie.h5").t())
print(ode.World("ode.h5").t())
print(spatiocyte.World("spatiocyte.h5").t())
print(bd.World("bd.h5").t())
print(meso.World("meso.h5").t())
print(egfrd.World("egfrd.h5").t())
In [15]:
w1 = gillespie.World()
w2 = ode.World()
w3 = spatiocyte.World()
w4 = bd.World()
w5 = meso.World()
w6 = egfrd.World()
First, you can place a molecule at the certain position with new_particle
.
In [16]:
sp1 = Species("A", 0.0025, 1)
pos = Real3(0.5, 0.5, 0.5)
(pid1, p1), suc1 = w1.new_particle(sp1, pos)
(pid2, p2), suc2 = w2.new_particle(sp1, pos)
pid3 = w3.new_particle(sp1, pos)
(pid4, p4), suc4 = w4.new_particle(sp1, pos)
(pid5, p5), suc5 = w5.new_particle(sp1, pos)
(pid6, p6), suc6 = w6.new_particle(sp1, pos)
new_particle
returns a particle created and whether it's succeeded or not. The resolution in representation of molecules differs. For example, GillespieWorld
has almost no information about the coordinate of molecules. Thus, it simply ignores the given position, and just counts up the number of molecules here.
ParticleID
is a pair of Integer
s named lot
and serial
.
In [17]:
print(pid6.lot(), pid6.serial())
print(pid6 == ParticleID((0, 1)))
Particle simulators, i.e. spatiocyte
, bd
and egfrd
, provide an interface to access a particle by its id. has_particle
returns if a particles exists or not for the given ParticleID
.
In [18]:
# print(w1.has_particle(pid1))
# print(w2.has_particle(pid2))
print(w3.has_particle(pid3)) # => True
print(w4.has_particle(pid4)) # => True
# print(w5.has_particle(pid5))
print(w6.has_particle(pid6)) # => True
After checking the existency, you can get the partcle by get_particle
as follows.
In [19]:
# pid1, p1 = w1.get_particle(pid1)
# pid2, p2 = w2.get_particle(pid2)
pid3, p3 = w3.get_particle(pid3)
pid4, p4 = w4.get_particle(pid4)
# pid5, p5 = w5.get_particle(pid5)
pid6, p6 = w6.get_particle(pid6)
Particle
consists of species
, position
, radius
and D
.
In [20]:
# print(p1.species().serial(), tuple(p1.position()), p1.radius(), p1.D())
# print(p2.species().serial(), tuple(p2.position()), p2.radius(), p2.D())
print(p3.species().serial(), tuple(p3.position()), p3.radius(), p3.D())
print(p4.species().serial(), tuple(p4.position()), p4.radius(), p4.D())
# print(p5.species().serial(), tuple(p5.position()), p5.radius(), p5.D())
print(p6.species().serial(), tuple(p6.position()), p6.radius(), p6.D())
In the case of spatiocyte
, a particle position is automatically round to the center of the voxel nearest to the given position.
You can even move the position of the particle. update_particle
replace the particle specified with the given ParticleID
with the given Particle
and return False
. If no corresponding particle is found, create new particle and return True
. If you give a Particle
with the different type of Species
, the Species
of the Particle
will be also changed.
In [21]:
newp = Particle(sp1, Real3(0.3, 0.3, 0.3), 0.0025, 1)
# print(w1.update_particle(pid1, newp))
# print(w2.update_particle(pid2, newp))
print(w3.update_particle(pid3, newp))
print(w4.update_particle(pid4, newp))
# print(w5.update_particle(pid5, newp))
print(w6.update_particle(pid6, newp))
list_particles
and list_particles_exact
return a list of pairs of ParticleID
and Particle
in the World
. World
automatically makes up for the gap with random numbers. For example, GillespieWorld
returns a list of positions randomly distributed in the World
size.
In [22]:
print(w1.list_particles_exact(sp1))
# print(w2.list_particles_exact(sp1)) # ODEWorld has no member named list_particles
print(w3.list_particles_exact(sp1))
print(w4.list_particles_exact(sp1))
print(w5.list_particles_exact(sp1))
print(w6.list_particles_exact(sp1))
You can remove a specific particle with remove_particle
.
In [23]:
# w1.remove_particle(pid1)
# w2.remove_particle(pid2)
w3.remove_particle(pid3)
w4.remove_particle(pid4)
# w5.remove_particle(pid5)
w6.remove_particle(pid6)
# print(w1.has_particle(pid1))
# print(w2.has_particle(pid2))
print(w3.has_particle(pid3)) # => False
print(w4.has_particle(pid4)) # => False
# print(w5.has_particle(pid5))
print(w6.has_particle(pid6)) # => False
In [24]:
w = spatiocyte.World(Real3(1, 2, 3), voxel_radius=0.01)
w.bind_to(m)
The size of a single lattice, called Voxel
, can be obtained by voxel_radius()
. SpatiocyteWorld
has methods to get the numbers of rows, columns, and layers. These sizes are automatically calculated based on the given edge_lengths
at the construction.
In [25]:
print(w.voxel_radius()) # => 0.01
print(tuple(w.shape())) # => (64, 152, 118)
# print(w.col_size(), w.row_size(), w.layer_size()) # => (64, 152, 118)
print(w.size()) # => 1147904 = 64 * 152 * 118
A position in the lattice-based space is treated as an Integer3
, column, row and layer, called a global coordinate. Thus, SpatiocyteWorld
provides the function to convert the Real3
into a lattice-based coordinate.
In [26]:
# p1 = Real3(0.5, 0.5, 0.5)
# g1 = w.position2global(p1)
# p2 = w.global2position(g1)
# print(tuple(g1)) # => (31, 25, 29)
# print(tuple(p2)) # => (0.5062278801751902, 0.5080682368868706, 0.5)
In SpatiocyteWorld
, the global coordinate is translated to a single integer. It is just called a coordinate. You can also treat the coordinate as in the same way with a global coordinate.
In [27]:
# p1 = Real3(0.5, 0.5, 0.5)
# c1 = w.position2coordinate(p1)
# p2 = w.coordinate2position(c1)
# g1 = w.coord2global(c1)
# print(c1) # => 278033
# print(tuple(p2)) # => (0.5062278801751902, 0.5080682368868706, 0.5)
# print(tuple(g1)) # => (31, 25, 29)
With these coordinates, you can handle a Voxel
, which represents a Particle
object. Instead of new_particle
, new_voxel
provides the way to create a new Voxel
with a coordinate.
In [28]:
# c1 = w.position2coordinate(Real3(0.5, 0.5, 0.5))
# ((pid, v), is_succeeded) = w.new_voxel(Species("A"), c1)
# print(pid, v, is_succeeded)
A Voxel
consists of species
, coordinate
, radius
and D
.
In [29]:
# print(v.species().serial(), v.coordinate(), v.radius(), v.D()) # => (u'A', 278033, 0.0025, 1.0)
Of course, you can get a voxel and list voxels with get_voxel
and list_voxels_exact
similar to get_particle
and list_particles_exact
.
In [30]:
# print(w.num_voxels_exact(Species("A")))
# print(w.list_voxels_exact(Species("A")))
# print(w.get_voxel(pid))
You can move and update the voxel with update_voxel
corresponding to update_particle
.
In [31]:
# c2 = w.position2coordinate(Real3(0.5, 0.5, 1.0))
# w.update_voxel(pid, Voxel(v.species(), c2, v.radius(), v.D()))
# pid, newv = w.get_voxel(pid)
# print(c2) # => 278058
# print(newv.species().serial(), newv.coordinate(), newv.radius(), newv.D()) # => (u'A', 278058, 0.0025, 1.0)
# print(w.num_voxels_exact(Species("A"))) # => 1
Finally, remove_voxel
remove a voxel as remove_particle
does.
In [32]:
# print(w.has_voxel(pid)) # => True
# w.remove_voxel(pid)
# print(w.has_voxel(pid)) # => False
In [33]:
w1 = gillespie.World()
w2 = ode.World()
w3 = spatiocyte.World()
w4 = bd.World()
w5 = meso.World()
w6 = egfrd.World()
By using a Shape
object, you can confine initial positions of molecules to a part of World
. In the case below, 60 molecules are positioned inside the given Sphere
. Diffusion of the molecules placed here is NOT restricted in the Shape
. This Shape
is only for the initialization.
In [34]:
sp1 = Species("A", 0.0025, 1)
sphere = Sphere(Real3(0.5, 0.5, 0.5), 0.3)
w1.add_molecules(sp1, 60, sphere)
w2.add_molecules(sp1, 60, sphere)
w3.add_molecules(sp1, 60, sphere)
w4.add_molecules(sp1, 60, sphere)
w5.add_molecules(sp1, 60, sphere)
w6.add_molecules(sp1, 60, sphere)
A property of Species
, 'location'
, is available to restrict diffusion of molecules. 'location'
is not fully supported yet, but only supported in spatiocyte
and meso
. add_structure
defines a new structure given as a pair of Species
and Shape
.
In [35]:
membrane = SphericalSurface(Real3(0.5, 0.5, 0.5), 0.4) # This is equivalent to call `Sphere(Real3(0.5, 0.5, 0.5), 0.4).surface()`
w3.add_structure(Species("M"), membrane)
w5.add_structure(Species("M"), membrane)
After defining a structure, you can bind molecules to the structure as follows:
In [36]:
sp2 = Species("B", 0.0025, 0.1, "M") # `'location'` is the fourth argument
w3.add_molecules(sp2, 60)
w5.add_molecules(sp2, 60)
The molecules bound to a Species
named B
diffuse on a structure named M
, which has a shape of SphericalSurface
(a hollow sphere). In spatiocyte
, a structure is represented as a set of particles with Species
M
occupying a voxel. It means that molecules not belonging to the structure is not able to overlap the voxel and it causes a collision. On the other hand, in meso
, a structure means a list of subvolumes. Thus, a structure doesn't avoid an incursion of other particles.
In [37]:
rng1 = GSLRandomNumberGenerator()
print([rng1.uniform_int(1, 6) for _ in range(20)])
With no argument, the random number generator is always initialized with a seed, 0
.
In [38]:
rng2 = GSLRandomNumberGenerator()
print([rng2.uniform_int(1, 6) for _ in range(20)]) # => same as above
You can initialize the seed with an integer as follows:
In [39]:
rng2 = GSLRandomNumberGenerator()
rng2.seed(15)
print([rng2.uniform_int(1, 6) for _ in range(20)])
When you call the seed
function with no input, the seed is drawn from the current time.
In [40]:
rng2 = GSLRandomNumberGenerator()
rng2.seed()
print([rng2.uniform_int(1, 6) for _ in range(20)])
GSLRandomNumberGenerator
provides several ways to get a random number.
In [41]:
print(rng1.uniform(0.0, 1.0))
print(rng1.uniform_int(0, 100))
print(rng1.gaussian(1.0))
World
accepts a random number generator at the construction. As a default, GSLRandomNumberGenerator()
is used. Thus, when you don't give a generator, behavior of the simulation is always same (determinisitc).
In [42]:
rng = GSLRandomNumberGenerator()
rng.seed()
w1 = gillespie.World(Real3(1, 1, 1), rng=rng)
You can access the GSLRandomNumberGenerator
in a World
through rng
function.
In [43]:
print(w1.rng().uniform(0.0, 1.0))
rng()
returns a shared pointer to the GSLRandomNumberGenerator
. Thus, in the example above, rng
and w1.rng()
point exactly the same thing.