[PYTHON] Easy to use E-Cell 4 Intermediate

Preparation

Please refer to the separate article for installation.

-Try using E-Cell 4 on Windows 7 or Mac OS X

This article is an intermediate edition. It is assumed that you understand the beginner's edition.

-Easy to use E-Cell 4 for beginners

In addition, this article assumes the use of IPython Notebook.

How to use

Space in E-Cell 4

In the beginner's edition, we explained the three elements Model, World, and Simulator in E-Cell 4. We also tried simple calculations using the ordinary differential equation solver ode and the stochastic method gillespie.

I created a World by giving a size when performing calculations with ode and gilespie, but what is the handling of space in E-Cell 4?

from ecell4 import *

w1 = ode.ODEWorld(Real3(1, 1, 1))
w2 = gillespie.GillespieWorld(Real3(1, 1, 1))

In ode and gillespie, we calculated in a cube with one side as in the above example, but in reality only that volume matters.

w3 = ode.ODEWorld(Real3(2, 0.5, 1))  # is almost equivalent to 'w1'
w4 = gillespie.GillespieWorld(Real3(2, 2, 0.25))  # is almost equivalent to 'w2'

However, since 1 is given as the volume after all, the same calculation result will be shown.

This seems reasonable in a system that is actually well agitated and spatially uniform, such as in vitro. However, the cells are clearly not spatially uniform. To take into account the localization of these molecules. Requires space-based biochemical calculations.

Various spatial representations and corresponding calculation techniques are available in E-Cell 4. Below, we will first look at the handling of space in E-Cell 4 using the spatial Gillespie method, which is one of them, as an example.

Spatial Gillespie method

In E-Cell 4, the spatial Gillespie method [^ 1] is included in the'meso'module. First, let's calculate using'run_simulation' as we saw in the beginner's edition without thinking about anything.

[^ 1]: More correctly, it is closer to the Next Subvolume Method. For details, refer to the paper.

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}, solver='meso')

\frac{d\mathrm{C}}{dt}=0.01\frac{\mathrm{A}}{V}\frac{\mathrm{B}}{V}-0.3\frac{\mathrm{C}}{V}=0\\
0.01\left(60-\mathrm{C}\right)^2=0.3\mathrm{C}\times V\\
\mathrm{C}=30

It may have taken some time compared to ode and gillespie, but you'll get almost the same result. That's right, the meso module does almost the same calculation as gilespie without any spatial settings. Let's expand'run_simulation'.

equilibrium3.py

from ecell4 import *

with reaction_rules():
    A + B == C | (0.01, 0.3)

m = get_model()

w = meso.MesoscopicWorld(Real3(1, 1, 1), Integer3(1, 1, 1))  # XXX: Point2
w.bind_to(m)  # XXX: Point1
w.add_molecules(Species('C'), 60)

sim = meso.MesoscopicSimulator(w)  # XXX: Point1
obs = FixedIntervalNumberObserver(0.1, ('A', 'B', 'C'))
sim.run(10, obs)

viz.plot_number_observer(obs)

It's almost the same except that it uses Mesoscopic World and Mesoscopic Simulator, but let's take a look at some new elements.

First, in'w.bind_to (m)', we associate a Model with World. We didn't have to do this in the beginner's edition, but don't forget that the attributes of Species are very important in spatial techniques. Instead, it is enough to give the Mesoscopic Simulator only World.

The next important difference is that when creating the Mesoscopic World, we give'Integer3 (1, 1, 1)'as the second argument. [^ 2]. This wasn't seen in ODE World or Gillespie World. Before explaining what this is, let's change it and see what happens.

[^ 2]: ‘Integer3’ is the same 3D vector as ‘Real3’, but the element is not a real number Real but an integer value Integer.

w = meso.MesoscopicWorld(Real3(1, 1, 1), Integer3(4, 4, 4))  # <- Integer3(1, 1, 1)

I think you got a different result than before. The higher the number, the more noticeable the difference will be.

