notebook.community

Edit and run

Lecture 5. その他の高度な機能

海津一成





In [1]:



    


%matplotlib inline
from ecell4 import *
from ecell4_base.core import *
from ecell4_base import *



    











In [2]:



    


import matplotlib as mpl
inline_rc = dict(mpl.rcParams)
import seaborn
seaborn.set(font_scale=1.5)
# mpl.rc("figure", figsize=(8.0, 5.5))
mpl.rc("figure", figsize=(6, 4))



    


















    






/home/kaizu/src/ecell4/local/lib/python3.6/site-packages/matplotlib-3.0.2-py3.6-linux-x86_64.egg/matplotlib/__init__.py:886: MatplotlibDeprecationWarning: 
examples.directory is deprecated; in the future, examples will be found relative to the 'datapath' directory.
  "found relative to the 'datapath' directory.".format(key))

構造体の形をつくる

E-Cell4は基本的な構造体の形をいくつかサポートしている。

  • PlanarSurface: 平面
  • Sphere: 球体
  • Cylinder: 円柱
  • Rod: ロッド
  • AABB: 直方体




In [3]:



    


help(Rod.__init__)



    


















    






Help on wrapper_descriptor:

__init__(...)
    Constructor.
    
    Parameters
    ----------
    length : float
        The length of a cylinder part of a rod.
    radius : float
        The radius of a cylinder and sphere caps.
    origin : Real3, optional
        The center position of a rod.


















In [4]:



    


w = spatiocyte.World(Real3(2, 1, 1), 0.01)
w.add_structure(Species('C'), Rod(1, 0.5, Real3(1, 0.5, 0.5)))
viz.plot_world(w)

これらの基本要素に対してその表面だけを構造体とすることもできる。





In [5]:



    


w = spatiocyte.World(Real3(2, 1, 1), 0.01)
obj = Rod(1, 0.5, Real3(1, 0.5, 0.5))
w.add_structure(Species('C'), obj.surface())
show(w)

これらを組み合せてより複雑な形状をつくることもできる。

  • Union: ふたつのShapeの和で新しいShapeをつくる
  • Complement: ふたつのShape、XとYを受けとり、その差（X-Y）で新しいShapeをつくる
  • AffineTransformation: 受けとったShapeに対して変形操作を行う。今は平行移動とリサイズ、軸周りの回転が可能




In [6]:



    


w = spatiocyte.World(Real3(1.75, 1, 1), 0.01)
obj1 = Sphere(Real3(0.5, 0.5, 0.5), 0.5)
obj2 = Sphere(Real3(1.25, 0.5, 0.5), 0.5)
obj = Union(obj1, obj2)
w.add_structure(Species('C'), obj.surface())
show(w)

ルールベース・モデリング

分子種に構造を持たせる

これまでの分子種の名前は単なる記号であった。

例えば、A > Bという反応を考えたとき、多くの場合、分種Aがなんらの修飾状態や構造の変化によって分子種Bになるのであって分子としては同じ分子であると考えられる。 しかし分子種AやBといった場合、その分子種間の関係性は明らかではない。

A + A > Bという反応ではおそらく分子種BはA分子が2分子結合した二量体だろうと予想できるがBという名前からそれはわからない。

こうした問題に対して、分子種に構造を与えるのがルールベース・モデリングの第一歩である。





In [7]:



    


with reaction_rules():
    # corresponding to 'A > B | 1'
    X(s=u) > X(s=p) | 1



    











In [8]:



    


m = get_model()

E-Cell4では括弧や等号を使って分子種の構造を決める。

  • 分子種X(s=u)はXという分子の修飾sが状態uである。
  • 分子種X(s=p)はXという分子の修飾sが状態'p'である。
  • この反応は、分子Xの修飾sの状態をuからpにする。

分子種の構造をWorldで確かめる。





In [9]:



    


w = gillespie.World()
w.add_molecules(Species('X(s=u)'), 60)
w.add_molecules(Species('X(s=p)'), 120)
print(w.num_molecules(Species('X(s=u)')))
print(w.num_molecules(Species('X(s=p)')))
print(w.num_molecules(Species("X")))



    


















    






60
120
180

Xという分子種は1つも追加していないにもかかわらずnum_moleculesで180個という数が得られた。

これはXという分子を修飾状態に関わらず全て数え上げているためである。 加えた2つの分子種はどちらもXという分子で修飾状態だけが異なる。従って、Xを数えた場合どちらも当てはまる。

