[PYTHON] Easy to use E-Cell 4 Advanced Edition
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 advanced edition. It is assumed that you understand the beginner and intermediate editions.
-Easy to use E-Cell 4 for beginners -Easy to use E-Cell 4 Intermediate
In addition, this article assumes the use of IPython Notebook.
How to use
I feel that it is no longer "easy" to use, but I would like to explain a little more advanced features at the end.
Rule-based modeling
Give the molecular species a structure
Previous molecular species names were just symbols.
For example, when considering the reaction A> B, in many cases, the species A becomes the molecule type B due to some modification state or structural change, and it is the same molecule as the molecule. It is possible. However, in the case of molecular species A and B, the relationship between the molecular species is not clear. Let me give another example. In the reaction A + A> B, molecular species B is probably It can be expected that the A molecule is a dimer with two molecules bonded, but it is not known from the name B, and even if the name'B'is changed to the name'A2', the computer cannot judge it.
To solve these problems, giving a structure to the molecular species is the first step in rule-based modeling. Below, let's actually experience rule-based modeling using E-Cell 4.
with reaction_rules():
X(s=u) > X(s=p) | 1 # corresponding to 'A > B | 1'
In E-Cell 4, the structure of the molecular species is determined using parentheses and equal signs. In the above example, `X (s = u)` and `X (s = p)` Each represents a molecular species name. Here, the molecular species `X (s = u)` means that the modification's' of the molecule'X'is in the state'u'. Similarly, the molecular species` `` X (s = p) ``` is that the modified's' of the molecule'X'is the state'p'. Thus, in this reaction, the molecule'X' is the modified's' state of'u'. Means to change from'p'. This is the structure of the molecular species.
Let's use World to confirm this.
from ecell4 import *
w = gillespie.GillespieWorld(Real3(1, 1, 1))
w.add_molecules(Species('X(s=u)'), 60)
w.add_molecules(Species('X(s=p)'), 120)
print(w.num_molecules(Species('X(s=u)'))) # will print 60
print(w.num_molecules(Species('X(s=p)'))) # will print 120
The number is counted by adding 60 and 120 molecular species, respectively. Even if the writing method of the molecular species is changed, the added number can be obtained as usual. Then, what happens if the following is done?
print(w.num_molecules(Species('X'))) # will print 180
I counted the number of the molecular species'X', but we did not add any of the molecular species'X'. Nevertheless, we got the number of 180. This is called'X'. This is because all the molecules named'X'are counted regardless of the modified state, not the molecular species named. The first two molecular species added are both named'X', only the modified state is different. Therefore, both are true when counting'X'.
By the way, if you want to count only the molecular species named'X', use `` `num_molecules_exact```. This time, there is no such molecular species, so you should get 0.
Next, consider the bond. The molecular species representing the complex in which the molecule'X'and the molecule'Y' are bonded can be written as follows.
with reaction_rules():
X(b1) + Y(b2) > X(b1^1).Y(b2^1) | 1
Look at the right-hand side of this reaction. You can see that'X'and'Y' are connected by dots'.'. In this way, the complex is represented by connecting multiple molecular species with dots. ..
By the way, X (b1 ^ 1) .Y (b2 ^ 1)` `` represents one molecular species as a whole, but some of them are'X (b1 ^ 1)'and'Y (b2 ^ 1). You can see that it contains two elements called)'. In E-Cell 4, these are called the unit molecular species Unit Species. X (b1 ^ 1) .Y (b2 ^ 1) `` is two It is a molecular species containing one unit molecular species. Of course, if it is ``` X (p = u) ``, it is made up of one unit molecular species.
Returning to the story, the unit molecule species of ``` X (b1 ^ 1) .Y (b2 ^ 1)` `` are modified as'b1'and'b2', respectively. The connection points between species are determined.'X' and'Y' are connected by the respective modifiers'b1' and'b2'. The number after the symbol'^' is a number indicating the correspondence of the connection points. Yes. This can be any number greater than 0, as long as it corresponds. Also, this time I gave the modifiers different names such as'b1'and'b2', but they can be the same.
Qualified with X (b1)` `` on the left side The fact that neither the state by the equal sign nor the connection by'^'is written in'b1' means that the connection point'b1'is not connected to anyone. To clarify. X (b1)` `` means'X' which is not bound to anyone through the modifier'b1'.
Now that we've explained the representation of the complex, let's try it in the same way as before.
w.add_molecules(Species('X(b1)'), 10)
w.add_molecules(Species('Y(b2)'), 20)
w.add_molecules(Species('X(b1^1).Y(b2^1)'), 30)
print(w.num_molecules(Species('X'))) # will print 40
print(w.num_molecules(Species('Y'))) # will print 50
print(w.num_molecules(Species('X.Y'))) # will print 30
Last but not least, in the case of a molecular species that has both a binding site and a modification site, it can be written as `X (a, b = c) ```, but X (b = c, a) `` It cannot be written as . Even if there are multiple parts, it must be written as ``` X (a, b, c, d = e, f = g) `. However, the order is changed` ` `X (c, a, b, f = g, d = e)` can be used with the same meaning as before.
Reaction rule
In the previous example, it became possible to count the number of molecules by giving a structure to the molecular species. For example, the molecular species'X'given when counting does not mean a specific molecular species, but counts. It is a kind of pattern that determines the specific molecular species conditions that should be, that is, those that include X in the unit molecular species. The rule-based modeling is to make use of such patterns when describing reactions. The next step.
Consider a simple example consisting of one unit molecule species X.
For example, the unit molecular species X has three phosphorylation sites's1','s2', and's3', each of which can independently take either of two states,'u' or'p'. At this time, the number of states that the unit molecule species X can take is $ 2 ^ 3 = 8 $.
Now, how should we write the phosphorylation reaction of such molecule X? Actually, it is as follows.
with reaction_rules():
X(s1=u, s2=u, s3=u) > X(s1=p, s2=u, s3=u) | 1
X(s1=u, s2=p, s3=u) > X(s1=p, s2=p, s3=u) | 1
X(s1=u, s2=u, s3=p) > X(s1=p, s2=u, s3=p) | 1
X(s1=u, s2=p, s3=p) > X(s1=p, s2=p, s3=p) | 1
X(s1=u, s2=u, s3=u) > X(s1=u, s2=p, s3=u) | 2
X(s1=p, s2=u, s3=u) > X(s1=p, s2=p, s3=u) | 2
X(s1=u, s2=u, s3=p) > X(s1=u, s2=p, s3=p) | 2
X(s1=p, s2=u, s3=p) > X(s1=p, s2=p, s3=p) | 2
X(s1=u, s2=u, s3=u) > X(s1=u, s2=u, s3=p) | 3
X(s1=u, s2=p, s3=u) > X(s1=u, s2=p, s3=p) | 3
X(s1=p, s2=u, s3=u) > X(s1=p, s2=u, s3=p) | 3
X(s1=p, s2=p, s3=u) > X(s1=p, s2=p, s3=p) | 3
It is not easy to grasp what this means and write it without mistakes. The reason why it is so troublesome is that each phosphorylation reaction is not independent, and one phosphorylation site is considered. Because all other states must be counted.
Therefore, let us consider the pattern of molecular species. When focusing on the phosphorylation site's1', the phosphorylation reaction is that the state of's1' of molecule X changes from'u'to'p'. That is, It suffices to take a molecule that fits the pattern of the molecule X whose state is'u'and change its state to'p'.
with reaction_rules():
X(s1=u) > X(s1=p) | 1
X(s2=u) > X(s2=p) | 2
X(s3=u) > X(s3=p) | 3
These are not reactions that have individual specific molecular species on the left and right, but are called reaction rules because they give a pattern of molecular species. In order to use this model, reactions are different from the past. It must be indicated that it is a model by law.
Furthermore, since the reaction rule is only a rule based on the pattern, it is not possible to know what kind of molecular species actually exists. In the above example, the molecular species'X (s1 = u)' appears on the left side. That doesn't mean it really exists. Remember that the first thing we envisioned was a molecular species with three modifications, such as'X (s1 = u, s2 = u, s3 = u)'. want.
In the rule-based model, by giving a molecular species that will be the seed seed, all possible molecular species will be generated in a worm-like manner by applying the reaction rule repeatedly. Think of it as a molecular species to add to.
m = get_model(seeds=(Species('X(s1=u, s2=u, s3=u)'),))
print(len(m.species_attributes())) # will print 8
print(len(m.reaction_rules())) # will print 12
From the species'X (s1 = u, s2 = u, s3 = u)' given here, it can be seen that 8 molecular species and 12 reactions can be obtained as desired.
Finally, let's calculate using this model.
y = run_simulation(
numpy.linspace(0, 5, 100),
{'X(s1=u,s2=u,s3=u)': 120},
model=m,
species_list=('X', 'X(s1=p)', 'X(s2=p)', 'X(s3=p)'))
The molecular species to be recorded are determined by specifying'species_list'. Actually, if you look only at the phosphorylation states of the modified's1','s2', and's3', you can see that they increase by the specified rate constants. Of course, the total amount of'X'molecules is constant.
Wildcard
Now that you know the basics, there's one more thing to remember: wildcards, represented by the symbol'\ _'in E-Cell 4.
The first and most important use is for binding. For example, when the molecular species'A'has a binding site'b', if you want to know the number of A's that are not bound to anyone, ```A (b) `` `. So how do we express A, which is bound to someone?
with reaction_rules():
A(b) + B(b) > A(b^1).B(b^1) | 1
A(b^_, s=u) > A(b^_, s=p) | 2
m = get_model(seeds=(Species('A(b, s=u)'),))
Notice the second reaction rule. That's what ```A (b ^ , s = u) `` `is written.'' Is a wildcard, that is, what. The second reaction rule changes the modification's' from'u'to'p', which means that this reaction only occurs when A is bound to someone. Yes.
If it is A that is simply combined with someone, write ```A (b ^ _) `` `.
This wildcard can of course be used for qualifications as well as joins. A simple example is ```A (s = _) `` `, which means that the qualification's' can be in any state. I don't think it's necessary to write it if it's okay, but I guarantee that A has the qualifier's'.
Another important thing is when it is used for the unit molecule species name. If you write `` `_ (s = u) ```, it has the modifier's', and any molecule in the state'u' Let's give it a try.
wildcard1.py
from ecell4 import *
with reaction_rules():
_(s=u) > _(s=p) | 1
m = get_model(seeds=(
Species('A(s=u)'), Species('B(s=u)'), Species('C(b^1,s=u).C(b^1,s=u)')))
y = run_simulation(
numpy.linspace(0, 5, 100),
{'A(s=u)': 30, 'B(s=u)': 40, 'C(b^1,s=u).C(b^1,s=u)': 60},
model=m,
species_list=('A(s=u)', 'B(s=u)', 'C(s=u)', '_(s=u)'))
It can be seen that the unit molecule species name applies to either'A',' B', or'C'. If a complex is formed, it also applies to each of them. Of course, it can also be used with Observer and num_molecules.
Wildcards can also be used to determine the attributes of a molecular species.
with species_attributes():
_ | {'D': '1', 'radius': '0.005'}
You can specify the attributes of all molecular species by using wildcards without specifying the same attributes for each molecular species. If you want to add an exception, specify it before this.
Make a video
I could use `` `viz.plot_world``` to easily visualize the state of the World. Here's how to make a video. There are several ways to make a video, but this time Let's use ParaView.
First, download and install ParaView from the above site. Below, ParaView-4.3.1-Windows-32bit.exe is used.
Record the position of the particles
First, record the data for visualization. Observer can also be used for this.
This time, let's use the following model.
from ecell4 import *
with species_attributes():
A | {'D': '1', 'location': 'C'}
B | {'D': '1', 'location': 'N'}
N | {'location': 'C'}
with reaction_rules():
A + N > B | 1e-4
B + C > A | 1e-4
m = get_model()
w = lattice.LatticeWorld(Real3(1, 1, 1), 0.005)
w.bind_to(m)
w.add_structure(Species('C'), Sphere(Real3(0.5, 0.5, 0.5), 0.441))
w.add_structure(Species('N'), Sphere(Real3(0.5, 0.5, 0.5), 0.350))
w.add_molecules(Species('A'), 720)
I am thinking of a model in which a nested sphere is created as a structure and flows from the outside to the inside. The inner sphere N is placed so as to overlap the outer sphere C.
To perform the calculation for the next 1 second, select the molecular species'A'and'B' every 0.01 seconds and record their positions in CSV format. However, the procedure is basically the same as before.
sim = lattice.LatticeSimulator(w)
obs1 = NumberObserver(('A', 'B'))
obs2 = FixedIntervalCSVObserver(0.01, 'sample%03d.csv', ('A', 'B'))
sim.run(1, (obs1, obs2))
This time, we are using two Observers at the same time. Of these, we use `` `FixedIntervalCSVObserverfor recording in CSV format. The arguments are the time width to record, the file name, and the list of molecular species. You can include the file number in the file name by using% 03d```.
After doing the above, you should be able to create from'sample000.csv'to'sample100.csv'. Save the file in the location where you ran Python, or if you are using IPython Notebook, under the directory where you started IPython Notebook on the server side. Should be done.
Visualization with ParaView
For the explanation of ParaView itself, let's directly refer to the following site, but here I will explain only the visualization of the result of E-Cell 4.
First, select the CSV format data saved with File-> Open. The saved serial number files are displayed together as shown in the figure below. After loading, click the Apply button in Properties at the bottom left. Press to complete loading.
After reading, visualize the data with Pipeline. First convert from CSV table to coordinates. Select Filters-> Alphabetical-> Table To Points. From the bottom left, X, Y, Z After selecting x, y, z as Column, press the Apply button.
Finally, select Filters-> Common-> Glyph. As before, change `Glyph Type``` from` `Properties``` to` `Sphere``` in the lower left. Set Scale Mode``` to `` scalarand press the button with the hint Reset using current data values `next to ``` Scale Factor``` below. I think that some numerical value will be set, but ignore it and press the Apply button after entering `4``` this time.
Then, the item "Coloring" will be added to the previous Properties, so select sid. You should be able to visualize as below. If you can not see anything, the upper left Pipeline Browser Make sure that the left eyeball of `Glyph1``` is active in`.
After that, press the `Play``` button above. If you can see how the red color increases properly, you are done. You can save the video from File-> Save Animation ... `` must.
Other tips
Use of Factory
I will explain how to write to make it easier to switch modules of E-Cell 4. If you are not good at programming, you can skip it.
You may have noticed that the difference when using different modules of ode, gillespie, meso, and lattice appears almost only when creating World. However, in reality, the name of Simulator is also different, and when creating World, 2 There may or may not be a second argument. To solve these problems, each module provides a convenience class called Factory.
f = ode.ODEFactory()
# f = gillespie.GillespieFactory()
# f = meso.MesoscopicFactory(Integer3(4, 4, 4))
# f = lattice.LatticeFactory(0.005)
w = f.create_world(Real3(1, 1, 1)) # XXX: Factory
w.bind_to(m)
w.add_molecules(Species('C'), 60)
sim = f.create_simulator(w) # XXX: Factory
obs = FixedIntervalNumberObserver(0.1, ('A', 'B', 'C'))
sim.run(10, obs)
By doing this, you can switch modules by switching Factory. Furthermore, let's make this a function so that you can easily compare.
comparison.py
from ecell4 import *
def run(m, f):
w = f.create_world(Real3(1, 1, 1)) # XXX: Factory
w.bind_to(m)
w.add_molecules(Species('C'), 60)
sim = f.create_simulator(w) # XXX: Factory
obs = FixedIntervalNumberObserver(0.1, ('A', 'B', 'C'))
sim.run(10, obs)
return obs
with species_attributes():
A | B | C | {'D': '1'}
with reaction_rules():
A + B == C | (0.01, 0.3)
m = get_model()
obs1 = run(m, lattice.LatticeFactory(0.005))
obs2 = run(m, ode.ODEFactory())
viz.plot_number_observer(obs1, '-', obs2, '--')
Of course, not all functions of World can be used in common, but using Factory makes it easier to separate the algorithm from the problem.
Exercises
-
Practice of molecular species in rule-based modeling. How to express a dimer consisting of two X molecules. Also, try how to handle when counting the number of X.
-
Practice of reaction rules in rule-based modeling. When there are two types of unit molecular species'X'and'Y' and the following reactions are established independently, how many molecular species exist?
- Xs form a dimer.
- X and Y combine.
- Y has one modification site and takes one of two states.
However, pay attention to the symmetry.
-
Practice on recording the position of the molecule. This time, I recorded it in a CSV format file with `` `FixedIntervalCSVObserver```, but when I opened it with a text editor or spreadsheet software, what was actually written Let's check if there is.
-
Practice on visualization with ParaView. ParaView can change its appearance in various ways depending on the settings, so make your own video.
Recommended Posts