In fact, this second argument represents the number of divisions in space. Meso is almost the same as gillespie, but while gillespie performs calculations in one uniform closed space, meso Divides the space into small rectangular parallelepipeds called Subvolumes, allowing each Subvolume to have different molecular concentrations. Therefore, in the above example, the cube with one side is divided into 64 cubes with one side of 0.25.

Since 60 molecular species C are input at the beginning, at most 1 C is contained in one Subvolume. At this time, the concentration of C is the same as the first case, but A after separation. The concentration of one A molecule that has just been separated is originally $ 1/1 = 1 $, but in Subvolume it becomes $ 1 / (1/64) = 64 $. With this number of divisions, the effect so far is not obtained, but as a result, the binding reaction becomes faster, so the amount of bound state C increases as shown above. Let's explain from another perspective. In meso Immediately after separation, A and B are always in the same Subvolume. On the other hand, gillespie assumes a well-stirred and uniform system, so A and B molecules immediately after separation are included in one of the 64 Subvolumes. At this time, the probability that A and B are included in the same Subvolume is $ 1/64 $. Therefore, the probability of combining is also high.

Gives the diffusion coefficient of the molecule

The reason why such a difference is made is that we have obtained space by using meso, but on the other hand, we have not yet considered the diffusion motion of molecules.

To set this, use the attribute value'D'of the numerator species Species in the way we did in the beginner's edition. Last time we used add_species_attribute, but like we did in ReactionRule, here is the special notation for E-Cell 4. Is available.

with species_attributes():
    A | {'D': '1'}
    B | {'D': '1'}
    C | {'D': '1'}

    # A | B | C | {'D': '1'}  # means the same as above

By using'species_attributes', you can set attributes with dict type [^ 3] following separator'|'. Here,'D' which means diffusion coefficient is set to 1.

Now let's add this to the above model and recalculate. I think it will take longer than last time, but I think we got the same result as gilespie again.

So why does the diffusion movement of molecules solve the previous problem?

Consider the free diffusion motion of a molecule with a diffusion coefficient of $ D $ in three-dimensional space. The unit of the diffusion coefficient is $ \ mathrm {\ mu m} ^ 2 / s $ or $ \ mathrm {nm} ^ 2 / . It is the square of the length divided by the time like mu s $. At this time, it is known that the square of the distance from the first position to the position after the time $ t $ seconds is equal to $ 6Dt $ on average. Has been.

Conversely, the average time scale for a molecule to stay in a space with a length scale of $ l $ is about $ l ^ 2 / 6D $. In the previous example, one side of the Subvolume is 0.25. Since the diffusion coefficient is 1, it is about 0.01 seconds. On the other hand, even if A and B molecules are in the same Subvolume, it takes about 1.5 seconds to react, so even if they are separated in the same Subvolume, it is almost the same. In the case of, it diffuses before the reaction and moves to another Subvolume. The smaller the $ l $, the smaller the subvolume volume $ l ^ 3 $, and therefore the higher the reaction rate after separation. On the other hand, it takes less time to get out of the space due to diffusion.

Molecular localization

So far, the'add_molecules' function has been used to add molecules to the World, as with ode and gillespie. On the other hand, Mesoscopic World allows the arrangement of molecules according to their spatial representation.

w = meso.MesoscopicWorld(Real3(1, 1, 1), Integer3(3, 3, 3))
w.add_molecules(Species('A'), 120)
w.add_molecules(Species('B'), 120, Integer3(1, 1, 1))

In Mesoscopic World'add_molecules', you can specify which Subvolume to place the molecule in by giving only one Integer3 as the third argument.

In the above example, the molecular species A is scattered throughout the space, but B is located just in the central Subvolume [^ 4]. To confirm this, use'num_molecules'.

[^ 4]: Since the coordinates start at 0, Integer3 (1, 1, 1) becomes the second Subvolume in each axis direction.

print(w.num_molecules(Species('B')))  # will print 120
print(w.num_molecules(Species('B'), Integer3(0, 0, 0)))  # will print 0
print(w.num_molecules(Species('B'), Integer3(1, 1, 1)))  # will print 120

Furthermore, if you are using IPython Notebook, you can grasp the localization of molecules more intuitively by using the ecell4.viz library [^ 5].

