3. How to Setup the Initial Condition

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.

from ecell4_base.core import *

3.1. Common APIs of World

Even though World describes the spatial representation specific to the corresponding algorithm, it has compatible APIs. In this section, we introduce the common interfaces of the six World classes.

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 Reals. In E-Cell4, all 3-dimensional positions are treated as a Real3. See also 8. More about 1. Brief Tour of E-Cell4 Simulations.

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.

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())

(1.0, 2.0, 3.0) 6.0
(1.0, 2.0, 3.0) 6.0
(1.0124557603503803, 2.0091789367798976, 3.0) 6.102614364352381
(1.0, 2.0, 3.0) 6.0
(1.0, 2.0, 3.0) 6.0
(1.0, 2.0, 3.0) 6.0

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.

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).

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.

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)

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")))

10 10 10 10 10 10
15 15 15 15 15 15

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.

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")))

10 10 10 10 10 10
15 15 15 15 15 15

World holds its simulation time.

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())

0.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 1.0

Finally, you can save and load the state of a World into/from a HDF5 file.

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

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")))

0.0 (1.0, 1.0, 1.0) 1.0 0 0
0.0 (1.0, 1.0, 1.0) 1.0 0 0
0.0 (1.0124557603503803, 1.0045894683899488, 1.0) 1.0171023940587303 0 0
0.0 (1.0, 1.0, 1.0) 1.0 0 0
0.0 (1.0, 1.0, 1.0) 1.0 0 0
0.0 (1.0, 1.0, 1.0) 1.0 0 0

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

1.0 (1.0, 2.0, 3.0) 6.0 10 15
1.0 (1.0, 2.0, 3.0) 6.0 10 15
1.0 (1.0124557603503803, 2.0091789367798976, 3.0) 6.102614364352381 10 15
1.0 (1.0, 2.0, 3.0) 6.0 10 15
1.0 (1.0, 2.0, 3.0) 6.0 10 15
1.0 (1.0, 2.0, 3.0) 6.0 10 15

All the World classes also accept a HDF5 file path as an unique argument of the constructor.

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())

1.0
1.0
1.0
1.0
1.0
1.0

3.2. How to Get Positions of Molecules

World also has the common functions to access the coordinates of the molecules.

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.

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 Integers named lot and serial.

print(pid6.lot(), pid6.serial())
print(pid6 == ParticleID((0, 1)))

0 1
True

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.

# 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

True
True
True

After checking the existency, you can get the partcle by get_particle as follows.

# 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.

# 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())

A (0.5062278801751902, 0.5080682368868706, 0.5) 0.0025 1.0
A (0.5, 0.5, 0.5) 0.0025 1.0
A (0.5, 0.5, 0.5) 0.0025 1.0

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.

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))

False
False
False

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.

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))

[(<ecell4_base.core.ParticleID object at 0x150e376a79a8>, <ecell4_base.core.Particle object at 0x150e376a79f0>)]
[(<ecell4_base.core.ParticleID object at 0x150e376a79a8>, <ecell4_base.core.Particle object at 0x150e376a79c0>)]
[(<ecell4_base.core.ParticleID object at 0x150e376a79a8>, <ecell4_base.core.Particle object at 0x150e376a7930>)]
[(<ecell4_base.core.ParticleID object at 0x150e376a79a8>, <ecell4_base.core.Particle object at 0x150e376a79f0>)]
[(<ecell4_base.core.ParticleID object at 0x150e376a79a8>, <ecell4_base.core.Particle object at 0x150e376a79c0>)]

You can remove a specific particle with remove_particle.

# 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

False
False
False

3.3. Lattice-based Coordinate

In addition to the common interface, each World can have their own interfaces. As an example, we explain methods to handle lattice-based coordinate here. SpatiocyteWorld is based on a space discretized to hexiagonal close packing lattices, LatticeSpace.

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.

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

0.01
(64, 152, 118)
1147904

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.

# 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.

# 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.

# 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.

# 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.

# 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.

# 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.

# print(w.has_voxel(pid))  # => True
# w.remove_voxel(pid)
# print(w.has_voxel(pid))  # => False

3.4 Structure

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.

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.

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:

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.

3.5. Random Number Generator

A random number generator is also a part of World. All World except ODEWorld store a random number generator, and updates it when the simulation needs a random value. On E-Cell4, only one class GSLRandomNumberGenerator is implemented as a random number generator.

rng1 = GSLRandomNumberGenerator()
print([rng1.uniform_int(1, 6) for _ in range(20)])

[6, 1, 2, 6, 2, 3, 6, 5, 4, 5, 5, 4, 2, 5, 4, 2, 3, 3, 2, 2]

With no argument, the random number generator is always initialized with a seed, 0.

rng2 = GSLRandomNumberGenerator()
print([rng2.uniform_int(1, 6) for _ in range(20)])  # => same as above

[6, 1, 2, 6, 2, 3, 6, 5, 4, 5, 5, 4, 2, 5, 4, 2, 3, 3, 2, 2]

You can initialize the seed with an integer as follows:

rng2 = GSLRandomNumberGenerator()
rng2.seed(15)
print([rng2.uniform_int(1, 6) for _ in range(20)])

[6, 5, 2, 4, 1, 1, 3, 5, 2, 6, 4, 1, 2, 5, 2, 5, 1, 2, 2, 6]

When you call the seed function with no input, the seed is drawn from the current time.

rng2 = GSLRandomNumberGenerator()
rng2.seed()
print([rng2.uniform_int(1, 6) for _ in range(20)])

[5, 3, 1, 4, 3, 2, 4, 3, 1, 5, 2, 2, 2, 5, 3, 4, 3, 5, 2, 1]

GSLRandomNumberGenerator provides several ways to get a random number.

print(rng1.uniform(0.0, 1.0))
print(rng1.uniform_int(0, 100))
print(rng1.gaussian(1.0))

0.03033520421013236
33
0.8935555455208181

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).

rng = GSLRandomNumberGenerator()
rng.seed()
w1 = gillespie.World(Real3(1, 1, 1), rng=rng)

You can access the GSLRandomNumberGenerator in a World through rng function.

print(w1.rng().uniform(0.0, 1.0))

0.6781027028337121

rng() returns a shared pointer to the GSLRandomNumberGenerator. Thus, in the example above, rng and w1.rng() point exactly the same thing.