きっちりXという名前の分子種だけを数えたい場合はnum_molecules_exactを使う。





In [10]:



    


print(w.num_molecules_exact(Species('X(s=u)')))
print(w.num_molecules_exact(Species('X(s=p)')))
print(w.num_molecules_exact(Species("X")))



    


















    






60
120
0

次に、Xという分子とYという分子が結合した複合体を表す分子種は以下のように書ける。





In [11]:



    


with reaction_rules():
    X(b1) + Y(b2) > X(b1^1).Y(b2^1) | 1



    











In [12]:



    


m = get_model()

XとYがドット.で結合されている。このように複合体は複数の分子種をドットでつないで表現する。

X(b1^1).Y(b2^1)は全体として一つの分子種を表すが、中にX(b1^1)とY(b2^1)という二つの要素が含まれる。E-Cell4ではこれらを単位分子種UnitSpeciesと呼ぶ。

X(b1^1).Y(b2^1)の単位分子種はb1、b2という修飾を持つ。

これが単位分子種間の結合箇所を決める。XとYは各々の修飾b1、b2を介して結合している。

記号^の後の数字は結合箇所の対応を表す数字である。これは対応さえ取れていれば、0より大きなどんな数字でも良い。

左辺のX(b1)で修飾b1に等号による状態も^による結合も記されていないのは、結合箇所b1が誰とも結合していないことを明らかにするため。

X(b1)は修飾b1を介して誰とも結合していないXを意味している。

複合体の表現を説明したところで先ほどと同じように試してみる。





In [13]:



    


w = gillespie.World()



    











In [14]:



    


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')))
print(w.num_molecules(Species('Y')))
print(w.num_molecules(Species('X.Y')))



    


















    






40
50
30

XとYは単体と複合体の和になっているのに対して、X.Yは複合体だけを数えているのがわかる。

注意点だが、結合部位と修飾部位を両方もつ分子種の場合、

X(a, b=c)とは書けるがX(b=c, a)とは書けない。

部位が複数ある場合もX(a, b, c, d=e, f=g)のように書く必要がある。ただし、順番を入れ替えたX(c, a, b, f=g, d=e)は先ほどと同じ意味になる。

反応則

先の例では分子種に構造を与えることで分子数を数えることができるようになった。

例えば、数えるときに与えた分子種Xは具体的な分子種ではなく、数えるべき分子種の条件（「単位分子種にXを含むもの」）を決めた一種のパターンである。

このようなパターンを反応に用いるのが、ルールベース・モデリングの次のステップである。

単位分子種Xには3つのリン酸化部位s1、s2、s3があり、各々uかpの二状態のいずれかを独立に取るものとする。

このとき単位分子種Xのとりえる状態数は$2^3=8$状態となる。

このような分子Xのリン酸化反応は次のように書ける。





In [15]:



    


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



    











In [16]:



    


m = get_model()

このように面倒になるのは、各リン酸化反応が独立であることが活かされず、1つのリン酸化部位を考えるのに他の全ての状態を数え上げなければならないからである。

そこで分子種のパターンで考える。

リン酸化部位s1に着目したとき、リン酸化反応は分子Xのs1の状態がuからpから変わることである。 すなわち、s1の状態がuである分子Xのパターンに当てはまる分子種をとってきて、その状態をpに変えると考える。





In [17]:



    


with reaction_rules():
    X(s1=u) > X(s1=p) | 1
    X(s2=u) > X(s2=p) | 2
    X(s3=u) > X(s3=p) | 3

これらは個々の具体的な分子種を左右に持った反応ではなく、分子種のパターンを与えるものであるから反応則ReactionRuleと呼ぶ。

反応則はあくまでもパターンに基づいたルールでしかないため、これだけでは実際にどんな分子種があるのかはわからない。

上の例で言えば左辺にX(s1=u)という分子種があるからといってそれが実在するとは限らない。

実際に最初に想定したのはX(s1=u, s2=u, s3=u)といった3つの修飾をもった分子種であった。

ルールベース・モデルでは種seedとなる分子種を与え、そこから繰り返し反応則を適用して出来うる全ての分子種を芋づる的に生成する。ここで種となる分子種は後でWorldに追加する分子種だと思えば良い。





In [18]:



    


m = get_model(seeds=[Species('X(s1=u, s2=u, s3=u)')])
print(len(m.reaction_rules()))



    


















    






12

ここで与えたX(s1=u, s2=u, s3=u)という種から12の反応が望み通りに得られることがわかる。

このモデルを使って計算してみると、





In [19]:



    


run_simulation(5, y0={'X(s1=u,s2=u,s3=u)': 120}, model=m,
               species_list=('X', 'X(s1=p)', 'X(s2=p)', 'X(s3=p)'),
               opt_kwargs={'ylim': (0, 130)})

ワイルドカード

分子種Aに結合部位bがあったとき「誰とも結合していないA」の数が知りたければA(b)とすれば良い。

では、「誰かと結合しているA」はどのように表現すればよいのだろうか？





In [20]:



    


with reaction_rules():
    A(b) + B(b) > A(b^1).B(b^1) | 1
    A(b^_, s=u) > A(b^_, s=p) | 2



    











In [21]:



    


m = get_model()

2つ目の反応則は修飾sをuからpに変えるが、Aがbを介して誰かと結合しているときにしか起こらない.

_がワイルド・カード、すなわちどんなものとも合致する記号である。 誰かと結合しているAはA(b^_)と表現される。

修飾にワイルドカードを使うと、

A(s=_)は修飾sの状態が何であっても良いことを意味する。

何であっても良いならば書く必要ない気もするが、Aにsという修飾があることを保証する。

単位分子種名にワイルドカードを使うと、

_(s=u)と書けば状態uの修飾sを持つどんな分子種にも当てはまる。試すと、





In [22]:



    


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

run_simulation(
    5, y0={'A(s=u)': 10, 'B(s=u)': 20, 'C(b^1,s=u).C(b^1,s=u)': 15}, model=m,
    species_list=('A(s=u)', 'B(s=u)', 'C(s=u)', '_(s=u)'))

分子種の属性を決める場合にもワイルドカードが使える。





In [23]:



    


with species_attributes():
    _ | {'D': 1, 'radius': 0.005}



    











In [24]:



    


m = get_model()

動画を作成する

viz.plot_movieを使う方法

Worldを可視化するのにはviz.plot_worldが使えた。

これらをつなげて動画にするにはどうすれば良いか？まずは各時刻における状態をObserverを使って記録する。





In [25]:



    


with species_attributes():
    A | {'D': 1, 'location': 'C', 'dimension': 3}
    B | {'D': 1, 'location': 'N', 'dimension': 3}
    N | {'location': 'C', 'dimension': 3}
    C | {'dimension': 3}

with reaction_rules():
    A + N > B | 1e-4
    B + C > A | 1e-4

m = get_model()

w = spatiocyte.World(ones(), 0.005)
w.bind_to(m)
w.add_structure(Species('C'), Sphere(0.5 * ones(), 0.441))
w.add_structure(Species('N'), Sphere(0.5 * ones(), 0.350))
w.add_molecules(Species('A'), 720)



    











In [26]:



    


viz.plot_world(w, species_list=('C', 'N'))

全ての状態を保存するにはFixedIntervalHDF5Observerを使う。





In [27]:



    


sim = spatiocyte.Simulator(w)
obs1 = NumberObserver(('A', 'B'))
obs2 = FixedIntervalHDF5Observer(0.01, 'data/sample%03d.h5')
sim.run(0.5, (obs1, obs2))
viz.plot_number_observer(obs1)

あとはこのObserverをviz.plot_movieで可視化する。





In [28]:



    


mpl.rcParams.update(inline_rc)



    


















    






/home/kaizu/src/ecell4/local/lib/python3.6/site-packages/matplotlib-3.0.2-py3.6-linux-x86_64.egg/matplotlib/__init__.py:855: MatplotlibDeprecationWarning: 
examples.directory is deprecated; in the future, examples will be found relative to the 'datapath' directory.
  "found relative to the 'datapath' directory.".format(key))
/home/kaizu/src/ecell4/local/lib/python3.6/site-packages/matplotlib-3.0.2-py3.6-linux-x86_64.egg/matplotlib/__init__.py:846: MatplotlibDeprecationWarning: 
The text.latex.unicode rcparam was deprecated in Matplotlib 2.2 and will be removed in 3.1.
  "2.2", name=key, obj_type="rcparam", addendum=addendum)

















In [29]:



    


viz.plot_movie(obs2, species_list=('A', 'B'))



    


















    























In [30]:



    


import seaborn
seaborn.set(font_scale=1.5)
mpl.rc("figure", figsize=(6, 4))

ParaViewを使う方法

ParaView（http://www.paraview.org/）というソフトウェアを使うと動画作成などより詳細な可視化が行える。

  • まず前回と同様にFixedIntervalCSVObserverで状態をCSV形式で保存する。




In [31]:



    


with species_attributes():
    A | {'D': 1, 'location': 'C', 'dimension': 3}
    B | {'D': 1, 'location': 'N', 'dimension': 3}
    N | {'location': 'C', 'dimension': 3}
    C | {'dimension': 3}

with reaction_rules():
    A + N > B | 1e-4
    B + C > A | 1e-4

m = get_model()

w = spatiocyte.World(ones(), 0.005)
w.bind_to(m)
w.add_structure(Species('C'), Sphere(0.5 * ones(), 0.441))
w.add_structure(Species('N'), Sphere(0.5 * ones(), 0.350))
w.add_molecules(Species('A'), 720)



    











In [32]:



    


sim = spatiocyte.Simulator(w)
obs = FixedIntervalCSVObserver(0.01, 'data/sample%03d.csv', ('A', 'B'))
sim.run(0.5, obs)
  • File->Openで保存したCSV形式のデータを選択する。保存した連番のファイルは下図のようにまとめて表示される。読み込んだら左下にあるPropertiesでApplyボタンを押す。
  • 読み込んだ後はPipelineによってデータを可視化する。最初にCSVテーブルから座標への変換する。Filters->Alphabetical->Table To Pointsを選択し、左下からX, Y, Z Columnとしてx, y, zを選択しApplyボタンを押す。
  • Filters->Common->Glyphを選択する。同様に左下のPropertiesからGlyph TypeをSphereに、Scale Modeをscalarにして下のScale Factor横のReset using current data valuesとヒントの出るボタンを押す。今回は4を入力してApplyボタンを押す。
  • 先ほどのPropertiesにColoringという項目が増えるのでsidを選ぶ。何も見えなければ左上のPipeline BrowserでGlyph1の左の目玉がアクティブにする。
  • 上のPlayボタンを押す。赤色が増える様子を確認したらFile->Save Animation...から動画が保存できる。

ParaViewのマクロを使う方法

以上の作業をあらかじめ登録されたマクロを使って省略できる。CSVファイルをParaViewで読み込むまでは同じ。

  • Macros->Add new macro...を選び、ecell4paraview.pyからダウンロードしたecell4paraview.pyを選択する。
  • Pipeline Browserで読み込んだCSVファイルを選択してからMacros->ecell4paraviewを選択する。

Factoryの利用

Factoryクラスを用いたソルバの切り替え

  • ode、gillespie、meso、latticeのモジュールを使いわけの違いはほぼWorldを作るときだけ。
  • しかし、Simulatorの名称やWorldを作る際のオプショナルな引数もある。
  • こうした問題を解消するため、各モジュールはFactoryと呼ばれる便利クラスがある。

Factoryだけを切り替えることでモジュールが切り替えられる。





In [33]:



    


factory = ode.Factory()
# f = gillespie.Factory()
# f = meso.Factory(Integer3(4, 4, 4))
# f = spatiocyte.Factory(0.005)

w = factory.world(ones())
w.bind_to(m)
w.add_molecules(Species('C'), 60)
sim = factory.simulator(w)
obs = FixedIntervalNumberObserver(0.1, ('A', 'B', 'C'))
sim.run(10, obs)

run_simulationでもFactoryを使うことでspatiocyteのVoxel半径のような特殊な引数を扱える。





In [34]:



    


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

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

run_simulation(10, {'C': 60}, solver=('spatiocyte', 0.01))

練習問題

  • 構造体の練習: Unionを使ってRodを自前で作ってみる。
  • ルールベース・モデリングにおける分子種の練習: 2つのX分子からなる二量体はどのように表現すればよいか。また、Xの数を数えるときの扱いに関して試す。

  • ルールベース・モデリングにおける反応則の練習: 二種類の単位分子種XとYがあり、以下の反応が独立に成立するとき、何通りの分子種が存在するか？

    1. X同士は二量体を形成する。
    2. XとYは結合する。
    3. Yには1つ修飾部位があり、2種類の状態のいずれかをとる。

    ただし、対称性に注意せよ。

  • 分子の位置の記録に関する練習: 今回、FixedIntervalCSVObserverでCSV形式のファイルに記録したが、テキストエディタや表計算ソフトなどで開いてみて実際にどのようなことが書き込まれているのか確認する。
  • ParaViewによる可視化に関する練習: ParaViewは設定によって色々と見栄えが変えられるので自分好みの動画を作ってみる。

Content source: kaizu/ecell4-lectures

Similar notebooks:

notebook.community | gallery | about