viz.plot_world(w, radius=0.01)

'viz.plot_world' displays the particle arrangement interactively on the IPython Notebook by giving World [^ 6]. Set radius as an argument to specify the size of the molecule.

[^ 6]: Molecules are displayed as particles, but in Mesoscopic World, there is no specific position where they actually exist in the Subvolume. Therefore, they are randomly placed in the Subvolume every time.

Since we were able to give localization to the initial arrangement of the molecule, let's actually perform a simulation based on this. In the above example, the diffusion coefficient was set to 1 and one side of the entire World was set to 1, so it is sufficient. 10 seconds will be enough to stir. Call'viz.plot_world'again after the calculation to make sure it is actually agitated.

Initial arrangement and reaction of molecules

Let's take an extreme example to see how the localization of molecules affects the reaction.

from ecell4 import *

with species_attributes():
    A | B | C | {'D': '1'}

with reaction_rules():
    A + B > C | 0.01

m = get_model()
w = meso.MesoscopicWorld(Real3(10, 1, 1), Integer3(10, 1, 1))
w.bind_to(m)

The model contains only simple bonds. World is a rectangular parallelepiped that is long only in the X-axis direction. After that, the molecules are arranged in a biased manner.

w.add_molecules(Species('A'), 1200, Integer3(2, 0, 0))
w.add_molecules(Species('B'), 1200, Integer3(7, 0, 0))
viz.plot_world(w, radius=0.025)

As an aside, there is a reason not to use'Integer3 (0, 0, 0)'and'Integer3 (9, 0, 0)' here. Basically, all Worlds in E-Cell 4 are periodic boundary conditions [^ Periodic Boundary]. Therefore, the previous two Subvolumes are actually next to each other.

[^ Periodic Boundary]: periodic boundary condition. Simply put, when focusing on a certain axis, when the molecule exceeds the right end, the left end comes out, and vice versa.

If the expected arrangement is obtained, let's calculate using the Mesoscopic Simulator.

sim = meso.MesoscopicSimulator(w)
obs1 = NumberObserver(('A', 'B', 'C'))  # XXX: saves the numbers after every steps
sim.run(5, obs1)
viz.plot_number_observer(obs1)

Did the reaction proceed properly? To check the result, try to distribute the arrangement uniformly with meso, or let the gillespie module calculate the same model.

When the solid line is biased and there is no dotted line. When there is a bias, there is a clear delay in the reaction. If you look closer, you may notice that the shape of the time series is also different between the solid line and the dotted line. No. This is probably because it takes time for A and B molecules to meet by diffusion. In fact, the time scale required to travel the distance between A and B (approximately 4) due to the initial arrangement is $ 4 ^ 2/2. (D_ \ mathrm {A} + D_ \ mathrm {B}) = 4 $ seconds.

Single-molecule particle size simulation with Spatiocyte

We have explained an example of spatial representation in E-Cell 4 using the spatial Gillespie method. Next, let's look at the higher resolution spatial representation '1 molecular grain size calculation'.

Microscopic lattice method Spatiocyte

In the spatial Gillespie method, the space normally handled by the Gillespie method is divided into smaller spaces, and the space is introduced by considering the diffusion of molecules between the subvolumes. However, the number of molecules in the subvolume is still large. It is treated and the position of each molecule is uncertain. In other words, the spatial resolution of this calculation method is equal to $ l $ on each side of Subvolume. To increase this resolution, $ l $ is added more and more. You can make it smaller, but in this method $ l $ must be large enough (at least 10 times larger) than the molecule diameter $ R $.

Then, how can we increase the spatial resolution to the size of the molecule? The answer is not the number of molecules, but the one-molecule particle size simulation that directly expresses the reaction-diffusion motion of each molecule with a position. E-Cell 4 can use multiple single-molecule particle size calculation techniques, but here we will explain using the microscopic lattice method Spatiocyte as an example [^ Arjunan10] [^ Spatiocyte].

Spatiocyte treats each molecule as one reactive rigid sphere and moves it in space represented by a hexagonal close-packed lattice [^ HCP]. Unlike the spatial Gillespie method, each molecule is given an identification number. It is possible to know where the molecule is by the particle size of the molecule. In terms of time scale, the time scale by diffusion is proportional to the square of the distance, so it is 100 times finer than the spatial Gillespie method. It will be the time width.

[^ HCP]: [http://ja.wikipedia.org/wiki/Hexagonal close-packed structure](http://ja.wikipedia.org/wiki/%E5%85%AD%E6%96%B9% E6% 9C% 80% E5% AF% 86% E5% 85% 85% E5% A1% AB% E6% A7% 8B% E9% 80% A0)

Anyway, let's use Spatiocyte. However, the basics are exactly the same.

Note: In the following examples, the name lattice is used for Spatiocyte, but it will be renamed to the module name spatiocyte in a future change .

equilibrium4.py

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)

The difference is the second argument given to LatticeWorld, which is called the Voxel radius. In Spatiocyte, the space is divided into molecular sizes to determine the position of the molecule, but the smallest unit of this space is Voxel. In most cases, this can be the size of the molecule. In this example, the radius of the molecule is given in the space of 1 $ \ mathrm {\ mu m} $ on each side. 5 $ \ mathrm {nm} $.

The calculation time will be longer than before, but the results will of course be consistent with ode and gillespie.

Diffusion movement of one molecule

Let's actually deal with one molecule to confirm the one-molecule particle size calculation.

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

The'new_particle'function puts one molecule in LatticeWorld at the specified location. At the same time, it returns the identification number'pid' of the created molecule, the tuple of the particle information'p', and'suc' whether it was successfully placed. If the molecule is already placed in that place, it cannot be placed on top of each other, so'suc' becomes False and fails. Particle'p' contains the position, molecular species, radius, and diffusion coefficient of the molecule. From now on, the particle'p'can be extracted and the information can be obtained by using the identification number'pid'.

Let's check it out.

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)

The'get_particle'function receives the identification number and returns the identification number and the particle. Naturally, the identification number is the same as the one given. From the particle, the position can be retrieved as Real3 by the'position' function. It is difficult to read as it is. It is displayed after being converted to tuple.

The returned position is slightly off the specified location because the Spatiocytes can only place the molecules on a grid. LatticeWorld will place the molecules in the grid closest to the given location.'Viz If you use the .plot_world'function, you can see that only one molecule is certainly centered.

The Observer can also be used to track the diffusion of this molecule.

sim = lattice.LatticeSimulator(w)
obs = FixedIntervalTrajectoryObserver(0.002, (pid,))
sim.run(1, obs)
viz.plot_trajectory(obs)

Here, it is displayed using the'viz.plot_trajectory'function, but it is also possible to extract the trajectory as a list of Real3 by the'data' function as in the case of NumberObserver.

print(len(obs.data())) # should return 1
print(len(obs.data()[0])) # should return 501

The data function returns a nested list. The index in the first list is the index of the particle. Since only one particle is created here, its index length is 1. The index in the next list is the Real3 list that corresponds to the particles in the first list. Here, the simulation was run for 1 second, and the trajectory was recorded in 0.002 second increments. Therefore, 501 Real3s are returned, including the initial value of one Real3 and 1 / 0.002 = 500 Real3s.

Also, although the identification number of the molecule added by the'add_molecules' function is unknown, particles can be collectively extracted from the molecular species by'list_particles'.

w.add_molecules(Species('A'), 5)

particles = w.list_particles(Species('A'))
for pid, p in particles:
    print(p.species().serial(), tuple(p.position()))

The'list_particles' function can be used in other worlds as well as the'add_molecules' function [^ list_particles] It is important to remember. As an aside, the formal function that handles one molecule in Spatiocyte is originally' It is list_voxels', and the coordinates are represented by one integer value instead of Real3.

[^ list_particles]: However, it is not available in ODE World.

Diffusion coefficient and secondary reaction

The binding reaction of the model that we have dealt with so far is also called a secondary reaction. Let's look at the relationship between this secondary reaction and the diffusion coefficient in Spatiocyte.

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

This time, we use a higher rate constant than before, but we get the same result as before. However, a detailed comparison with the result of the ordinary differential equation shows that the result is different (the solid line in the figure below is the Spatiocyte). , The dotted line is due to ode).

Is this something wrong? In fact, when the rate constant is set to a much higher value, the reaction rate stops endlessly with ode, but it is hardly faster with Spatiocyte.

This is due to the difference in the definition of the reaction rate constant between the solver and the single molecule particle size calculation technique. The former is called macroscopic or effective effective reaction rate constant, while the latter. Is called the microscopic, or intrinsic intrinsic reaction rate constant.

The macroscopic reaction rate constant represents the rate at which the reaction occurs when the molecules are mixed, whereas the microscopic reaction rate constant indicates the degree of reactivity when the molecules actually collide. Therefore, in a microscopic view, it must first collide before it reacts. In Spatiocyte, no matter how fast this microscopic rate constant is, it is more than the speed of collision caused by diffusion. This situation is called diffusion rate constant. This is very similar to the fact that it took time for the molecules arranged unevenly to react with each other in the spatial Gillespie method.

The following relationship is known in three-dimensional space between this macroscopic rate constant $ k_ \ mathrm {on} $ and the microscopic rate constant $ k_a $ [^ Collins Kimball].

[^ CollinsKimball]: Smoluchowski–Collins–Kimball formula.

\frac{1}{k_\mathrm{on}}=\frac{1}{k_a}+\frac{1}{4\pi RD_\mathrm{tot}}

Where $ R $ is the sum of the radii of the two colliding molecules and $ D_ \ mathrm {tot} $ is the sum of the diffusion coefficients. In the above example, the second term on the right side $ k_D = 4 \ pi RD_ \ Since mathrm {tot} $ is about 0.25, the macroscopic reaction rate constant is about 0.2 when combined with the microscopic reaction rate constant 1.0 [^ alpha].

[^ alpha]: However, in Spatiocyte, the secondary reaction rate constant must be slower than $ 3 \ sqrt {2} RD $, and the collision rate constant $ k_D $ is also $ 3 \ sqrt, unless special settings are made. {2} RD $.

Unlike ordinary differential equations and the Gillespie method, which assume such a well-stirred system (that is, the diffusion coefficient is infinite), the one-molecule particle size calculation can separate the diffusion and reaction of molecules properly. As can be seen from the equation of, if the microscopic velocity constant is sufficiently small with respect to $ k_D $, it almost matches the macroscopic velocity constant (reaction rate-determining).

Treatment of structures in Spatiocyte

Finally, I will explain how to handle structures such as cell membranes using Spatiocyte. In E-Cell 4, the structure is still in the development stage, but Spatiocyte can use its function to some extent.

Currently, the shapes that can be handled are limited, but let's first look at a sphere as an example. To limit the movement of molecules to a spherical space, we first create a spherical structure.

w = lattice.LatticeWorld(Real3(1, 1, 1), 0.005)
sph = Sphere(Real3(0.5, 0.5, 0.5), 0.45)
print(w.add_structure(Species('C'), sph))  # will print 539805
viz.plot_world(w)

'Sphere' means a sphere with the first argument as the center and the second argument as the radius. I gave it the Species'C' and added it to the World.

In Spatiocyte, the structure is represented by filling the area in space with a special Voxel. In the above example, all Voxels in the specified sphere are occupied by the molecular species C.'viz. If you try to display it with plot_world', you can see that it is actually distributed in a spherical shape. However, since the number of molecules is too large to display as about 540,000, only a part is plotted, but actually Is completely filled.

Once the sphere is created, create a molecule that moves around only in this sphere. To do this, specify the'location'attribute.

with species_attributes():
    A | {'D': '1', 'location': 'C'}

m = get_model()

All you have to do is'add_molecules' the molecules as usual.

w.bind_to(m)
w.add_molecules(Species('A'), 120)
viz.plot_world(w, species_list=('A',))  # visualize A-molecules only

Now that we have specified that molecular species A exists on top of the structure of molecular species C, even'add_molecules' is actually placed only in the sphere. Note that, of course, the structure before adding A. Need to be made.

Let's check if the movement of this molecular species A is really restricted in the sphere. Again,'FixedIntervalTrajectoryObserver' can be used.

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)

'pid_list' is a list of 10 suitable identification numbers out of 60 A molecules. Colored the diffusion trajectory of these 10 molecules. Certainly, movement is restricted in the sphere. I understand.

Structure and reaction

Finally, let me explain the transition of molecules between structures. First, molecular species without the'location'attribute do not belong to any structure. In the previous example, specify'location'. Otherwise it will be placed outside the sphere.

Let's take a plane as an example and try it. First, to make a plane, we can make it as follows using three Real3s, the origin'origin'and the two axis vectors'unit0' and'unit1'.

ps = PlanarSurface(origin, unit0, unit1)

Using this, assuming molecular species A on a plane and normal molecular species B,

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

It is a little difficult to understand, but B molecule is arranged only on the plane. Then, how should we express the reaction that this B molecule is adsorbed on the plane M and becomes A molecule?

with reaction_rules():
    B + M == A | (1e-3, 1.5)

As you can see, when the B molecule collides with the structure M, it is adsorbed and becomes the A molecule. On the contrary, the A molecule becomes the original plane and the B molecule when separated.

After that, you can calculate with Simulator as usual.

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

By the way, the molecular species of the structure can be omitted by separating it from the structure, but it cannot be omitted by binding. It is impossible to make A from B in the situation where there is no M around which the product A molecule should be placed. That's why. On the contrary, in this case, there is no problem because the A molecule can place the B molecule anywhere on the plane M. Furthermore, in the case of separation, the structure is primary or not omitted. It is still a reaction, but in binding, the secondary reaction becomes a primary reaction, and the meaning of the rate constant also changes.

with reaction_rules():
    B + M > A | 1e-3
    A > B | 3.0  # means the same as A > B + M

Exercises

1. Exercise on the spatial Gillespie method: We calculated the time scale of diffusion for the case of Subvolume with one side of $ l $ in the spatial Gillespie method. On the other hand, we also explained the time scale for the molecules immediately after separation to rejoin. Consider the dependence of these two time scales on $ l $ and clarify the conditions that $ l $ should satisfy. However, it is known that the spatial Gillespie method should satisfy the following conditions for $ l $. Yes.

R^2 \ll l^2 \ll 6D\tau_\mathrm{min}

 Where $ R $ is the diameter of the molecule (more precisely, the radius of reaction) and $ \ tau_ \ mathrm {min} $ is the shortest life-time for each reaction in the molecule [^ Elf04].

[^Elf04]: Elf J. and Ehrenberg M., "Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases", Syst. Biol., 1(2), 230-236 (2004)

 1. Exercises on the initial arrangement of molecules: Modify the above example to investigate the effects on the reaction by changing various parameters such as the initial arrangement of AB molecules, reaction rate constant, and diffusion coefficient. Is it possible to create a model that reproduces the time series like the above example using the uniform Gillespie method?

 1. Exercise on molecular diffusion: In the case of free diffusion in three dimensions, the mean square distance to the position of the molecule with the diffusion coefficient $ D $ after $ t $ seconds is Mean Square Displacement.

    ```math
\left<L^2\right>=6Dt

Is known to be the case. Is it actually the case with Spatiocyte?

1. Exercise on macroscopic and microscopic rate constants: Macroscopic and microscopic rate constants have the same unit of $ \ mathrm {nM} ^ {-1} s ^ {-1} $ and $ \ mathrm {\ mu m} ^ 3s ^ {-1} $. Two molecules with a diffusion rate of 1 $ \ mathrm {\ mu m} ^ 2 / s $ and a diameter of 5 $ \ mathrm {nm} $ What is the collision rate constant $ k_D $ $ \ mathrm {nM} ^ {-1} s ^ {-1} $?

2. Structure practice: Using a plane structure, check the relationship between the diffusion coefficient $ D $ on the plane and its mean square distance.