We showed an example of E-Cell4 spatial representation.
Next let's simulate the models with more detailed spatial representation called single molecule resolution.
In spatical Gillespie method, we divided the Space into smaller Space, then we diffuse the molecules in the Subvolume. However we simulated the molecules in the Subvolume as the number of the molecules, and the location of the molecules are NOT determinated.
In other words, the spatical resolution of spatical Gillespie method is equal to the side of the Subvolume $l$. To improve this resolution we need to make the size of $l$ small, but in this method the $l$ must be larger than the (at least) 10 times the diameter of molecule $R$.
How can we improve the spatical resolution to the size of the molecule? The answer is the simulation with single-molecule resolution. This method simulate the molecule not with the number of the molecules, but with the spatical reaction diffusion of each molecule.
E-Cell4 has multiple single-molecule resolution method, here we explain about Spatiocyte lattice-based method.
Spatiocyte treat each molecule as hard sphere and diffuse the molecules in hexagonal close-packed lattice.
Spatiocyte has ID for each molecule and the position of the molecule with single-molecule resolution. To use the time scale, Spatiocyte has 100 times smaller time-step than spatical Gillespie, because the time scale of diffusion increases with the square of the distance.
Next, let's use Spatiocyte.
In [1]:
from ecell4 import *
with species_attributes():
A | B | C | {'D': '1'}
with reaction_rules():
A + B == C | (0.01, 0.3)
m = get_model()
w = lattice.LatticeWorld(Real3(1, 1, 1), 0.005) # The second argument is 'voxel_radius'.
w.bind_to(m)
w.add_molecules(Species('C'), 60)
sim = lattice.LatticeSimulator(w)
obs = FixedIntervalNumberObserver(0.1, ('A', 'B', 'C'))
sim.run(10, obs)
There is a distinct difference in the second argument for LatticeWorld. This is called Voxel radius. Spatiocyte defines the locations of the molecules with dividing the Space with molecule size, and call the minimum unit for this Space as Voxel.
In most cases the size of the molecule would be good for Voxel radius.
In this example, we set 5 $\mathrm{nm}$ as the radius of the molecule in the Space with the side 1 $\mathrm{\mu m}$ .
It takes more time to simulate, but the result is same with ODE or Gillespie.
In [2]:
from ecell4 import *
with species_attributes():
A | {'D': '1'}
m = get_model()
w = lattice.LatticeWorld(Real3(1, 1, 1), 0.005)
w.bind_to(m)
(pid, p), suc = w.new_particle(Species('A'), Real3(0.5, 0.5, 0.5))
new_particle method places a particle to a coordinate in LatticeWorld, and returns the particle's pid, the information about the particle p, and verify whether the particle is cooked with suc. If a particle is already placed in the coordinate you can NOT place a particle over it and suc will be False and fail.
p contains the particle position, species type, radius, and diffusion coefficient. You can inspect the p with the particle's ID pid.
Let's check p.
In [3]:
pid, p = w.get_particle(pid)
print(p.species().serial()) # will print: A
print(p.radius(), p.D()) # will print: (0.005, 1.0)
print(tuple(p.position())) # will print: (0.49806291436591293, 0.49652123150307814, 0.5)
get_particle method receives the particle ID and returns the ID and particle (of cource the IDs are same). You can inspect the coordinate of the particle as Real3 with position method. It is hard to directly read the coordinate, here we printed it after converting to tuple. As you can see the tuple coodinate is slightly different from original Real3. This is because Spatiocyte can place the molecule only on the lattice. LatticeWorld places the molecule the nearest lattice for the argument Real3.
You can visualize the coordinate of the molecule with viz.plot_world method, and check the molecule in the center of the World.
In [4]:
viz.plot_world(w, save_image=True)
Out[4]:
And you can use Observer to track the trajectory of molecular diffusion process.
In [5]:
sim = lattice.LatticeSimulator(w)
obs = FixedIntervalTrajectoryObserver(0.002, (pid,))
sim.run(1, obs)
viz.plot_trajectory(obs)
Out[5]:
Here we visualized the trajectory with viz.plot_trajectory method, you can also obtain it as Real3 list with data method.
In [6]:
print(len(obs.data()))
print(len(obs.data()[0]))
data method returns nested list. First index means the index of the particle. Second index means the index of the Real3. In this case we threw just one particle, so the first result is 1, and next 501 means time-series coordinate of the only one particle (initial coordinate and the coordinates in 1/0.002 = 500 time points).
Also you can obtain the particles in bulk with list_particles method and species type.
In [7]:
w.add_molecules(Species('A'), 5)
particles = w.list_particles(Species('A'))
for pid, p in particles:
print(p.species().serial(), tuple(p.position()))
Please remember list_particles method, this method can be used for other World as well as add_molecules method.
On a different note, in Spatiocyte proper method to inspect the single molecule is list_voxels, and the coordinate is described with index of voxel (not Real3).
The models we have addressed are called second-order reaction. Let's look at the relationship between this second-order reaction and the diffusion coefficient in Spatiocyte.
In [8]:
from ecell4 import *
with species_attributes():
A | B | C | {'D': '1'}
with reaction_rules():
A + B > C | 1.0
m = get_model()
w = lattice.LatticeWorld(Real3(2, 1, 1), 0.005)
w.bind_to(m)
w.add_molecules(Species('A'), 120)
w.add_molecules(Species('B'), 120)
obs = FixedIntervalNumberObserver(0.005, ('A', 'B', 'C'))
sim = lattice.LatticeSimulator(w)
sim.run(1.0, obs)
%matplotlib inline
odew = ode.ODEWorld(Real3(2, 1, 1))
odew.bind_to(m)
odew.add_molecules(Species('A'), 120)
odew.add_molecules(Species('B'), 120)
odeobs = FixedIntervalNumberObserver(0.005, ('A', 'B', 'C'))
odesim = ode.ODESimulator(odew)
odesim.run(1.0, odeobs)
viz.plot_number_observer(obs, "-", odeobs, "--")
Although we used faster kinetic constant than before, the result is same. But by contrast with ODE simulation, you can find the difference between them (solid line is Spatiocyte, dash line is ODE).
Is this fault of Spatiocyte? (No) Actually Spatiocyte reaction rate couldn't be faster, while ODE reaction rate can be faster infinitely.
This is caused by the difference between the definition of reaction rate constant in ODE solver and single molecule simulation method.
The former is called macroscopic or effective reaction rate constant, the latter is called microscopic or intrinsic reaction rate constant.
The macroscopic rate represents the reaction rate in mixed molecular state, meanwhile microscopic rate represents the reactivity in molecule collision. So in microscopic perspective, the first thing molecules need to react is collision.
In Spatiocyte however you make this microscopic rate faster, you can NOT make the reaction rate faster than diffusion rate. This is called diffusion-limited. This is similar to what the molecules coordinated disproportionately need time to react.
It is known that there is a relationship between this macroscopic rate constant $k_\mathrm{on}$ and microscopic rate constant $k_a$ in 3D space.
$ \frac{1}{k_\mathrm{on}}=\frac{1}{k_a}+\frac{1}{4\pi RD_\mathrm{tot}} $
Here, $R$ is the sum of two molecule's radius in collision, $D_\mathrm{tot}$ is the sum of diffusion coefficients. In the case of the above IPython Notebook cell, $k_D=4\pi RD_\mathrm{tot}$ is almost 0.25 and microscopic rate constant is 1.0. So the macroscopic rate constant is almost 0.2.
(However unless you specify the configuration for Spatiocyte, the second order reaction rate must be slower than $3\sqrt{2} RD$, and the dissociation constant $k_D$ is also $3\sqrt{2} RD$.
The single molecule simulation method can separate molecular diffusion and reaction in accurate manner contrary to ODE or Gillespie method supposed well mixed system (that is diffusion coefficient is infinite).
However if the microscopic rate constant $k_D$ is small enough, the macroscopic rate constant is almost equal to microscopic one (reaction late-limit).
Next we explain a way to create a structure like cell membrane. Although The structure feature in E-Cell4 is still in development, Spatiocyte supports the structure on some level. Let's look a sphere structure as an example.
To restrict the molecular diffusion inside of the sphere, first we create it.
In [9]:
from ecell4 import *
with species_attributes():
A | {'D': '1', 'location': 'C'}
m = get_model()
In [10]:
w = lattice.LatticeWorld(Real3(1, 1, 1), 0.005)
w.bind_to(m)
sph = Sphere(Real3(0.5, 0.5, 0.5), 0.45)
print(w.add_structure(Species('C'), sph)) # will print 539805
In [11]:
viz.plot_world(w, save_image=True)
Out[11]:
The Sphere class first argument is the center of the sphere, and second argument is the radius. Then we created and added a Species named C.
The structure in the Spatiocyte method is described by filling the space with the Voxel. In the example above, the Voxels in the sphere are occupied with Species named C.
You can see those distribution with viz.plot_world. (However, the number of the species is too large to visualize. So we plot only a part of it, but actually the sphere is fully occupied with the species.)
Next we create Species moving inside this sphere. To that end we give location attribute to the Species.
and add_molecules to the World.
In [12]:
w.add_molecules(Species('A'), 120)
viz.plot_world(w, species_list=('A',), save_image=True) # visualize A-molecules only
Out[12]:
Now we restricted the trajectories of Species A on the structure of Species C, and add_molecules works like that.
As a note, you need to create the structure before add_molecule.
We can use FixedIntervalTrajectoryObserver to check the restriction of the diffusion area.
In [13]:
pid_list = [pid for pid, p in w.list_particles(Species('A'))[: 10]]
obs = FixedIntervalTrajectoryObserver(1e-3, pid_list)
sim = lattice.LatticeSimulator(w)
sim.run(1, obs)
viz.plot_trajectory(obs, save_image=True)
Out[13]:
pid_list is a list for 10 IDs of 60 A species. The trajectories are colored by this 10 species. Certainly the trajectories are restricted in the sphere.
At the end we explain about molecular translocation among the structures.
A species without location attribute is not an member of any structures.
In the example above, if you do NOT write location attribute with Species A, A is placed outside of the sphere.
Next let's create a surface structure. To create a surface we need to use three Real3, those are original point (origin) and two axis vector (unit0, unit1).
ps = PlanarSurface(origin, unit0, unit1)
Use this ps and suppose Species A on the surface and a normal Species B.
In [1]:
from ecell4 import *
with species_attributes():
A | {'D': '0.1', 'location': 'M'}
B | {'D': '1'}
m = get_model()
w = lattice.LatticeWorld(Real3(1, 1, 1))
w.bind_to(m)
origin = Real3(0, 0, 0.5)
unit0 = Real3(1, 0, 0)
unit1 = Real3(0, 1, 0)
w.add_structure(
Species('M'), PlanarSurface(origin, unit0, unit1)) # Create a structure first
w.add_molecules(Species('B'), 480) # Throw-in B-molecules
viz.plot_world(w, species_list=('A', 'B'))
Out[1]:
It might be hard to understand, but actually the species B are placed only on a surface. Then how can we make absorbed this species B to a surface M and synthesize a species A.
with reaction_rules():
B + M == A | (1e-3, 1.5)
This means that a species B becomes A when B collides with a structure M.
And a species A dissociates and becomes M and B on in the reverse reaction direction.
Now you can simulate a model with structure.
In [2]:
%matplotlib inline
with reaction_rules():
B + M == A | (1e-3, 1.5)
sim = lattice.LatticeSimulator(w)
obs = NumberObserver(('A', 'B'))
sim.run(2, obs)
viz.plot_number_observer(obs)
viz.plot_world(w, species_list=('A', 'B'))
Out[2]:
In the dissociation from a structure, you can skip to write the structure. But in the binding, you can NOT. Because it is impossible to create A from B without any M around there.
By contrast the species A wherever on the sphere M can create the species B.
What is more, in the case of the diffusion the first order reaction occurs in the presence or absence of the structure.
But in the case of the binding the second order reaction turns into the first order reaction and the meaning of rate constant also changes.
with reaction_rules():
B + M > A | 1e-3
A > B | 3.0 # means the same as A > B + M
In [ ]: