E-Cell System version 4 日本語ドキュメント

E-Cell Systemは、細胞のような複雑で異種のマルチスケールシステムのモデリング、シミュレーション、解析のためのソフトウェアプラットフォームです。

E-Cell4は、GNU General Public Licenseバージョン2の下でライセンスされたフリーでオープンソースのソフトウェアです。ソースコードは GitHub で入手できます。

インストール手順 については、https://github.com/ecell/ecell4 を参照してください。

Tutorials

1. E-Cell4を用いたシミュレーションの概要

まず下記のようにプロットライブラリ(matplotlib)とE-Cell4ライブラリを読み込む必要があります

%matplotlib inline
from ecell4 import *

1.1. クイックデモ

E-Cell Systemバージョン4には Model World Simulator の3つの基本コンポーネントがあります． これらのコンポーネントはシミュレーションにおける下記の概念を記述します．

• Model はその名の通りシミュレートする内容を記述します．
• World はある状態を記述します．例えば初期条件やあるタイムポイントでの状態などです．
• Simulator はソルバーについて記述します.

Model はソルバーから独立しています。 すべてのソルバーは、単一の Model インスタンスを共有することができます。 各ソルバーのアルゴリズムには対応する World と Simulator のペアがあります（これらのペアは Factory クラスにカプセル化されています）． World は必ずしも Model と Simulator に結びついている必要はありませんが，Simulator は Model と World の両方が必要です．

シミュレーションを実行する前に Model を作成する必要があります． E-Cell4は Model を構築するための複数の方法をサポートしています． (モデルの構築方法も参照してください．) ここでは reaction_rules 関数と with ステートメントを使用する最も簡単な方法を説明します．

with reaction_rules():
    A + B > C | 0.01  # equivalent to create_binding_reaction_rule
    C > A + B | 0.3   # equivalent to create_unbinding_reaction_rule

m1 = get_model()
print(m1)

<ecell4.core.NetworkModel object at 0x7fac0c0b4fb0>

reaction_rules の後に括弧（）を書くことを覚えておいてください． ここでは、m1 という名前の2つの ReactionRule を持つモデルが作成されました． with ブロックの行は、それぞれ結合と乖離の ReactionRule を表しています． 質量作用の法則を示す反応の速度はセパレーター | の後に定義します． すなわち 0.01 もしくは 0.3 が反応速度になります． 通常の微分方程式の形では，このモデルは次のように記述できます．

よりコンパクトな記述として A + B == C | (0.01, 0.3) と書くことも可能です．

E-Cell4には、run_simulation という特定のモデルでシミュレーションを実行するためのシンプルなインターフェイスがあります。 これにより、自分で World とSimulator をインスタンス化せずにシミュレーションを実行することができます。 (ただモデルを解くには、ボリューム、各種の初期値と時間は与えなければなりません。)

run_simulation(10.0, model=m1, y0={'A': 60, 'B': 60}, volume=1.0)

シミュレーションアルゴリズムを切り替えるのは、次のようにソルバの型を指定するだけでできます。 ( デフォルトでは ode が使われるようになっています)

run_simulation(10.0, model=m1, y0={'A': 60, 'B': 60}, solver='gillespie')

1.2. 空間シミュレーションとその可視化

E-Cell4は、複数の空間シミュレーションアルゴリズム、egfrd、spatiocyte、meso をサポートするようになりました。 これらの空間ソルバーは、非空間ソルバー ( ode および gillespie ) で使用されるモデルに加えて、各 Species についての追加の情報 (拡散係数および半径) を必要とします。

species_attributes を持つ with ステートメントは、これらのプロパティ情報を記述するために利用できます：

with species_attributes():
    A | B | C | {'radius': '0.005', 'D': '1'}  # 'D' is for the diffusion coefficient

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

m2 = get_model()

プロパティーは浮動小数点数を表しますが、各属性は文字列として指定する必要があります。

今度は、上記と同じ方法で空間シミュレーションを実行できます (egfrd の場合、シミュレーションし終わるまでにかかる時間は長くなります)。

run_simulation(10.0, model=m2, y0={'A': 60, 'B': 60}, solver='meso')

Structure (例えば 膜、細胞質および核といったもの) は、今のところ spatiocyte および meso でのみのサポートとなっています。 シミュレーションのための各 Species が属する location は、その属性で指定しなければなりません。

with species_attributes():
    A | {'D': '1', 'location': 'S'}  # 'S' is a name of the structure

m3 = get_model()  # with no reactions

E-Cell4はプリミティブなシェイプを Sphere のような構造体としてサポートしています：

sphere = Sphere(Real3(0.5, 0.5, 0.5), 0.48)  # a center position and radius

E-Cell4は、シミュレーション中に状態を記録する様々な種類の Observer を提供しています。 下記では2つの Observer が分子の位置を記録すると宣言されています。 FixedIntervalTrajectoryObserver は分子の軌跡を記録し、FixedIntervalHDF5Observer は World を与えられた時間間隔でHDF5ファイルに保存します：

obs1 = FixedIntervalTrajectoryObserver(1e-3)
obs2 = FixedIntervalHDF5Observer(0.1, 'test%02d.h5')

run_simulation は構造体とオブザーバを引数 structures と observers として受け取りることができます（help（run_simulation）も参照してください）：

run_simulation(1.0, model=m3, y0={'A': 60}, structures={'S': sphere},
               solver='spatiocyte', observers=(obs1, obs2), return_type=None)

E-Cell4には、viz.plot_world というインタラクティブな視覚化が可能な World を視覚化する機能があります。 viz.plot_world は、分子の位置を3Dでプロットします。 さらに、load_world を使用することで、HDF5ファイルからWorldの状態を簡単に復元することができます。

# viz.plot_world(load_world('test00.h5'), species_list=['A'])
viz.plot_world(load_world('test00.h5'), species_list=['A'], interactive=True)

また、 FixedIntervalTrajectoryObserver の場合、viz.plot_trajectory は軌道のプロットを生成します。(ここでもインタラクティブなプロットの生成が可能です。)：

# viz.plot_trajectory(obs1)
viz.plot_trajectory(obs1, interactive=True)

より詳しくは 5. How to Log and Visualize Simulations を参照してください。

2. モデルを構築する方法

ModelはSpecies と ReactionRule のセットで構成されます。

• Species はモデル中の分子エンティティ（例えば、タンパク質のタイプまたは状態）を表します。また Species にはサイズなどの属性もあります。
• ReactionRule は Species 間のインタラクションを表します。(例えば結合や乖離など。)

%matplotlib inline
from ecell4 import *

2.1. Species

Speciesはその名前を与えることで生成することができます

sp1 = Species("A")
print(sp1.serial())

Speciesの名前にはいくつかの命名規則があります。 この命名規約では特殊記号（例：括弧 （）、ドット .、アンダーバー_）、数字と空白の使用に注意する必要があります。 Speciesには、その属性を扱うための一連のAPIがあります。

sp1.set_attribute("radius", "0.005")
sp1.set_attribute("D", "1")
sp1.set_attribute("location", "cytoplasm")
print(sp1.get_attribute("radius"))
sp1.remove_attribute("radius")
print(sp1.has_attribute("radius"))

set_attributeの引数は、属性の名前とその値です。 どちらも文字列でなければなりません。 radius、D（拡散係数）、locationは頻繁に使われるので、一度に設定するためのショートカットがあります。

sp1 = Species("A", "0.005", "1", "cytoplasm")  # serial, radius, D, location

Species が同じものかどうかはそれらの serial に基づいて評価されます。

print(Species("A") == Species("B"), Species("A") == Species("A"))

False True

Speciesは、1つ以上の UnitSpeciesから成ります。

sp1 = Species()
usp1 = UnitSpecies("C")
print(usp1.serial())
sp1.add_unit(usp1)
sp1.add_unit(UnitSpecies("A"))
sp1.add_unit(UnitSpecies("B"))
print(sp1.serial(), sp1.num_units())

C
C.A.B 3

Speciesは、serialから複製することができます。 serialでは、UnitSpeciesのすべての連番が区切り文字「。」で結合されます。 UnitSpeciesの順序はSpeciesの比較時に影響してきます。

sp1 = Species("C.A.B")
print(sp1.serial())
print(Species("A.B.C") == Species("C.A.B"))
print(Species("A.B.C") == Species("A.B.C"))

C.A.B
False
True

UnitSpeciesはサイトを持つことができます。 サイトは name、state、 bondからなり、UnitSpeciesで自動的にソートされます。 nameはUnitSpeciesで一意でなければなりません。 すべての値は文字列でなければなりません。 括弧、ドット、空白は含めず、 bond以外の数字から始めることはできません。

usp1 = UnitSpecies("A")
usp1.add_site("us", "u", "")
usp1.add_site("ps", "p", "_")
usp1.add_site("bs", "", "_")
print(usp1.serial())

A(bs^_,ps=p^_,us=u)

UnitSpecies can be also reproduced from its serial. Please be careful with the order of sites where a site with a state must be placed after sites with no state specification:

usp1 = UnitSpecies()
usp1.deserialize("A(bs^_, us=u, ps=p^_)")
print(usp1.serial())

A(bs^_,ps=p^_,us=u)

Of course, a site of UnitSpecies is available even in Species』 serial.

sp1 = Species("A(bs^1, ps=u).A(bs, ps=p^1)")
print(sp1.serial())
print(sp1.num_units())

A(bs^1,ps=u).A(bs,ps=p^1)
2

The information (UnitSpecies and its site) is used for rule-based modeling. The way of rule-based modeling in E-Cell4 will be explained in 7. Introduction of Rule-based Modeling local ipynb readthedocs.

2.2. ReactionRule

ReactionRule は reactants、 productsとkで構成されます。 reactantsとproductsは Speciesのリストで、kは浮動小数点として与えられた反応速度です。

rr1 = ReactionRule()
rr1.add_reactant(Species("A"))
rr1.add_reactant(Species("B"))
rr1.add_product(Species("C"))
rr1.set_k(1.0)

これは、AとBからCへの結合反応になります。 この反応定義では、Speciesに属性を設定する必要はありません。 上記の一連の操作はユーティリティ関数を使うと一行で記述できます。 ReactionRuleを検査するためにはas_string関数が利用できます：

rr1 = create_binding_reaction_rule(Species("A"), Species("B"), Species("C"), 1.0)
print(rr1.as_string())

A+B>C|1

You can also provide components to the constructor:

rr1 = ReactionRule([Species("A"), Species("B")], [Species("C")], 1.0)
print(rr1.as_string())

A+B>C|1

基本的にReactionRuleはレートkを用いた質量作用反応を表します。

ode solver also supports rate laws thought it’s under development yet.

ode.ODERatelaw is explained in 6. How to Solve ODEs with Rate Law Functions local ipynb readthedocs.

2.3. NetworkModel

いくつかのModelコンポーネントを作成する方法を学びました。 次に、これらのコンポーネントをModelに登録してみましょう。

sp1 = Species("A", "0.005", "1")
sp2 = Species("B", "0.005", "1")
sp3 = Species("C", "0.01", "0.5")

rr1 = create_binding_reaction_rule(Species("A"), Species("B"), Species("C"), 0.01)
rr2 = create_unbinding_reaction_rule(Species("C"), Species("A"), Species("B"), 0.3)

SpeciesとReactionRuleはadd_species_attributeとadd_reaction_ruleで登録できます。

m1 = NetworkModel()
m1.add_species_attribute(sp1)
m1.add_species_attribute(sp2)
m1.add_species_attribute(sp3)
m1.add_reaction_rule(rr1)
m1.add_reaction_rule(rr2)

これで結合と乖離の反応を持つ単純なモデルができました。 Modelの内容を確認するには species_attributesとreaction_rulesが利用できます。

print([sp.serial() for sp in m1.species_attributes()])
print([rr.as_string() for rr in m1.reaction_rules()])

['A', 'B', 'C']
['A+B>C|0.01', 'C>A+B|0.3']

モデルの種の属性は、空間シミュレーションには不可欠ですが、非空間アルゴリズム（gillespieとode）では必ずしも必要ではありません。 最初にプッシュされた属性は、後でプッシュされる属性よりも優先度が高くなります。 また、モデルの属性に基づいて種を属性付けすることもできます。

Species』 attributes in Model are indispensable for spatial simulations, but not necessarily needed for non-spatial algorithms, i.e. gillespie and ode. The attribute pushed first has higher priority than one pushed later. You can also attribute a Species based on the attributes in a Model.

sp1 = Species("A")
print(sp1.has_attribute("radius"))
sp2 = m1.apply_species_attributes(sp1)
print(sp2.has_attribute("radius"))
print(sp2.get_attribute("radius"))

False
True
0.005

Species, ReactionRule, NetworkModel に関連するすべての関数は、C++でも同じように利用できます。

run_simulationでこのモデルを解くことができます。1.「E-Cell4シミュレーションの概要」 You can solve this model with run_simulation as explained in 1. Brief Tour of E-Cell4 Simulations local ipynb readthedocs :

run_simulation(10.0, model=m1, y0={'C': 60})

2.4. モデル構築のためのPythonユーティリティ

1. E-Cell4を用いたシミュレーションの概要 local ipynb readthedocs で示したように, E-Cell4 は with ステートメントを用いた簡単なモデルの構築方法を提供しています。

with species_attributes():
    A | B | {'radius': '0.005', 'D': '1'}
    C | {'radius': '0.01', 'D': '0.5'}

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

m1 = get_model()

可逆反応の場合、 <>はPython 2では ==の代わりに利用できますが、Python3では非推奨です。 withステートメントでは、宣言されていない変数は自動的に Speciesとみなされます。 Pythonの変数、関数、および文は、 withブロックでも使用できます。

from math import log

ka, kd, kf = 0.01, 0.3, 0.1
tau = 10.0

with reaction_rules():
    E0 + S == ES | (ka, kd)

    if tau > 0:
        ES > E1 + P | kf
        E1 > E0 | log(2) / tau
    else:
        ES > E0 + P | kf

m1 = get_model()
del ka, kd, kf, tau

一方、ブロックの外側でも変数が宣言されると、以下のようにその変数を Speciesとして使用することはできません：

A = 10

try:
    with reaction_rules():
        A + B == C | (0.01, 0.3)
except Exception as e:
    print(repr(e))

del A

RuntimeError('invalid expression; "10" given',)

where A + B == C exactly means 10 + B == C.

合成や分解のような ReactionRuleの左辺や右辺がない場合は、次のように記述したいかもしれません。

with reaction_rules():
    A > | 1.0  # XXX: will raise SyntaxError
    > A | 1.0  # XXX: will raise SyntaxError

しかし、これはPythonの SyntaxErrorのために受け入れられません。

このケースを記述するために、特別な演算子、チルダ 〜が利用可能です。 〜は、以下の Speciesの化学量論係数をゼロに設定します。 つまり、Speciesは ReactionRuleで無視されます。

with reaction_rules():
    ~A > A | 2.0  # equivalent to `create_synthesis_reaction_rule`
    A > ~A | 1.0  # equivalent to `create_degradation_reaction_rule`

m1 = get_model()
print([rr.as_string() for rr in m1.reaction_rules()])

['>A|2', 'A>|1']

The following Species』 name is not necessarily needed to be the same as others. The model above describes :

from math import exp
run_simulation(10.0, model=m1, opt_args=['-', lambda t: 2.0 * (1 - exp(-t)), '--'])

一連の反応は1行で記述することができます。 行を2つ以上の物理行に分割するには、括弧で囲みます。

with reaction_rules():
    (E + S == ES | (0.5, 1.0)
         > E + P | 1.5)

m1 = get_model()
print([rr.as_string() for rr in m1.reaction_rules()])

['E+S>ES|0.5', 'ES>E+S|1', 'ES>E+P|1.5']

The method uses global variables in ecell4.util.decorator (e.g. REACTION_RULES) to cache objects created in the with statement:

import ecell4.util.decorator

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

print(ecell4.util.decorator.REACTION_RULES)  #XXX: Only for debugging
get_model()
print(ecell4.util.decorator.REACTION_RULES)  #XXX: Only for debugging

[<ecell4.core.ReactionRule object at 0x000001A123D02BA0>, <ecell4.core.ReactionRule object at 0x000001A123D02CD8>]
[]

Modelを構築する際のモジュール性のために、デコレータ関数も便利です。 Pythonのデコレータ関数も利用できます。 デコレータ関数は、モデルのモジュール性を高めます。

@species_attributes
def attrgen1(radius, D):
    A | B | {'radius': str(radius), 'D': str(D)}
    C | {'radius': str(radius * 2), 'D': str(D * 0.5)}

@reaction_rules
def rrgen1(kon, koff):
    A + B == C | (kon, koff)

attrs1 = attrgen1(0.005, 1)
rrs1 = rrgen1(0.01, 0.3)
print(attrs1)
print(rrs1)

[<ecell4.core.Species object at 0x000001A123D02CD8>, <ecell4.core.Species object at 0x000001A123D02BA0>, <ecell4.core.Species object at 0x000001A123D02DC8>]
[<ecell4.core.ReactionRule object at 0x000001A123D02D20>, <ecell4.core.ReactionRule object at 0x000001A123D02DE0>]

withステートメントの場合とは対照的で、ここではデコレータの後に括弧を追加しないでください。 reaction_rulesとspecies_attributesでデコレートされた関数はそれぞれReactionRulesとSpeciesのリストを返します。 add_reaction_rulesとadd_species_attributesを使用して、リストをModelに一度に登録できます。

m1 = NetworkModel()
m1.add_species_attributes(attrs1)
m1.add_reaction_rules(rrs1)
print(m1.num_reaction_rules())

2

This method is modular and reusable relative to the way using with statement.

!pip install ecell

Requirement already satisfied (use --upgrade to upgrade): ecell in /home/nbcommon/anaconda3_410/lib/python3.5/site-packages
You are using pip version 8.1.2, however version 9.0.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.

3. 初期条件を設定する方法

ここでは、World クラスの基礎を説明します。 E-Cell4では、spatiocyte.SpatiocyteWorld, egfrd.EGFRDWorld, bd.BDWorld, meso.MesoscopicWorld, gillespie.GillespieWorld, そして ode.ODEWorldの6種類のWorldクラスがサポートされています。

ほとんどのソフトウェアでは、初期状態は Modelの一部であると考えられています。 しかし、E-Cell4では、初期条件は Modelとは別にWorldとして設定する必要があります。 World には、現在時刻、分子数、分子座標、構造、乱数発生器などの状態に関する情報が格納されています。 一方、Model には分子間の相互作用のタイプと分子の共通の性質が含まれています。

import ecell4

3.1. Worldの共通API

Worldは、対応するアルゴリズムに固有の空間表現を表していますが、互換性のあるAPIを持っています。 このセクションでは、6つの Worldクラスの共通インターフェースを紹介します。

from ecell4.core import *
from ecell4.gillespie import GillespieWorld
from ecell4.ode import ODEWorld
from ecell4.spatiocyte import SpatiocyteWorld
from ecell4.bd import BDWorld
from ecell4.meso import MesoscopicWorld
from ecell4.egfrd import EGFRDWorld

Worldクラスは、アルゴリズムに固有のパラメータを決定するコンストラクタ内に異なるセットの引数を受け取ります。 しかし、少なくとも、すべてのWorldクラスは、edge_lengthsという名前のサイズでしかインスタンス化できません。 edge_lengthsのタイプはReal3であり、これはRealsのトリプレットを表します。 E-Cell4では、すべての3次元位置がReal3として扱われます。

edge_lengths = Real3(1, 2, 3)
w1 = GillespieWorld(edge_lengths)
w2 = ODEWorld(edge_lengths)
w3 = SpatiocyteWorld(edge_lengths)
w4 = BDWorld(edge_lengths)
w5 = MesoscopicWorld(edge_lengths)
w6 = EGFRDWorld(edge_lengths)

Worldにはサイズと体積のゲッターメソッドがあります。

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.0, 2.0, 3.0) 6.0
(1.0, 2.0, 3.0) 6.0
(1.0, 2.0, 3.0) 6.0
(1.0, 2.0, 3.0) 6.0

次に、分子をWorldに加えましょう。 ここでは、分子の形状を伝えるために、「半径」と「拡散係数」の属性付きのSpeciesを与える必要があります。 下記の例では「0.0025」が半径に，「1」が拡散係数に相当します。 分子の位置は、必要に応じてWorldでランダムに決定されます。 add_molecules関数における10は追加する分子の数を表しています。

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)

ModelがWorldにバインドされた後、一度Speciesで設定した半径と拡散係数を再び書く必要はありません（変更したくない限り）。

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)

同様に、remove_moleculesとnum_molecules_exactも利用できます。

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

num_molecules_exactと異なりnum_moleculesは与えられたSpeciesと一致するすべての数をルールベースの方法で返します。 WorldのすべてのSpeciesが分子間相互作用を持たない場合、num_moleculesは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 はシミュレーションの時間を保持しています。

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

最後に、HDF5ファイルにWorldの状態を保存したりロードすることができます。

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 = GillespieWorld()
w2 = ODEWorld()
w3 = SpatiocyteWorld()
w4 = BDWorld()
w5 = MesoscopicWorld()
w6 = EGFRDWorld()
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.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
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.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
1.0 (1.0, 2.0, 3.0) 6.0 10 15

すべてのWorldクラスでは、コンストラクターの唯一の引数としてHDF5ファイルパスも使用できます。

print(GillespieWorld("gillespie.h5").t())
print(ODEWorld("ode.h5").t())
print(SpatiocyteWorld("spatiocyte.h5").t())
print(BDWorld("bd.h5").t())
print(MesoscopicWorld("meso.h5").t())
print(EGFRDWorld("egfrd.h5").t())

1.0
1.0
1.0
1.0
1.0
1.0

3.2. 分子の位置を取得する方法

World には、分子の座標にアクセスする共通の機能もあります。

w1 = GillespieWorld()
w2 = ODEWorld()
w3 = SpatiocyteWorld()
w4 = BDWorld()
w5 = MesoscopicWorld()
w6 = EGFRDWorld()

まず、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, p3), suc3 = 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は作成されたパーティクルとそれが成功したかどうかを返します。 分子の表現における分解能は異なっています。 たとえば、GillespieWorld には分子の座標に関する情報はほとんどありません。 それ故 GillespieWorld は与えられた位置を無視し、分子の数をカウントアップすることのみ行います。

ParticleID は lot と serial と名付けられた Integer のペアになります。

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

0 1
True

パーティクルシミュレータ、すなわちspatiocyte、bdおよびegfrdは、idでパーティクルにアクセスするためのインターフェイスを提供しています。 has_particle は、指定された 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

存在するかをチェックした後で、次のように get_particle でパーティクルを取得できます。

# 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 は species, position, radius そして 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

Spatiocyteの場合、粒子の位置は与えられた位置に最も近いボクセルの中心に自動的に丸められます。

パーティクルの位置を移動させることもできます。 update_particle は指定された ParticleID で指定されたパーティクルを指定された Particle に置き換え、Falseを返します。 ParticleID に対応するパーティクルが見つからない場合は、新しいパーティクルを作成してTrueを返します。 異なるタイプの Species でパーティクルを与えると、パーティクルの Species も変更されます。

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 と list_particles_exact は、World の ParticleID と Particle のペアのリストを返します。

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.core.ParticleID object at 0x7fef945e3e40>, <ecell4.core.Particle object at 0x7fef945e3ea0>)]
[(<ecell4.core.ParticleID object at 0x7fef945e3e40>, <ecell4.core.Particle object at 0x7fef945e3e88>)]
[(<ecell4.core.ParticleID object at 0x7fef945e3e40>, <ecell4.core.Particle object at 0x7fef945e3e70>)]
[(<ecell4.core.ParticleID object at 0x7fef945e3e40>, <ecell4.core.Particle object at 0x7fef945e3ea0>)]
[(<ecell4.core.ParticleID object at 0x7fef945e3e40>, <ecell4.core.Particle object at 0x7fef945e3e88>)]

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 シミュレータでの座標情報

共通のインターフェースに加えて、各Worldは独自のインターフェースを持つことができます。 ここでは、Lattice シミュレータで格子ベースの座標を扱う方法を例として説明します。 SpatiocyteWorldは、六方最密充填格子で構成される LatticeSpace に離散化された空間に基づいています。

w = SpatiocyteWorld(Real3(1, 2, 3), voxel_radius=0.01)
w.bind_to(m)

Voxel と呼ばれる単一の格子のサイズは、voxel_radius メソッドによって取得できます。 SpatiocyteWorld には、行数、列数、およびレイヤ数を取得するメソッドがあります。 これらのサイズは、World構築時に与えられた edge_lengths に基づいて自動的に計算されます。

print(w.voxel_radius())  # => 0.01
print(tuple(w.shape()))  # => (62, 152, 116)
# print(w.col_size(), w.row_size(), w.layer_size())  # => (62, 152, 116)
print(w.size())  # => 1093184 = 62 * 152 * 116

0.01
(64, 154, 118)
1163008

格子ベースの空間内の位置は、グローバル座標と呼ばれるInteger3 (列、行および層)として扱われます。 このようにして、SpatiocyteWorld は、Real3を格子ベースの座標に変換する機能を提供します。

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

SpatiocyteWorldでは、グローバル座標は単一の整数に変換されます。 それは単に座標と呼ばれています。 グローバル座標で同じ方法で座標を扱うこともできます。

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

これらの座標を使用すると、Particleオブジェクトを表す Voxel を処理できます。 new_particleの代わりに、new_voxelは新しいボクセルを座標で作成する方法を提供します。

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)

<ecell4.core.ParticleID object at 0x7fef945e3ee8> <ecell4.core.Voxel object at 0x7fef945e3f18> True

Voxel は、species、座標、半径およびD からなります。

print(v.species().serial(), v.coordinate(), v.radius(), v.D())  # => (u'A', 278033, 0.0025, 1.0)

A 300634 0.0025 1.0

もちろん、get_particleとlist_particles_exact同様にget_voxelとlist_voxels_exactを使ってVoxelとVoxelのリストを取得することができます。

print(w.num_voxels_exact(Species("A")))
print(w.list_voxels_exact(Species("A")))
print(w.get_voxel(pid))

1
[(<ecell4.core.ParticleID object at 0x7fef945e3f60>, <ecell4.core.Voxel object at 0x7fef945e3f90>)]
(<ecell4.core.ParticleID object at 0x7fef945e3f60>, <ecell4.core.Voxel object at 0x7fef945e3f48>)

update_particleに対応するupdate_voxelでボクセルを移動および更新することができます。

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

300659
A 300659 0.0025 1.0
1

最後に、remove_voxelは、remove_particleのようにVoxelを削除します。

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

True
False

3.4. Structure

w1 = GillespieWorld()
w2 = ODEWorld()
w3 = SpatiocyteWorld()
w4 = BDWorld()
w5 = MesoscopicWorld()
w6 = EGFRDWorld()

Shapeオブジェクトを使用すると、分子の初期位置をWorldの一部に限定できます。 以下の場合、60個の分子が与えられた球の内部に配置されます。 ここに置かれた分子の拡散はShapeに制限されていません。 このShapeは、初期化専用です。

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)

分子の拡散を制限するためにはSpeciesのプロパティ、「location」が利用可能です。 locationはspatiocyteとmesoでのみサポートされています。 add_structureは、SpeciesとShapeのペアとして与えられる新しいStructureを定義します。

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)

Structureを定義した後、次のように分子をStructureにバインドできます。

sp2 = Species("B", "0.0025", "0.1", "M")  # `'location'` is the fourth argument
w3.add_molecules(sp2, 60)
w5.add_molecules(sp2, 60)

Bという名前のSpeciesにバインドした分子sp2は、SphericalSurface（中空の球）の形をしたMという名前の構造上に拡散します。 Spatiocyteでは、Structureは、Species Mが Voxelを占める粒子の集合として表されます。 これは、Structureに属さない分子がVoxelと重なり合うことができず、衝突を引き起こすことを意味します。

一方、mesoでは、Structureとはサブボリュームのリストを意味します。 したがって、構造は他の粒子の侵入を避けることはできません。

3.5. Random Number Generator

Random Number GeneratorもWorldの一部です。 ODEWorld以外のすべてのWorldでは、Random Number Generatorが格納されており、シミュレーションにランダムな値が必要なときに更新されます。 E-Cell4では、Random Number GeneratorとしてGSLRandomNumberGeneratorクラスが1つだけ実装されています。

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]

引数を指定しない場合、Random Number Generatorは常に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]

次のように、シードを整数で初期化することができます。

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]

入力なしでシード関数を呼び出すと、シードは現在時刻の情報から生み出されます。

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

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

GSLRandomNumberGeneratorは、乱数を得るためのいくつかの方法を提供します。

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

0.03033520421013236
33
0.8935555455208181

Worldは、構築時にRandom Number Generatorを受け入れます。 デフォルトでは、GSLRandomNumberGenerator()が使用されます。 そのため、Generatorを与えないと、シミュレーションの振る舞いは常に同じです（これをdeterminisiticと呼びます）。

rng = GSLRandomNumberGenerator()
rng.seed()
w1 = GillespieWorld(Real3(1, 1, 1), rng=rng)

World中のGSLRandomNumberGeneratorへはrng関数を使ってアクセスできます。

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

0.37904000049456954

rng()は、GSLRandomNumberGeneratorへの共有ポインタを返します。 したがって、上記の例では、rngとw1.rng()は全く同じことを指しています。

!pip install ecell

Collecting ecell
  Downloading ecell-4.1.2-cp36-cp36m-manylinux1_x86_64.whl (45.5MB)
    100% |################################| 45.5MB 29kB/s  eta 0:00:011 day, 23:59:59    68% |#####################           | 31.0MB 52.8MB/s eta 0:00:01    91% |#############################   | 41.9MB 51.9MB/s eta 0:00:01    98% |############################### | 44.8MB 59.9MB/s eta 0:00:01
Installing collected packages: ecell
Successfully installed ecell-4.1.2

4. シミュレーションの走らせ方

セクション2と3では、モデルを構築し、初期状態を設定する方法について説明しました。 それでは、シミュレーションを実行してみましょう。 「Spatiocyte.SpatiocyteSimulator」、「egfrd.EGFRDSimulator」、「bd.BDSimulator」、「meso.MesoscopicSimulator」、「gillespie.GillespieSimulator」、および「ode.ODESimulator」の6つの「シミュレータ」クラスが存在します。 それぞれの SimulatorクラスはWorldの対応する型だけを受け入れますが、それらのすべてが同じ Modelを受け入れます。

import ecell4

4.1. シミュレーターのセットアップの方法

アルゴリズム固有の引数を持つ初期化（いわゆるコンストラクタ関数）を除いて、すべてのシミュレータは同じAPIを持っています。

from ecell4.core import *
from ecell4.gillespie import GillespieWorld, GillespieSimulator
from ecell4.ode import ODEWorld, ODESimulator
from ecell4.spatiocyte import SpatiocyteWorld, SpatiocyteSimulator
from ecell4.bd import BDWorld, BDSimulator
from ecell4.meso import MesoscopicWorld, MesoscopicSimulator
from ecell4.egfrd import EGFRDWorld, EGFRDSimulator

Simulatorを構築する前に、Simulatorのタイプに対応するModelとWorldを用意してください。

from ecell4 import species_attributes, reaction_rules, get_model

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

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

m = get_model()

w1 = GillespieWorld()
w2 = ODEWorld()
w3 = SpatiocyteWorld()
w4 = BDWorld()
w5 = MesoscopicWorld()
w6 = EGFRDWorld()

シミュレータは、構築時に下記の順序でモデルとワールドの両方を必要とします。

sim1 = GillespieSimulator(m, w1)
sim2 = ODESimulator(m, w2)
sim3 = SpatiocyteSimulator(m, w3)
sim4 = BDSimulator(m, w4)
sim5 = MesoscopicSimulator(m, w5)
sim6 = EGFRDSimulator(m, w6)

ModelをWorldにバインドした場合は、シミュレータを作成するために要るのはWorldだけです。

w1.bind_to(m)
w2.bind_to(m)
w3.bind_to(m)
w4.bind_to(m)
w5.bind_to(m)
w6.bind_to(m)

sim1 = GillespieSimulator(w1)
sim2 = ODESimulator(w2)
sim3 = SpatiocyteSimulator(w3)
sim4 = BDSimulator(w4)
sim5 = MesoscopicSimulator(w5)
sim6 = EGFRDSimulator(w6)

もちろん、SimulatorにバインドされたModelとWorldは、Simulatorから以下のように引き出すことができます。

print(sim1.model(), sim1.world())
print(sim2.model(), sim2.world())
print(sim3.model(), sim3.world())
print(sim4.model(), sim4.world())
print(sim5.model(), sim5.world())
print(sim6.model(), sim6.world())

(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.gillespie.GillespieWorld object at 0x7f6b329d1a08>)
(<ecell4.ode.ODENetworkModel object at 0x7f6b329d1af8>, <ecell4.ode.ODEWorld object at 0x7f6b329d1ae0>)
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.spatiocyte.SpatiocyteWorld object at 0x7f6b329d1a08>)
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.bd.BDWorld object at 0x7f6b329d1ae0>)
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.meso.MesoscopicWorld object at 0x7f6b329d1ac8>)
(<ecell4.core.Model object at 0x7f6b329d1af8>, <ecell4.egfrd.EGFRDWorld object at 0x7f6b329d1ae0>)

Worldを自分で更新した後は、シミュレーションを実行する前にSimulatorの内部状態を初期化する必要があります。

w1.add_molecules(Species('C'), 60)
w2.add_molecules(Species('C'), 60)
w3.add_molecules(Species('C'), 60)
w4.add_molecules(Species('C'), 60)
w5.add_molecules(Species('C'), 60)
w6.add_molecules(Species('C'), 60)

sim1.initialize()
sim2.initialize()
sim3.initialize()
sim4.initialize()
sim5.initialize()
sim6.initialize()

一定のステップ間隔を有するアルゴリズムの場合は上記に加えdtも必要です。

sim2.set_dt(1e-6)  # ODESimulator. This is optional
sim4.set_dt(1e-6)  # BDSimulator

4.2. シミュレーションの実行

シミュレーションを実行するために、Simulatorはstepとrunの2つのAPIを提供します。 step()は、シミュレータが期待する時間(next_time())の間シミュレーションを進めます。

print(sim1.t(), sim1.next_time(), sim1.dt())
print(sim2.t(), sim2.next_time(), sim2.dt())  # => (0.0, 1e-6, 1e-6)
print(sim3.t(), sim3.next_time(), sim3.dt())
print(sim4.t(), sim4.next_time(), sim4.dt())  # => (0.0, 1e-6, 1e-6)
print(sim5.t(), sim5.next_time(), sim5.dt())
print(sim6.t(), sim6.next_time(), sim6.dt())  # => (0.0, 0.0, 0.0)

(0.0, 0.032171880483933143, 0.032171880483933143)
(0.0, 1e-06, 1e-06)
(0.0, 1.6666666666666667e-05, 1.6666666666666667e-05)
(0.0, 1e-06, 1e-06)
(0.0, 0.029999537310460653, 0.029999537310460653)
(0.0, 0.0, 0.0)

sim1.step()
sim2.step()
sim3.step()
sim4.step()
sim5.step()
sim6.step()

print(sim1.t())
print(sim2.t())  # => 1e-6
print(sim3.t())
print(sim4.t())  # => 1e-6
print(sim5.t())
print(sim6.t())  # => 0.0

0.0321718804839
1e-06
1.66666666667e-05
1e-06
0.0299995373105
0.0

last_reactions()は、最後のステップで発生するReactionRuleとReactionInfoのペアのリストを返します。 各アルゴリズムにはReactionInfoの独自の実装があります。 詳細については、help(module.ReactionInfo)を参照してください。

print(sim1.last_reactions())
# print(sim2.last_reactions())
print(sim3.last_reactions())
print(sim4.last_reactions())
print(sim5.last_reactions())
print(sim6.last_reactions())

[(<ecell4.core.ReactionRule object at 0x7f6b329d1ae0>, <ecell4.gillespie.ReactionInfo object at 0x7f6b329d1af8>)]
[]
[]
[(<ecell4.core.ReactionRule object at 0x7f6b329d1ae0>, <ecell4.meso.ReactionInfo object at 0x7f6b329d1af8>)]
[]

step(upto)は、next_timeがuptoより小さい場合はnext_timeのシミュレーションを進め、それ以外の場合はuptoまで実行します。 またstep(upto)は時間が上限(upto)に達していないかどうかを返します。

print(sim1.step(1.0), sim1.t())
print(sim2.step(1.0), sim2.t())
print(sim3.step(1.0), sim3.t())
print(sim4.step(1.0), sim4.t())
print(sim5.step(1.0), sim5.t())
print(sim6.step(1.0), sim6.t())

(True, 0.11785771174833615)
(True, 2e-06)
(True, 3.3333333333333335e-05)
(True, 2e-06)
(True, 0.0709074540101994)
(True, 0.0)

upto の時間ちょうどまでシミュレーションをは走らせるには, それがTrueを返す間 step(upto) を呼びます。

while sim1.step(1.0): pass
while sim2.step(0.001): pass
while sim3.step(0.001): pass
while sim4.step(0.001): pass
while sim5.step(1.0): pass
while sim6.step(0.001): pass

print(sim1.t())  # => 1.0
print(sim2.t())  # => 0.001
print(sim3.t())  # => 0.001
print(sim4.t())  # => 0.001
print(sim5.t())  # => 1.0
print(sim6.t())  # => 0.001

1.0
0.001
0.001
0.001
1.0
0.001

これはちょうど run の動作と同じものです。 run(tau)はt()+tauの時間までシミュレーションを進めます。

sim1.run(1.0)
sim2.run(0.001)
sim3.run(0.001)
sim4.run(0.001)
sim5.run(1.0)
sim6.run(0.001)

print(sim1.t())  # => 2.0
print(sim2.t())  # => 0.002
print(sim3.t())  # => 0.002
print(sim4.t())  # => 0.002
print(sim5.t())  # => 2.0
print(sim6.t())  # => 0.02

2.0
0.002
0.002
0.002
2.0
0.002

num_stepsは、シミュレーション中のステップ数を返します。

print(sim1.num_steps())
print(sim2.num_steps())
print(sim3.num_steps())
print(sim4.num_steps())
print(sim5.num_steps())
print(sim6.num_steps())

37
2001
120
2001
34
581

4.3. アルゴリズムをカプセル化してファクトリクラスに

Modelの移植性とWorldと Simulatorの一貫したAPIのおかげで、アルゴリズムに共通するスクリプトを書くのは非常に簡単です。 しかし、アルゴリズムを切り替えるときには、コード内のクラスの名前を1つずつ書き直す必要があります。

この問題を回避するために、E-Cell4はアルゴリズムごとに Factoryクラスも提供しています。 FactoryはWorldと Simulatorを構造化に必要な引数でカプセル化します。 Factoryクラスを使うことで、あなたのスクリプトはアルゴリズムの変更に対してポータブルかつロバストになります。

from ecell4.gillespie import GillespieFactory
from ecell4.ode import ODEFactory
from ecell4.spatiocyte import SpatiocyteFactory
from ecell4.bd import BDFactory
from ecell4.meso import MesoscopicFactory
from ecell4.egfrd import EGFRDFactory

Factoryはcreate_worldと create_simulatorの2つの関数を提供します。

def singlerun(f, m):
    w = f.create_world(Real3(1, 1, 1))
    w.bind_to(m)
    w.add_molecules(Species('C'), 60)

    sim = f.create_simulator(w)
    sim.run(0.01)
    print(sim.t(), w.num_molecules(Species('C')))

上記のsinglerunはアルゴリズムから独立しています。 したがって、単にFactoryを切り替えるだけで、結果を簡単に比較することができます。

singlerun(GillespieFactory(), m)
singlerun(ODEFactory(), m)
singlerun(SpatiocyteFactory(), m)
singlerun(BDFactory(bd_dt_factor=1), m)
singlerun(MesoscopicFactory(), m)
singlerun(EGFRDFactory(), m)

(0.01, 60)
(0.01, 59)
(0.01, 60)
(0.01, 60)
(0.01, 60)
(0.01, 60)

WorldやSimulatorを初期化するためにいくつかのパラメータを指定する必要があるとき、 run_simulationはsolverの代わりに Factoryを受け取ります。

from ecell4.util import run_simulation
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=GillespieFactory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=ODEFactory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=SpatiocyteFactory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=BDFactory(bd_dt_factor=1))[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=MesoscopicFactory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', factory=EGFRDFactory())[-1])

[0.01, 0.0, 0.0, 60.0]
[0.01, 0.17972919304001073, 0.17972919304001067, 59.82027080696036]
[0.01, 0.0, 0.0, 60.0]
[0.01, 0.0, 0.0, 60.0]
[0.01, 0.0, 0.0, 60.0]
[0.01, 0.0, 0.0, 60.0]

!pip install ecell

Collecting ecell
  Downloading ecell-4.1.2-cp36-cp36m-manylinux1_x86_64.whl (45.5MB)
    100% |################################| 45.5MB 29kB/s  eta 0:00:01
Installing collected packages: ecell
Successfully installed ecell-4.1.2

5. シミュレーションを記録して視覚化する方法

ここでは、シミュレーション結果のログを取る方法と、それを視覚化する方法について説明します。

%matplotlib inline
import math
from ecell4 import *

5.1. Observer を用いたシミュレーションのログの取り方

E-Cell4は、 Observerというロギングのためのクラスを提供しています。 ObserverクラスはSimulatorの run関数を呼び出すときに与えられます。

def create_simulator(f=gillespie.GillespieFactory()):
    m = NetworkModel()
    A, B, C = Species('A', '0.005', '1'), Species('B', '0.005', '1'), Species('C', '0.005', '1')
    m.add_reaction_rule(create_binding_reaction_rule(A, B, C, 0.01))
    m.add_reaction_rule(create_unbinding_reaction_rule(C, A, B, 0.3))
    w = f.create_world()
    w.bind_to(m)
    w.add_molecules(C, 60)
    sim = f.create_simulator(w)
    sim.initialize()
    return sim

最も一般的な Observerの1つはFixedIntervalNumberObserverです。これは与えられた時間間隔で分子の数を記録します。 FixedIntervalNumberObserverは、時間間隔と、ログの対象とするのSpeciesのシリアルのリストを必要とします。

obs1 = FixedIntervalNumberObserver(0.1, ['A', 'B', 'C'])
sim = create_simulator()
sim.run(1.0, obs1)

FixedIntervalNumberObserverのdata関数はログに記録されたデータを返します。

print(obs1.data())

[[0.0, 0.0, 0.0, 60.0], [0.1, 1.0, 1.0, 59.0], [0.2, 3.0, 3.0, 57.0], [0.30000000000000004, 4.0, 4.0, 56.0], [0.4, 6.0, 6.0, 54.0], [0.5, 8.0, 8.0, 52.0], [0.6000000000000001, 9.0, 9.0, 51.0], [0.7000000000000001, 10.0, 10.0, 50.0], [0.8, 10.0, 10.0, 50.0], [0.9, 10.0, 10.0, 50.0], [1.0, 11.0, 11.0, 49.0]]

targets()はコンストラクタの引数として指定した Speciesのリストを返します。

print([sp.serial() for sp in obs1.targets()])

['A', 'B', 'C']

NumberObserverは、反応が起こったときのすべてのステップの後に分子の数を記録します。 このオブザーバーはすべての反応を記録するのに便利ですが、 odeでは利用できません。

obs1 = NumberObserver(['A', 'B', 'C'])
sim = create_simulator()
sim.run(1.0, obs1)
print(obs1.data())

[[0.0, 0.0, 0.0, 60.0], [0.043181315145277815, 1.0, 1.0, 59.0], [0.0999354928781079, 2.0, 2.0, 58.0], [0.10051642633328928, 3.0, 3.0, 57.0], [0.1576666300428589, 4.0, 4.0, 56.0], [0.22160562023979907, 5.0, 5.0, 55.0], [0.24838212056833223, 6.0, 6.0, 54.0], [0.38083606360319644, 7.0, 7.0, 53.0], [0.44115742167963257, 8.0, 8.0, 52.0], [0.478654279832333, 9.0, 9.0, 51.0], [0.524776672724245, 10.0, 10.0, 50.0], [0.5831861079171226, 11.0, 11.0, 49.0], [0.6016715021229188, 12.0, 12.0, 48.0], [0.6458184477531915, 13.0, 13.0, 47.0], [0.7164454853598512, 14.0, 14.0, 46.0], [0.8389493519306314, 15.0, 15.0, 45.0], [0.9250934945924608, 16.0, 16.0, 44.0], [0.9314408219193694, 17.0, 17.0, 43.0], [1.0, 17.0, 17.0, 43.0]]

TimingNumberObserverでは、ログの時刻をそのコンストラクタの引数として与えることができます。

obs1 = TimingNumberObserver([0.0, 0.1, 0.2, 0.5, 1.0], ['A', 'B', 'C'])
sim = create_simulator()
sim.run(1.0, obs1)
print(obs1.data())

[[0.0, 0.0, 0.0, 60.0], [0.1, 1.0, 1.0, 59.0], [0.2, 3.0, 3.0, 57.0], [0.5, 8.0, 8.0, 52.0], [1.0, 13.0, 13.0, 47.0]]

run関数は複数のオブザーバを一度に受け付けます。

obs1 = NumberObserver(['C'])
obs2 = FixedIntervalNumberObserver(0.1, ['A', 'B'])
sim = create_simulator()
sim.run(1.0, [obs1, obs2])
print(obs1.data())
print(obs2.data())

[[0.0, 60.0], [0.05163820372508536, 59.0], [0.10222920269451034, 58.0], [0.14249526069798346, 57.0], [0.1709401154327075, 56.0], [0.1786522132631355, 55.0], [0.32345579486892295, 54.0], [0.3320212006041937, 53.0], [0.37295796722048336, 52.0], [0.4297993097403142, 51.0], [0.442315663109651, 50.0], [0.456452361786781, 49.0], [0.759999628617218, 48.0], [0.7601731295722447, 47.0], [0.847168492783278, 46.0], [0.856013633150855, 47.0], [0.8614538661854448, 46.0], [0.870506834401381, 45.0], [0.978711403259659, 44.0], [1.0, 44.0]]
[[0.0, 0.0, 0.0], [0.1, 1.0, 1.0], [0.2, 5.0, 5.0], [0.30000000000000004, 5.0, 5.0], [0.4, 8.0, 8.0], [0.5, 11.0, 11.0], [0.6000000000000001, 11.0, 11.0], [0.7000000000000001, 11.0, 11.0], [0.8, 13.0, 13.0], [0.9, 15.0, 15.0], [1.0, 16.0, 16.0]]

FixedIntervalHDF5ObservedrはWorldのデータ全体を一定の間隔で出力ファイルに記録します。 2番目の引数は、出力ファイル名の接頭辞です。 filename()は、次に保存するようにスケジュールされたファイルの名前を返します。 ％02dのような多くのフォーマット文字列では、ファイル名にステップカウントを使用することができます。 書式文字列を使用しないと、最新のデータがファイルに上書きされます。

obs1 = FixedIntervalHDF5Observer(0.2, 'test%02d.h5')
print(obs1.filename())
sim = create_simulator()
sim.run(1.0, obs1)
print(obs1.filename())

test00.h5
test06.h5

w = load_world('test05.h5')
print(w.t(), w.num_molecules(Species('C')))

1.0 44

FixedIntervalCSVObserverの使用方法は、FixedIntervalHDF5Observerの使用方法とほぼ同じです。 FixedIntervalCSVObserverは半径(r)とSpeciesのシリアル番号(sid)を持つパーティクルの位置(x、y、z)をCSVファイルに保存します。

obs1 = FixedIntervalCSVObserver(0.2, "test%02d.csv")
print(obs1.filename())
sim = create_simulator()
sim.run(1.0, obs1)
print(obs1.filename())

test00.csv
test06.csv

下記は出力CSVファイルの最初の10行です。

print(''.join(open("test05.csv").readlines()[: 10]))

x,y,z,r,sid
0.17508278810419142,0.12197625380940735,0.53627523058094084,0,0
0.24078186159022152,0.62765785818919539,0.7254820850212127,0,0
0.40445487247779965,0.16720660566352308,0.92957371729426086,0,0
0.24083335651084781,0.59580284473486245,0.51183571317233145,0,0
0.90290508698672056,0.4137285016477108,0.39729581261053681,0,0
0.51044316007755697,0.26319809840060771,0.40705726365558803,0,0
0.84177179471589625,0.085374288959428668,0.64473599148914218,0,0
0.54177056765183806,0.098819603677839041,0.17240164172835648,0,0
0.28669063560664654,0.71829434810206294,0.91102162329480052,0,0

粒子シミュレーションのために、E-Cell4は、分子の軌道を追跡するためのオブザーバ（Observer）を提供し、これはFixedIntervalTrajectoryObserverと名付けられています。 ParticleIDが特に指定されていないときは、すべての軌跡を記録します。 シミュレーション中にパーティクルIDが反応のために失われると、FixedIntervalTrajectoryObserverはパーティクルをそれ以上トレースするのを単に止めます。

sim = create_simulator(spatiocyte.SpatiocyteFactory(0.005))
obs1 = FixedIntervalTrajectoryObserver(0.01)
sim.run(0.1, obs1)

print([tuple(pos) for pos in obs1.data()[0]])

[(0.46540305112880387, 0.739008344562721, 0.48), (0.2122891110412088, 0.9266471820493494, 0.665), (0.49806291436591293, 0.8689121551303868, 0.685), (0.7838367176906171, 0.8313843876330611, 0.64), (0.9961258287318259, 0.7967433714816836, 0.5), (0.7920016834998943, 0.9122134253196088, 0.63), (0.7348469228349536, 0.7534421012924616, 0.515), (0.9716309313039941, 0.851591647054698, 0.525), (1.1349302474895393, 0.9122134253196088, 0.63), (0.9961258287318259, 0.9439676901250381, 0.725), (0.718516991216399, 0.8833459118601275, 0.93)]

ほとんどの場合、Worldは各平面の周期的境界を想定しています。 境界条件に起因するエッジでのパーティクルの大きなジャンプを避けるために、 FixedIntervalTrajectoryObserverは位置のシフトを維持しようとします。 従って、Observerに格納されたパーティクルの位置は、Worldに対して与えられた立方体に必ずしも限定されません。 境界条件の拡散を正確に追跡するには、ロギングのステップ間隔を十分小さくする必要があります。 もちろん、このオプションは無効にすることができます。 help(FixedIntervalTrajectoryObserver)を参照してください。

5.2. ログデータの可視化

このセクションでは、Observerによって記録されたログデータの可視化ツールについて説明します。

まず、時系列データのためにviz.plot_number_observerはNumberObserverとFixedIntervalNumberObserverおよびTimingNumberObserverによって提供されるデータをプロットします。

viz.plot_number_observerの詳しい使用法については、help(viz.plot_number_observer)を参照してください。

obs1 = NumberObserver(['C'])
obs2 = FixedIntervalNumberObserver(0.1, ['A', 'B'])
sim = create_simulator()
sim.run(10.0, [obs1, obs2])

viz.plot_number_observer(obs1, obs2)

プロットのスタイルを設定したり、プロットに任意の関数を追加することもできます。

viz.plot_number_observer(obs1, '-', obs2, ':', lambda t: 60 * math.exp(-0.3 * t), '--')

位相平面でのプロットは、x軸とy軸を指定することによっても利用できます。

viz.plot_number_observer(obs2, 'o', x='A', y='B')

空間シミュレーションでは、Worldの状態を可視化するために、viz.plot_worldが利用可能です。 この関数は、インタラクティブな方法で3次元ボリューム内のパーティクルの点をプロットします。 描画領域の右ボタンをクリックすると、画像を保存できます。

sim = create_simulator(spatiocyte.SpatiocyteFactory(0.005))
# viz.plot_world(sim.world())
viz.plot_world(sim.world(), interactive=False)

FixedIntervalHDF5Observerとして与えられた一連のHDF5ファイルからムービーを作ることもできます。

注： viz.plot_movieは、オプションinteractive=Falseの時に追加ライブラリ ffmpegを必要とします。

sim = create_simulator(spatiocyte.SpatiocyteFactory(0.005))
obs1 = FixedIntervalHDF5Observer(0.02, 'test%02d.h5')
sim.run(1.0, obs1)
viz.plot_movie(obs1)

最後に、 FixedIntervalTrajectoryObserverに対応して、viz.plot_trajectoryはパーティクル軌道の可視化を提供します。

sim = create_simulator(spatiocyte.SpatiocyteFactory(0.005))
obs1 = FixedIntervalTrajectoryObserver(1e-3)
sim.run(1, obs1)
# viz.plot_trajectory(obs1)
viz.plot_trajectory(obs1, interactive=False)

6. 反応速度式の関数でODEを解く方法

odeソルバーはNetworkModelを受け入れますが、 ode.ODESimulatorはそのモデルクラスode.ODENetworkModelと ode.ODEReactionRuleを機能拡張のために所有しています。 これらのクラスのインタフェースは ModelクラスとReactionRuleとほとんど同じです。 ここでは、特に ode.ODERatelawに関するodeに特有の使用法について説明します。

%matplotlib inline
from ecell4 import *

しかし、ode における反応速度式のサポートはまだ開発中の状態です。 一部の関数は、将来廃止される可能性があります。 現在、反応速度式を有効にするには、オプション ecell4.util.decorator.ENABLE_RATELAWを次のように有効にする必要があります：

util.decorator.ENABLE_RATELAW = True

6.1. ode.ODEReactionRule

ode.ODEReactionRuleはReactionRuleとほぼ同じメンバーを持っています。

rr1 = ReactionRule()
rr1.add_reactant(Species("A"))
rr1.add_reactant(Species("B"))
rr1.add_product(Species("C"))
rr1.set_k(1.0)
print(len(rr1.reactants()))  # => 2
print(len(rr1.products()))  # => 1
print(rr1.k())  # => 1.0
print(rr1.as_string())  # => A+B>C|1

2
1
1.0
A+B>C|1

rr2 = ode.ODEReactionRule()
rr2.add_reactant(Species("A"))
rr2.add_reactant(Species("B"))
rr2.add_product(Species("C"))
rr2.set_k(1.0)
print(len(rr2.reactants()))  # => 2
print(len(rr2.products()))  # => 1
print(rr2.k())  # => 1.0
print(rr2.as_string())  # => A+B>C|1

2
1
1.0
A+B>C|1

共通のメンバーに加えて、 ode.ODEReactionRuleは各Speciesの化学量論的係数を格納することができます：

rr2 = ode.ODEReactionRule()
rr2.add_reactant(Species("A"), 1.0)
rr2.add_reactant(Species("B"), 2.0)
rr2.add_product(Species("C"), 2.5)
rr2.set_k(1.0)
print(rr2.as_string())

A+2*B>2.5*C|1

次のように係数にアクセスすることもできます。

print(rr2.reactants_coefficients())  # => [1.0, 2.0]
print(rr2.products_coefficients())  # => [2.5]

[1.0, 2.0]
[2.5]

6.2. ode.ODERatelaw

ode.ODEReactionRuleはode.ODERatelawにバインドされます。 ode.ODERatelawは与えられたSpeciesの値に基づいて導関数（フラックスまたは速度）を計算する関数を提供します。 ode.ODERatelawMassActionはode.ODEReactionRuleにバインドされたデフォルトのクラスです。

rr1 = ode.ODEReactionRule()
rr1.add_reactant(Species("A"))
rr1.add_reactant(Species("B"))
rr1.add_product(Species("C"))
rl1 = ode.ODERatelawMassAction(2.0)
rr1.set_ratelaw(rl1)  # equivalent to rr1.set_k(2.0)
print(rr1.as_string())

A+B>C|2

ode.ODERatelawCallbackは、フラックスを計算するためのユーザ定義関数を有効にします。

def mass_action(reactants, products, volume, t, rr):
    veloc = 2.0 * volume
    for value in reactants:
        veloc *= value / volume
    return veloc

rl2 = ode.ODERatelawCallback(mass_action)
rr1.set_ratelaw(rl2)
print(rr1.as_string())

A+B>C|mass_action

bound関数は5つの引数を受け入れ、浮動小数点を速度として返す必要があります。 第1および第2のリストは、それぞれ反応物および生成物のそれぞれの値を含みます。 化学量論係数にアクセスする必要があるときは、引数に rr（ode.ODEReactionRule）を使います。

ラムダ関数も利用可能です。

rl2 = ode.ODERatelawCallback(lambda r, p, v, t, rr: 2.0 * r[0] * r[1])
rr1.set_k(0)
rr1.set_ratelaw(rl2)
print(rr1.as_string())

A+B>C|<lambda>

6.3. ode.ODENetworkModel

ode.ODENetworkModelはReactionRuleと ode.ODEReactionRuleの両方を受け入れます。 ReactionRuleは暗黙的に変換され、ode.ODEReactionRuleとして保存されます。

m1 = ode.ODENetworkModel()
rr1 = create_unbinding_reaction_rule(Species("C"), Species("A"), Species("B"), 3.0)
m1.add_reaction_rule(rr1)
rr2 = ode.ODEReactionRule(create_binding_reaction_rule(Species("A"), Species("B"), Species("C"), 0.0))
rr2.set_ratelaw(ode.ODERatelawCallback(lambda r, p, v, t, rr: 0.1 * r[0] * r[1]))
m1.add_reaction_rule(rr2)

ode.ODENetworkModel中のode.ODEReactionRuleのリストには、モデルのメンバー reaction_rules（）を介してアクセスできます。

print([rr.as_string() for rr in m1.reaction_rules()])

['C>A+B|3', 'A+B>C|<lambda>']

最後に、次のように他のソルバーと同じ方法でシミュレーションを実行できます。

run_simulation(1.0, model=m1, y0={'A': 60, 'B': 60})

Python のデコレータを使用したモデリングは、レート（浮動小数）の代わりに関数を指定することによっても利用できます。 浮動小数点数が設定されている場合、それは質量作用反応の運動速度であると仮定されますが、一定速度ではありません。

with reaction_rules():
    A + B == C | (lambda r, *args: 0.1 * reduce(mul, r), 3.0)

m1 = get_model()

For the simplicity, you can directory defining the equation with Species names as follows:

with reaction_rules():
    A + B == C | (0.1 * A * B, 3.0)

m1 = get_model()

反応物や生成物としてリストアップされていないSpeciesを（速度則中で）呼び出すと、自動的に酵素としてリストに追加されます。

with reaction_rules():
    S > P | 1.0 * E * S / (30.0 + S)

m1 = get_model()
print(m1.reaction_rules()[0].as_string())

S+E>P+E|((1.0*E*S)/(30.0+S))

ここで式中のEは、反応物および生成物リストの両方に付加されます。

run_simulation(10.0, model=m1, y0={'S': 60, 'E': 30})

Speciesの名前の誤字に注意してください。 タイプミスをした場合、それは意図せず新しい酵素として認識されてしまいます：

with reaction_rules():
    A13P2G > A23P2G | 1500 * A13B2G  # typo: A13P2G -> A13B2G

m1 = get_model()
print(m1.reaction_rules()[0].as_string())

A13P2G+A13B2G>A23P2G+A13B2G|(1500*A13B2G)

酵素の自動宣言を無効にしたい場合は、 util.decorator.ENABLE_IMPLICIT_DECLARATIONを無効にしてください。 その値が Falseの場合、上記の場合はエラーが発生します：

util.decorator.ENABLE_IMPLICIT_DECLARATION = False

try:
    with reaction_rules():
        A13P2G > A23P2G | 1500 * A13B2G
except RuntimeError as e:
    print(repr(e))

util.decorator.ENABLE_IMPLICIT_DECLARATION = True

RuntimeError('unknown variable [A13B2G] was used.',)

E-Cell4は、生化学反応ネットワークのシミュレーションに特化していますが、合成反応ルールを使用することで、常微分方程式を直感的に変換することができます。 たとえば、Lotka-Volterraの式は次のようになります。

ここで は以下のように解かれます。

with reaction_rules():
    A, B, C, D = 1.5, 1, 3, 1

    ~x > x | A * x - B * x * y
    ~y > y | -C * y + D * x * y

run_simulation(10, model=get_model(), y0={'x': 10, 'y': 5})

6.4. 速度則における参照について

ここで、速度則の定義の詳細を説明します。

第1に、Species で速度則のより単純な定義を使用すると、限られた数の数学関数（例えば exp,log,sin,cos,tan,asin,acos,atan,やpiなど） はブロック外で関数を宣言してもそこで利用できます。

try:
    from math import erf

    with reaction_rules():
        S > P | erf(S / 30.0)
except TypeError as e:
    print(repr(e))

TypeError('a float is required',)

This is because erf is tried to be evaluated agaist S / 30.0 first, but it is not a floating number. In contrast, the following case is acceptable:

from math import erf

with reaction_rules():
    S > P | erf(2.0) * S

m1 = get_model()
print(m1.reaction_rules()[0].as_string())

S>P|(0.995322265019*S)

erf（2.0）、 0.995322265019の結果だけがレート法に渡されます。 したがって、上記のレート法は erf関数への参照を持ちません。 同様に、外部に宣言された変数の値は受け入れ可能ですが、参照としては使用できません。

kcat, Km = 1.0, 30.0

with reaction_rules():
    S > P | kcat * E * S / (Km + S)

m1 = get_model()
print(m1.reaction_rules()[0].as_string())
kcat = 2.0
print(m1.reaction_rules()[0].as_string())

S+E>P+E|((1.0*E*S)/(30.0+S))
S+E>P+E|((1.0*E*S)/(30.0+S))

変数の値を変更しても、速度則には影響しません。 一方、独自の関数を使用して速度則を定義すると、外部の変数への参照を保持することができます。

k1 = 1.0

with reaction_rules():
    S > P | (lambda r, *args: k1 * r[0])  # referring k1

m1 = get_model()

obs1 = run_simulation(2, model=m1, y0={"S": 60}, return_type='observer')
k1 = 2.0
obs2 = run_simulation(2, model=m1, y0={"S": 60}, return_type='observer')

viz.plot_number_observer(obs1, '-', obs2, '--')

ただし、この場合、パラメータのセットごとに新しいモデルを作成する方が良いです。

def create_model(k):
    with reaction_rules():
        S > P | k

    return get_model()

obs1 = run_simulation(2, model=create_model(k=1.0), y0={"S": 60}, return_type='observer')
obs2 = run_simulation(2, model=create_model(k=2.0), y0={"S": 60}, return_type='observer')
# viz.plot_number_observer(obs1, '-', obs2, '--')

6.5. odeについてより深く

ode.ODEWorldでは、各Speciesの値は浮動小数です。 しかし、互換性のために、共通メンバ num_moleculesとadd_moleculesはその値を整数とみなします。

w = ode.ODEWorld()
w.add_molecules(Species("A"), 2.5)
print(w.num_molecules(Species("A")))

2

実数を設定/取得するには、 set_valueとget_valueを使います：

w.set_value(Species("B"), 2.5)
print(w.get_value(Species("A")))
print(w.get_value(Species("B")))

2.0
2.5

デフォルトでは ode.ODESimulatorはROSENBROCK4_CONTROLLERというRosenblock法を使用してODEを解きます。 それに加えて、2つのソルバー、 EULERとRUNGE_KUTTA_CASH_KARP54が利用可能です。 ROSENBROCK4_CONTROLLERとRUNGE_KUTTA_CASH_KARP54はエラー制御により時間発展中のステップサイズを適応的に変更しますが、 EULERはこれを行いません。

with reaction_rules():
    A > ~A | 1.0

m1 = get_model()

w1 = ode.ODEWorld()
w1.set_value(Species("A"), 1.0)
sim1 = ode.ODESimulator(m1, w1, ode.EULER)
sim1.set_dt(0.01) # This is only effective for EULER
sim1.run(3.0, obs1)

ode.ODEFactoryは、ソルバ型とデフォルトのステップ間隔も受け入れます。

run_simulation(3.0, model=m1, y0={"A": 1.0}, factory=ode.ODEFactory(ode.EULER, 0.01))

下記の例も参照してください。

• Glycolysis of Human Erythrocytes
• Drosophila Circadian Clock
• Attractors

7. ルールベースモデリング入門

E-Cell4は、ルールベースモデリング環境を提供します。

%matplotlib inline
from ecell4 import *

7.1. Species.count

まず、Speciesはcount関数を提供します。

Species.count は Species の間でマッチした数を数えます。

sp1 = Species("A(b^1).B(b^1)")
sp2 = Species("A(b^1).A(b^1)")
pttrn1 = Species("A")
print(pttrn1.count(sp1))  # => 1
print(pttrn1.count(sp2))  # => 2

1
2

上記の場合、Species.countはサイトには関係なく Species の A という名前の UnitSpecies の数を返します。 結合の有無の状態を明確にする場合は下記のようになります。

pttrn1 = Species("A(b)")
pttrn2 = Species("A(b^_)")
print(pttrn1.count(sp1))  # => 0
print(pttrn2.count(sp1))  # => 1

0
1

ここでA(b)は結合bが空いていることを示唆し、A(b^_)は結合bが有ることを意味します。 アンダースコア _はワイルドカードを意味します。 同様に、サイトの状態を指定することもできます。

sp1 = Species("A(b=u)")
pttrn1 = Species("A(b)")
pttrn2 = Species("A(b=u)")
print(pttrn1.count(sp1))  # => 1
print(pttrn2.count(sp1))  # => 1

1
1

A(b)は状態については何も言いませんが、 A(b=u)は状態と結合の両方をその分子の条件として明記します。 A(b=u)は、Aという名前の UnitSpeciesが状態がuで結合が空である bというサイトを持つことを意味します。 ワイルドカード _は状態と名前両方で用いることが可能です。

sp1 = Species("A(b=u^1).B(b=p^1)")
pttrn1 = Species("A(b=_^_)")  # This is equivalent to `A(b^_)` here
pttrn2 = Species("_(b^_)")
print(pttrn1.count(sp1))  # => 1
print(pttrn2.count(sp1))  # => 2

1
2

ワイルドカード _ は何にでもマッチし、マッチしたパターンは Species でもワイルドカード間で一貫しません。

一方、名前付きのワイルドカード _1 、_2 などは、マッチ内の一貫性を与えます。

sp1 = Species("A(b^1).B(b^1)")
pttrn1 = Species("_._")
pttrn2 = Species("_1._1")
print(pttrn1.count(sp1))  # => 2
print(pttrn2.count(sp1))  # => 0

2
0

最初のパターンは2つの方法（A.B と B.A）でマッチしますが、2番目のパターンは何もマッチしません。

Species.countは対称の場合でも常に UnitSpecies の順序を区別します。 したがって、 _1._1 は二量体の数を意味しません。

sp1 = Species("A(b^1).A(b^1)")
pttrn1 = Species("_1._1")
print(pttrn1.count(sp1))  # => 2

2

7.2. ReactionRule.count と ReactionRule.generate

ReactionRule はまた、与えられたreactantsのリストとのマッチ数をカウントする機能を持っています。

rr1 = create_unimolecular_reaction_rule(Species("A(p=u)"), Species("A(p=p)"), 1.0)
sp1 = Species("A(b^1,p=u).B(b^1)")
print(rr1.count([sp1]))  # => 1

1

ReactionRule.generate は、マッチ情報に基づいて生成された ReactionRule のリストを返します。

print([rr.as_string() for rr in rr1.generate([sp1])])

[u'A(b^1,p=u).B(b^1)>A(b^1,p=p).B(b^1)|1']

ReactionRule.generate は与えられたリストの中のSpecies の順番で結果が違ってきます：

rr1 = create_binding_reaction_rule(Species("A(b)"), Species("B(b)"), Species("A(b^1).B(b^1)"), 1.0)
sp1 = Species("A(b)")
sp2 = Species("B(b)")
print([rr.as_string() for rr in rr1.generate([sp1, sp2])])
print([rr.as_string() for rr in rr1.generate([sp2, sp1])])

[u'A(b)+B(b)>A(b^1).B(b^1)|1']
[]

現在の実装では、 ReactionRule.generate は常に一意の ReactionRule のリストを返すわけではありません。

sp1 = Species("A(b,c^1).A(b,c^1)")
sp2 = Species("B(b,c^1).B(b,c^1)")
print(rr1.count([sp1, sp2]))  # => 4
print([rr.as_string() for rr in rr1.generate([sp1, sp2])])

4
[u'A(b,c^1).A(b,c^1)+B(b,c^1).B(b,c^1)>A(b^1,c^2).A(b,c^2).B(b^1,c^3).B(b,c^3)|1', u'A(b,c^1).A(b,c^1)+B(b,c^1).B(b,c^1)>A(b^1,c^2).A(b,c^2).B(b,c^3).B(b^1,c^3)|1', u'A(b,c^1).A(b,c^1)+B(b,c^1).B(b,c^1)>A(b,c^1).A(b^2,c^1).B(b^2,c^3).B(b,c^3)|1', u'A(b,c^1).A(b,c^1)+B(b,c^1).B(b,c^1)>A(b,c^1).A(b^2,c^1).B(b,c^3).B(b^2,c^3)|1']

上に挙げた ReactionRules は異なって見えますが、すべての結果は同じことを示しています。

print(set([unique_serial(rr.products()[0]) for rr in rr1.generate([sp1, sp2])]))

set([u'A(b,c^1).A(b^2,c^1).B(b^2,c^3).B(b,c^3)'])

これは、これらの ReactionRule は、同じ Species を生成するとはいえ、異なるマッチ情報に基づいて生成されるために起こります。

詳細については、以下のセクションを参照してください。

ワイルドカードは ReactionRule でも利用できます。

rr1 = create_unimolecular_reaction_rule(Species("A(p=u^_)"), Species("A(p=p^_)"), 1.0)
print([rr.as_string() for rr in rr1.generate([Species("A(p=u^1).B(p^1)")])])

[u'A(p=u^1).B(p^1)>A(p=p^1).B(p^1)|1']

もちろん、ワイルドカードは UnitSpeciesの名前として受け入れられます。

rr1 = create_unimolecular_reaction_rule(Species("_(p=u)"), Species("_(p=p)"), 1.0)
print([rr.as_string() for rr in rr1.generate([Species("A(p=u)")])])
print([rr.as_string() for rr in rr1.generate([Species("B(b^1,p=u).C(b^1,p=u)")])])

[u'A(p=u)>A(p=p)|1']
[u'B(b^1,p=u).C(b^1,p=u)>B(b^1,p=p).C(b^1,p=u)|1', u'B(b^1,p=u).C(b^1,p=u)>B(b^1,p=u).C(b^1,p=p)|1']

State中での名前付きワイルドカードは、サイト間の対応を指定するのに便利です。

rr1 = create_unbinding_reaction_rule(Species("AB(a=_1, b=_2)"), Species("B(b=_2)"), Species("A(a=_1)"), 1.0)
print([rr.as_string() for rr in rr1.generate([Species("AB(a=x, b=y)")])])
print([rr.as_string() for rr in rr1.generate([Species("AB(a=y, b=x)")])])

[u'AB(a=x,b=y)>B(b=y)+A(a=x)|1']
[u'AB(a=y,b=x)>B(b=x)+A(a=y)|1']

無名ワイルドカード _ は、マッチの間の等価性を気にしません。

Products は順番に生成されます。

rr1 = create_binding_reaction_rule(Species("_(b)"), Species("_(b)"), Species("_(b^1)._(b^1)"), 1.0)
print(rr1.as_string())
print([rr.as_string() for rr in rr1.generate([Species("A(b)"), Species("A(b)")])])
print([rr.as_string() for rr in rr1.generate([Species("A(b)"), Species("B(b)")])])

_(b)+_(b)>_(b^1)._(b^1)|1
[u'A(b)+A(b)>A(b^1).A(b^1)|1']
[u'A(b)+B(b)>A(b^1).B(b^1)|1']

その対称性のために、上記のセルの前者の場合は元の反応速度の半分となるように設定されることが望ましい場合があります。

これは 後者の分子の組み合わせの数は ではありますが，前者の分子の組み合わせの数が であるためです。 (ここでAとBの両方の分子の数は です。) またこれは gillespie と ode に当てはまります。 しかし、egfrd と spatiocyte では、反応速度は固有のものであり、組み合わせの数には影響されません (したがって、E-Cell4では、その場合でも速度の変更は行われません。)。

アルゴリズムの違いについては、 Homodimerization and Annihilation を参照してください。

名前のないワイルドカードとは対照的に、名前付けされたワイルドカードは名前付けされたワイルドカード同士での整合性を保持し、 ReactionRule で常に同じ値を示します。

rr1 = create_binding_reaction_rule(Species("_1(b)"), Species("_1(b)"), Species("_1(b^1)._1(b^1)"), 1.0)
print(rr1.as_string())
print([rr.as_string() for rr in rr1.generate([Species("A(b)"), Species("A(b)")])])
print([rr.as_string() for rr in rr1.generate([Species("A(b)"), Species("B(b)")])])  # => []

_1(b)+_1(b)>_1(b^1)._1(b^1)|1
[u'A(b)+A(b)>A(b^1).A(b^1)|1']
[]

名前付きのワイルドカードは、 UnitSpecies と site の名前の間でも厳密に整合性を保持しています。

rr1 = create_binding_reaction_rule(Species("A(b=_1)"), Species("_1(b)"), Species("A(b=_1^1)._1(b^1)"), 1.0)
print(rr1.as_string())
print([rr.as_string() for rr in rr1.generate([Species("A(b=B)"), Species("A(b)")])])  # => []
print([rr.as_string() for rr in rr1.generate([Species("A(b=B)"), Species("B(b)")])])

A(b=_1)+_1(b)>A(b=_1^1)._1(b^1)|1
[]
[u'A(b=B)+B(b)>A(b=B^1).B(b^1)|1']

7.3. NetfreeModel

NetfreeModelはルールベースモデリングのためのModelクラスです。

NetfreeModel のインターフェースは NetworkModel とほぼ同じですが、ルールとマッチを考慮しています。

rr1 = create_binding_reaction_rule(Species("A(r)"), Species("A(l)"), Species("A(r^1).A(l^1)"), 1.0)

m1 = NetfreeModel()
m1.add_reaction_rule(rr1)
print(m1.num_reaction_rules())

m2 = NetworkModel()
m2.add_reaction_rule(rr1)
print(m2.num_reaction_rules())

1
1

2. モデルを構築する方法 で説明したPythonの表記法はここでも利用可能です。

モデルを NetfreeModel として取得するには、get_model 関数で is_netfree オプションを True に設定します。

with reaction_rules():
    A(r) + A(l) > A(r^1).A(l^1) | 1.0

m1 = get_model(is_netfree=True)
print(repr(m1))

<ecell4.core.NetfreeModel object at 0x7fbb44dc7708>

Model.query_reaction_rules は、与えられた反応物に対する ReactionRule のリストを返します。

NetworkModel は Species の等価性に基づいて ReactionRule を返します。

print(len(m2.query_reaction_rules(Species("A(r)"), Species("A(l)"))))  # => 1
print(len(m2.query_reaction_rules(Species("A(l,r)"), Species("A(l,r)"))))  # => 0

1
0

一方、 NetfreeModel は、格納された ReactionRule をルールベースの方法に適用することによってリストを生成します。

print(len(m1.query_reaction_rules(Species("A(r)"), Species("A(l)"))))  # => 1
print(len(m1.query_reaction_rules(Species("A(l,r)"), Species("A(l,r)"))))  # => 1

1
1

NetfreeModel は生成されたオブジェクトをキャッシュしません。 したがって、 NetfreeModel.query_reaction_rules は遅いですが、より少ないメモリしか必要としません。

print(m1.query_reaction_rules(Species("A(l,r)"), Species("A(l,r)"))[0].as_string())
print(m1.query_reaction_rules(Species("A(l,r^1).A(l^1,r)"), Species("A(l,r)"))[0].as_string())
print(m1.query_reaction_rules(Species("A(l,r^1).A(l^1,r)"), Species("A(l,r^1).A(l^1,r)"))[0].as_string())

A(l,r)+A(l,r)>A(l,r^1).A(l^1,r)|2
A(l,r^1).A(l^1,r)+A(l,r)>A(l,r^1).A(l^1,r^2).A(l^2,r)|2
A(l,r^1).A(l^1,r)+A(l,r^1).A(l^1,r)>A(l,r^1).A(l^1,r^2).A(l^2,r^3).A(l^3,r)|2

NetfreeModel.expand は与えられたシードに対して繰り返し Reaction Rules を適用することによって NetfreeModelをNetworkModelに展開します。

with reaction_rules():
    _(b) + _(b) == _(b^1)._(b^1) | (1.0, 1.0)

m3 = get_model(True)
print(m3.num_reaction_rules())

m4 = m3.expand([Species("A(b)"), Species("B(b)")])
print(m4.num_reaction_rules())
print([rr.as_string() for rr in m4.reaction_rules()])

2
6
[u'A(b)+A(b)>A(b^1).A(b^1)|1', u'A(b)+B(b)>A(b^1).B(b^1)|1', u'B(b)+B(b)>B(b^1).B(b^1)|1', u'A(b^1).A(b^1)>A(b)+A(b)|1', u'A(b^1).B(b^1)>A(b)+B(b)|1', u'B(b^1).B(b^1)>B(b)+B(b)|1']

展開の無限ループを避けるために、繰り返しの上限と Species における UnitSpecies の数の上限を制限することができます。

m2 = m1.expand([Species("A(l, r)")], 100, {Species("A"): 4})
print(m2.num_reaction_rules())  # => 4
print([rr.as_string() for rr in m2.reaction_rules()])

4
[u'A(l,r)+A(l,r)>A(l,r^1).A(l^1,r)|2', u'A(l,r^1).A(l^1,r)+A(l,r^1).A(l^1,r)>A(l,r^1).A(l^1,r^2).A(l^2,r^3).A(l^3,r)|2', u'A(l,r)+A(l,r^1).A(l^1,r)>A(l,r^1).A(l^1,r^2).A(l^2,r)|2', u'A(l,r)+A(l,r^1).A(l^1,r^2).A(l^2,r)>A(l,r^1).A(l^1,r^2).A(l^2,r^3).A(l^3,r)|2']

7.4. Species, ReactionRule そして NetfreeModel の違い

上記で説明した機能は似ていますが、基準にはいくつかの違いがあります。

sp1 = Species("A(b^1).A(b^1)")
sp2 = Species("A(b)")
rr1 = create_unbinding_reaction_rule(sp1, sp2, sp2, 1.0)
print(sp1.count(sp1))
print([rr.as_string() for rr in rr1.generate([sp1])])

2
[u'A(b^1).A(b^1)>A(b)+A(b)|1']

Species.count はマッチングのための2つの異なる方向（左から右 と 右から左）を示唆していますが、 ReactionRule.generate はたった1つの ReactionRuleしか返しません。なぜならその順は生成物に影響しないからです。

sp1 = Species("A(b^1).B(b^1)")
rr1 = create_unbinding_reaction_rule(
    sp1, Species("A(b)"), Species("B(b)"), 1.0)
sp2 = Species("A(b^1,c^2).A(b^3,c^2).B(b^1).B(b^3)")
print(sp1.count(sp2))
print([rr.as_string() for rr in rr1.generate([sp2])])

2
[u'A(b^1,c^2).A(b^3,c^2).B(b^1).B(b^3)>A(b,c^1).A(b^2,c^1).B(b^2)+B(b)|1', u'A(b^1,c^2).A(b^3,c^2).B(b^1).B(b^3)>A(b^1,c^2).A(b,c^2).B(b^1)+B(b)|1']

この場合、 ReactionRule.generate は Species.count と同様に動作します。

しかし、 Netfree.query_reaction_rules はより速い反応速度を持つたった一つの ReationRule を返します：

m1 = NetfreeModel()
m1.add_reaction_rule(rr1)
print([rr.as_string() for rr in m1.query_reaction_rules(sp2)])

[u'A(b^1,c^2).A(b^3,c^2).B(b^1).B(b^3)>A(b,c^1).A(b^2,c^1).B(b^2)+B(b)|2']

NetfreeModel.query_reaction_rules は、生成された各 ReactionRule が他と同じかどうかをチェックし、可能ならばそれを集約します。

上で説明したように、 ReactionRule.generate は Species の順序次第で結果が異なりますが、 Netfree.query_reaction_rulesはそうではありません。

sp1 = Species("A(b)")
sp2 = Species("B(b)")
rr1 = create_binding_reaction_rule(sp1, sp2, Species("A(b^1).B(b^1)"), 1.0)
m1 = NetfreeModel()
m1.add_reaction_rule(rr1)
print([rr.as_string() for rr in rr1.generate([sp1, sp2])])
print([rr.as_string() for rr in m1.query_reaction_rules(sp1, sp2)])
print([rr.as_string() for rr in rr1.generate([sp2, sp1])])  # => []
print([rr.as_string() for rr in m1.query_reaction_rules(sp2, sp1)])

[u'A(b)+B(b)>A(b^1).B(b^1)|1']
[u'A(b)+B(b)>A(b^1).B(b^1)|1']
[]
[u'B(b)+A(b)>A(b^1).B(b^1)|1']

名前のないワイルドカードは、一致していなければなりませんが、名前のないワイルドカードは必ずしも一致している必要はありません。

sp1 = Species("_(b)")
sp2 = Species("_1(b)")
sp3 = Species("A(b)")
sp4 = Species("B(b)")
rr1 = create_binding_reaction_rule(sp1, sp1, Species("_(b^1)._(b^1)"), 1)
rr2 = create_binding_reaction_rule(sp2, sp2, Species("_1(b^1)._1(b^1)"), 1)
print(sp1.count(sp2))  # => 1
print(sp1.count(sp3))  # => 1
print(sp2.count(sp2))  # => 1
print(sp2.count(sp3))  # => 1
print([rr.as_string() for rr in rr1.generate([sp3, sp3])])
print([rr.as_string() for rr in rr1.generate([sp3, sp4])])
print([rr.as_string() for rr in rr2.generate([sp3, sp3])])
print([rr.as_string() for rr in rr2.generate([sp3, sp4])])  # => []

1
1
1
1
[u'A(b)+A(b)>A(b^1).A(b^1)|1']
[u'A(b)+B(b)>A(b^1).B(b^1)|1']
[u'A(b)+A(b)>A(b^1).A(b^1)|1']
[]

8. 「1. E-Cell4を用いたシミュレーションの概要」についてより詳しく

一度1. E-Cell4を用いたシミュレーションの概要を読んだ後は、World と Simulator を使うのは難しくありません。 volume と {『C』 : 60} は World に相当し、ソルバーは以下の Simulator に相当します。

%matplotlib inline
from ecell4 import *

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

y = run_simulation(10.0, {'C': 60}, volume=1.0)

ここでは run_simulation がその内部で行っていることを分解して示します。 run_simulation はデフォルトでODEシミュレータを使用するので、ODEWorld を段階的に作成してみましょう。

8.1. ODEWorld の作成

World はこのようにして作ることができます。

w = ode.ODEWorld(Real3(1, 1, 1))

Real3 は座標ベクトルです。

この例では、 ODEWorld コンストラクタの最初の引数はキューブです。

run_simulation 引数のように、ode.ODEWorld 引数にはvolumeを使用できないことに注意してください。

今度はシミュレーション用のキューブボックスを作成した後、分子をキューブに投げ込みましょう。

w = ode.ODEWorld(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
print(w.t(), w.num_molecules(Species('C')))  # must return (0.0, 60)

(0.0, 60)

分子を追加するには add_molecules、分子を削除するにはremove_molecules、分子の数を知るには num_molecules を使います。

各メソッドの第一引数は、あなたが知りたい Species です。

t メソッドで現在時刻を取得できます。

しかし、ODEソルバーの分子数は実数であり、これらの _molecules 関数では整数に対してのみ動作します。

ODEで実数を扱うときは、 set_value と get_value を使います。

8.2. Real3 の使い方

Simulator の詳細の前に、Real3 について詳しく説明します。

pos = Real3(1, 2, 3)
print(pos)  # must print like <ecell4.core.Real3 object at 0x7f44e118b9c0>
print(tuple(pos))  # must print (1.0, 2.0, 3.0)

<ecell4.core.Real3 object at 0x7f05b0e6f630>
(1.0, 2.0, 3.0)

Real3 オブジェクトの内容を直接 print することはできません。

Real3 をPythonタプルまたはリストに一度変換する必要があります。

pos1 = Real3(1, 1, 1)
x, y, z = pos[0], pos[1], pos[2]
pos2 = pos1 + pos1
pos3 = pos1 * 3
pos4 = pos1 / 5
print(length(pos1))  # must print 1.73205080757
print(dot_product(pos1, pos3))  # must print 9.0

1.73205080757
9.0

dot_product のような基本的な関数も使えます。

もちろん、Real3 をnumpyの配列に変換することもできます。

import numpy
a = numpy.asarray(tuple(Real3(1, 2, 3)))
print(a)  # must print [ 1.  2.  3.]

[ 1.  2.  3.]

Integer3 は整数のトリプレットになります。

g = Integer3(1, 2, 3)
print(tuple(g))

(1, 2, 3)

もちろん、単純な算術演算を Integer3 に適用することもできます。

print(tuple(Integer3(1, 2, 3) + Integer3(4, 5, 6)))  # => (5, 7, 9)
print(tuple(Integer3(4, 5, 6) - Integer3(1, 2, 3)))  # => (3, 3, 3)
print(tuple(Integer3(1, 2, 3) * 2))  # => (2, 4, 6)
print(dot_product(Integer3(1, 2, 3), Integer3(4, 5, 6)))  # => 32
print(length(Integer3(1, 2, 3)))  # => 3.74165738677

(5, 7, 9)
(3, 3, 3)
(2, 4, 6)
32
3.74165738677

8.3. ODESimulator の作成と実行

下記のように Model と World を使って Simulatorを作成することができます。

with reaction_rules():
    A + B > C | 0.01  # equivalent to create_binding_reaction_rule
    C > A + B | 0.3   # equivalent to create_unbinding_reaction_rule

m = get_model()

sim = ode.ODESimulator(m, w)
sim.run(10.0)

run メソッドを呼び出すと、シミュレーションが実行されます。

この例では、シミュレーションは10秒間実行されます。

また World の状態は下記のようにしてチェックできます。

print(w.t(), w.num_molecules(Species('C')))  # must return (10.0, 30)

(10.0, 30)

Species C の数が60から30に減少していることがわかります。

World はある時点での状態を表しているため、World でシミュレーションの遷移を見ることはできません。

時系列の結果を得るには、 Observer を使います。

w = ode.ODEWorld(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
sim = ode.ODESimulator(m, w)

obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)
print(obs.data())  # must return [[0.0, 0.0, 60.0], ..., [10.0, 29.994446899691276, 30.005553100308752]]

[[0.0, 0.0, 60.0], [0.1, 1.7722206098711988, 58.227779390128795], [0.2, 3.4860124889661757, 56.51398751103382], [0.30000000000000004, 5.1376332715496495, 54.862366728450354], [0.4, 6.724090809612153, 53.27590919038786], [0.5, 8.243129756755453, 51.75687024324456], [0.6000000000000001, 9.69320376680592, 50.3067962331941], [0.7000000000000001, 11.073435590968808, 48.92656440903121], [0.8, 12.383567691608423, 47.616432308391595], [0.9, 13.62390591657045, 46.37609408342957], [1.0, 14.795258681171735, 45.20474131882829], [1.1, 15.898873899780316, 44.10112610021971], [1.2000000000000002, 16.936375633755194, 43.06362436624483], [1.3, 17.90970211303103, 42.090297886969], [1.4000000000000001, 18.821046466966358, 41.17895353303366], [1.5, 19.67280118133671, 40.32719881866331], [1.6, 20.46750699887015, 39.53249300112988], [1.7000000000000002, 21.207806714875233, 38.792193285124796], [1.8, 21.896404094818397, 38.10359590518163], [1.9000000000000001, 22.536027939969333, 37.463972060030684], [2.0, 23.129401186565122, 36.8705988134349], [2.1, 23.679214801318086, 36.320785198681946], [2.2, 24.188106157761755, 35.81189384223828], [2.3000000000000003, 24.658641522840725, 35.34135847715931], [2.4000000000000004, 25.093302252968826, 34.9066977470312], [2.5, 25.49447428958285, 34.50552571041718], [2.6, 25.86444054777012, 34.13555945222991], [2.7, 26.205375806607968, 33.79462419339207], [2.8000000000000003, 26.51934373464207, 33.48065626535796], [2.9000000000000004, 26.80829570797557, 33.19170429202446], [3.0, 27.07407111744664, 32.92592888255339], [3.1, 27.318398885233666, 32.68160111476636], [3.2, 27.542899949227504, 32.457100050772524], [3.3000000000000003, 27.749090501388586, 32.25090949861144], [3.4000000000000004, 27.938385796523626, 32.0616142034764], [3.5, 28.112104375214862, 31.88789562478517], [3.6, 28.271472567998387, 31.728527432001645], [3.7, 28.417629170003387, 31.582370829996645], [3.8000000000000003, 28.551630196373424, 31.44836980362661], [3.9000000000000004, 28.674453642480778, 31.32554635751925], [4.0, 28.787004191633095, 31.212995808366934], [4.1000000000000005, 28.890117821767866, 31.109882178232162], [4.2, 28.984566277309003, 31.015433722691025], [4.3, 29.071061377146354, 30.928938622853675], [4.4, 29.150259141926327, 30.8497408580737], [4.5, 29.22276372627192, 30.777236273728107], [4.6000000000000005, 29.289131148912805, 30.710868851087223], [4.7, 29.349872817459094, 30.650127182540935], [4.800000000000001, 29.405458847152726, 30.594541152847302], [4.9, 29.456321176884504, 30.543678823115524], [5.0, 29.50285648638057, 30.49714351361946], [5.1000000000000005, 29.545428921493947, 30.45457107850608], [5.2, 29.58437263404594, 30.415627365954087], [5.300000000000001, 29.61999414522251, 30.38000585477752], [5.4, 29.652574540327944, 30.347425459672085], [5.5, 29.682371504475235, 30.317628495524794], [5.6000000000000005, 29.70962120751236, 30.29037879248767], [5.7, 29.734540047958276, 30.265459952041752], [5.800000000000001, 29.757326264061593, 30.242673735938435], [5.9, 29.778161421010864, 30.221838578989164], [6.0, 29.797211782455758, 30.20278821754427], [6.1000000000000005, 29.814629574218237, 30.18537042578179], [6.2, 29.8305541480645, 30.169445851935528], [6.300000000000001, 29.84511305246321, 30.15488694753682], [6.4, 29.858423017246, 30.141576982754028], [6.5, 29.87059085870597, 30.129409141294058], [6.6000000000000005, 29.881714310971407, 30.11828568902862], [6.7, 29.89188278945169, 30.10811721054834], [6.800000000000001, 29.90117809152673, 30.098821908473298], [6.9, 29.90967503947212, 30.09032496052791], [7.0, 29.917442070090726, 30.082557929909303], [7.1000000000000005, 29.92454177537786, 30.075458224622167], [7.2, 29.931031398107503, 30.068968601892525], [7.300000000000001, 29.936963286002193, 30.063036713997835], [7.4, 29.942385307810326, 30.057614692189702], [7.5, 29.947341234455514, 30.052658765544514], [7.6000000000000005, 29.951871088073982, 30.048128911926046], [7.7, 29.956011461640173, 30.043988538359855], [7.800000000000001, 29.959795811548137, 30.04020418845189], [7.9, 29.963254725453346, 30.036745274546682], [8.0, 29.966416167388665, 30.033583832611363], [8.1, 29.969305702114955, 30.030694297885073], [8.200000000000001, 29.971946700366555, 30.028053299633473], [8.3, 29.97436052665158, 30.02563947334845], [8.4, 29.9765667110567, 30.023433288943327], [8.5, 29.978583106421084, 30.021416893578944], [8.6, 29.98042603207785, 30.019573967922177], [8.700000000000001, 29.982110405343864, 30.017889594656165], [8.8, 29.983649861761148, 30.01635013823888], [8.9, 29.985056865065513, 30.014943134934516], [9.0, 29.98634280775227, 30.01365719224776], [9.1, 29.9875181030258, 30.01248189697423], [9.200000000000001, 29.988592268884133, 30.011407731115895], [9.3, 29.98957400499266, 30.010425995007367], [9.4, 29.990471262976794, 30.009528737023235], [9.5, 29.99129131068186, 30.00870868931817], [9.600000000000001, 29.992040790927835, 30.007959209072194], [9.700000000000001, 29.99272577521901, 30.007274224781018], [9.8, 29.993351812847337, 30.00664818715269], [9.9, 29.993923975779538, 30.00607602422049], [10.0, 29.994446899691276, 30.005553100308752]]

E-Cell4にはいくつかのタイプの Observer があります。

FixedIntervalNumberObserver は、時系列の結果を得るための最もシンプルな Observer です。

その名前が示唆するように、この Observer は各時間ステップの分子の数を記録します。

第1引数は時間ステップ、第2引数は分子種です。

その結果は data メソッドで確認できますが、これにはショートカットがあります。

viz.plot_number_observer(obs)

上記は時系列の結果を簡単にプロットします。

ここでは run_simulation 関数の内部について説明しました。

Simulator を作成した後に World を変更した場合は Simulator にそれを示す必要があります。

またそのあとで sim.initialize() を呼び出すことを忘れないでください。

8.4. ソルバーの差し替えについて

run_simulation の内部を上記で示したので、ソルバーを確率論的手法に替えることは難しくないでしょう。

from ecell4 import *

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

m = get_model()

# ode.ODEWorld -> gillespie.GillespieWorld
w = gillespie.GillespieWorld(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)

# ode.ODESimulator -> gillespie.GillespieSimulator
sim = gillespie.GillespieSimulator(m, w)
obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)

viz.plot_number_observer(obs)

World と Simulator は Model 自体を決して変更しないので、一つの Model に対して 複数の Simulator を差し替えることができます。

9. Spatial Gillespie Method

9.1. E-Cell4 における空間表現

E-Cell4における空間表現とはどのようなものでしょうか?

from ecell4 import *

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

ODEWorld と GillespieWorld にサイズ 1 の立方体を作成しました。 この場合、重要なのはボリュームのみです。

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'

これは同じ結果を返します。 ボリュームが 同じ 1 であるためです。

これは均一な系では妥当ですが、細胞は均一な系ではありません。 したがって、我々は分子の局在化のために空間表現を考慮する必要があります。

E-Cell4では、いくつかのタイプの空間表現とシミュレーション方法を使用できます。 ここではSpatial Gillespie Methodの例を示します。

9.2. Spatial Gillespie Method

E-Cell4では、Spatial Gillespie Method が meso モジュールに含まれています。

まずは ode でのように run_simulationから始めてみましょう。

%matplotlib inline
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')

定常状態においては、 C の数は以下のように与えられます。

あなたは ode や gillespie とほとんど同じ結果を得るでしょう (ode や gillespie よりも時間はかかるかもしれません)。

これは当然のことです。なぜなら、meso モジュールは、追加の空間パラメータを与えない限り、Gillespieとほぼ同じだからです。

次に、 run_simulationを分解して見てみましょう。

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)

これは MesoscopicWorld と MesoscopicSimulator を除いて普段のものと変わりませんが、いくつかの新しい要素が見えるかと思います。

まず w.bind_to(m) で、Model を World に関連付けました。 以前の基本的な演習では、私たちはこれをしませんでした。

spatial methodでは、Species の属性が必要です。 これを忘れないでください。 その後、 MesoscopicSimulatorを作るためには World だけが必要です。

次に、重要な違いは、MesoscopicWorld の2番目の引数、すなわち Integer3(1,1,1) です。 ODEWorldとGillespieWorldはこの第二引数を持っていません。 これを説明する前に、この引数を変更してシミュレーションを再実行してみましょう。

from ecell4 import *

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

m = get_model()

w = meso.MesoscopicWorld(Real3(1, 1, 1), Integer3(4, 4, 4))  # 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)

先程とは異なるプロットが出ているはずです。 Integer3 の値を大きくすると、より大きく異なっているプロットになるでしょう。

実は、この第2引数は、空間パーティションの数を意味します。

meso は gillespie とほぼ同じですが、meso は空間を立方体に分割します（我々はこの立方体をサブボリュームと呼んでいます）。 そして gillespie がただ1つの均一な閉鎖空間を持つのとは対照的に、各サブボリュームは異なる分子濃度を有します。

したがって、前の例では、我々は 1 の辺を持つ1つの立方体を，辺 0.25 の64(444) 個の立方体に分割したわけです。

そして次に60個の C 分子を World に投げました。したがって、各サブボリュームにはせいぜい1つのspeciesしかないわけです。

9.3. 分子の拡散係数の定義

先程の違いはどこから来ているのでしょうか?

これは meso で space を得たにもかかわらず，分子の拡散係数を考慮していないためです。

拡散係数を設定するには、前に(2. モデルの構築方法) で述べた方法で Species の attribute D を使います。

1. E-Cell4シミュレーションの概要 でお見せしたように，ここでは E-Cell4 の特殊記法を用います。

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

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

get_model()

<ecell4.core.NetworkModel at 0x7ff0add53c48>

with species_attributes(): ステートメントで拡散係数を設定することができます。

ここでは、すべての拡散係数を 1 に設定します。このモデルをもう一度シミュレートしましょう。

今度は、大きな Integer3 の値でも gillespie とほぼ同じ結果を得るはずです（シミュレーションには gillespie よりはるかに時間がかかります）。

分子拡散は問題に対してどのように働いたのでしょうか?

3D空間で自由拡散 (Speciesの拡散係数は )を考えてみましょう。

拡散係数の単位は、

$ :raw-latex:`\mathrm{\mu m}`^2/s $ or $ :raw-latex:`mathrm{nm}`^2/:raw-latex:`mu `s $

のように(空間の)長さの2乗を時間で割ったものです。

時間$ 0

System Message: WARNING/2 (から)

latex exited with error [stdout] This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) (preloaded format=latex) restricted \write18 enabled. entering extended mode (./math.tex LaTeX2e <2016/02/01> Babel <3.9q> and hyphenation patterns for 81 language(s) loaded. (/usr/share/texlive/texmf-dist/tex/latex/base/article.cls Document Class: article 2014/09/29 v1.4h Standard LaTeX document class (/usr/share/texlive/texmf-dist/tex/latex/base/size12.clo)) (/usr/share/texlive/texmf-dist/tex/latex/base/inputenc.sty (/usr/share/texlive/texmf-dist/tex/latex/ucs/utf8x.def)) (/usr/share/texlive/texmf-dist/tex/latex/ucs/ucs.sty (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-global.def)) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsmath.sty For additional information on amsmath, use the `?' option. (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amstext.sty (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsgen.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsbsy.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsopn.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amscls/amsthm.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amssymb.sty (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amsfonts.sty)) (/usr/share/texlive/texmf-dist/tex/latex/tools/bm.sty) (./math.aux) (/usr/share/texlive/texmf-dist/tex/latex/ucs/ucsencs.def) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsa.fd) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsb.fd) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12363 = U+304B, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $か ら$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12425 = U+3089, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $から $ [1] (./math.aux) ) (see the transcript file for additional information) Output written on math.dvi (1 page, 240 bytes). Transcript written on math.log.

t

System Message: WARNING/2 (までの点の距離の平方の平均が)

latex exited with error [stdout] This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) (preloaded format=latex) restricted \write18 enabled. entering extended mode (./math.tex LaTeX2e <2016/02/01> Babel <3.9q> and hyphenation patterns for 81 language(s) loaded. (/usr/share/texlive/texmf-dist/tex/latex/base/article.cls Document Class: article 2014/09/29 v1.4h Standard LaTeX document class (/usr/share/texlive/texmf-dist/tex/latex/base/size12.clo)) (/usr/share/texlive/texmf-dist/tex/latex/base/inputenc.sty (/usr/share/texlive/texmf-dist/tex/latex/ucs/utf8x.def)) (/usr/share/texlive/texmf-dist/tex/latex/ucs/ucs.sty (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-global.def)) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsmath.sty For additional information on amsmath, use the `?' option. (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amstext.sty (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsgen.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsbsy.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsopn.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amscls/amsthm.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amssymb.sty (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amsfonts.sty)) (/usr/share/texlive/texmf-dist/tex/latex/tools/bm.sty) (./math.aux) (/usr/share/texlive/texmf-dist/tex/latex/ucs/ucsencs.def) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsa.fd) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsb.fd) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12414 = U+307E, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $ま での点の距離の平方の平均が$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12391 = U+3067, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $まで の点の距離の平方の平均が$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12398 = U+306E, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $までの 点の距離の平方の平均が$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-112.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-112.def) ! Package ucs Error: Unknown Unicode character 28857 = U+70B9, (ucs) possibly declared in uni-112.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $までの点 の距離の平方の平均が$ Package ucs Warning: Unknown character 12398 = 0x306E appeared again. on input line 12. (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-141.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-141.def) ! Package ucs Error: Unknown Unicode character 36317 = U+8DDD, (ucs) possibly declared in uni-141.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $までの点の距 離の平方の平均が$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-150.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-150.def) ! Package ucs Error: Unknown Unicode character 38626 = U+96E2, (ucs) possibly declared in uni-150.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $までの点の距離 の平方の平均が$ Package ucs Warning: Unknown character 12398 = 0x306E appeared again. on input line 12. (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-94.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-94.def) ! Package ucs Error: Unknown Unicode character 24179 = U+5E73, (ucs) possibly declared in uni-94.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $までの点の距離の平 方の平均が$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-101.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-101.def) ! Package ucs Error: Unknown Unicode character 26041 = U+65B9, (ucs) possibly declared in uni-101.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $までの点の距離の平方 の平均が$ Package ucs Warning: Unknown character 12398 = 0x306E appeared again. on input line 12. Package ucs Warning: Unknown character 24179 = 0x5E73 appeared again. on input line 12. (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-87.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-87.def) ! Package ucs Error: Unknown Unicode character 22343 = U+5747, (ucs) possibly declared in uni-87.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $までの点の距離の平方の平均 が$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12364 = U+304C, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $までの点の距離の平方の平均が $ Overfull \hbox (319.22516pt too wide) in paragraph at lines 12--13 []$[][][][][][][][][][][][][][][][][][][][][][][][][][][][]$ [1] (./math.aux) ) (see the transcript file for additional information) Output written on math.dvi (1 page, 484 bytes). Transcript written on math.log.

6Dt $に等しいことが知られています。

逆に、長さスケール$ l

System Message: WARNING/2 (の空間における時間スケールの平均は約)

latex exited with error [stdout] This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) (preloaded format=latex) restricted \write18 enabled. entering extended mode (./math.tex LaTeX2e <2016/02/01> Babel <3.9q> and hyphenation patterns for 81 language(s) loaded. (/usr/share/texlive/texmf-dist/tex/latex/base/article.cls Document Class: article 2014/09/29 v1.4h Standard LaTeX document class (/usr/share/texlive/texmf-dist/tex/latex/base/size12.clo)) (/usr/share/texlive/texmf-dist/tex/latex/base/inputenc.sty (/usr/share/texlive/texmf-dist/tex/latex/ucs/utf8x.def)) (/usr/share/texlive/texmf-dist/tex/latex/ucs/ucs.sty (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-global.def)) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsmath.sty For additional information on amsmath, use the `?' option. (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amstext.sty (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsgen.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsbsy.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsopn.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amscls/amsthm.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amssymb.sty (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amsfonts.sty)) (/usr/share/texlive/texmf-dist/tex/latex/tools/bm.sty) (./math.aux) (/usr/share/texlive/texmf-dist/tex/latex/ucs/ucsencs.def) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsa.fd) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsb.fd) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12398 = U+306E, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の 空間における時間スケールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-122.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-122.def) ! Package ucs Error: Unknown Unicode character 31354 = U+7A7A, (ucs) possibly declared in uni-122.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空 間における時間スケールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-149.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-149.def) ! Package ucs Error: Unknown Unicode character 38291 = U+9593, (ucs) possibly declared in uni-149.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間 における時間スケールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12395 = U+306B, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間に おける時間スケールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12362 = U+304A, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間にお ける時間スケールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12369 = U+3051, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間におけ る時間スケールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12427 = U+308B, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間における 時間スケールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-102.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-102.def) ! Package ucs Error: Unknown Unicode character 26178 = U+6642, (ucs) possibly declared in uni-102.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間における時 間スケールの平均は約$ Package ucs Warning: Unknown character 38291 = 0x9593 appeared again. on input line 12. (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12473 = U+30B9, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間における時間ス ケールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12465 = U+30B1, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間における時間スケ ールの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12540 = U+30FC, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間における時間スケー ルの平均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12523 = U+30EB, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 $の空間における時間スケール の平均は約$ Package ucs Warning: Unknown character 12398 = 0x306E appeared again. on input line 12. (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-94.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-94.def) ! Package ucs Error: Unknown Unicode character 24179 = U+5E73, (ucs) possibly declared in uni-94.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 ...空間における時間スケールの平 均は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-87.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-87.def) ! Package ucs Error: Unknown Unicode character 22343 = U+5747, (ucs) possibly declared in uni-87.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 ...間における時間スケールの平均 は約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-48.def) ! Package ucs Error: Unknown Unicode character 12399 = U+306F, (ucs) possibly declared in uni-48.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 ...における時間スケールの平均は 約$ (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-125.def) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uninames.dat) (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-125.def) ! Package ucs Error: Unknown Unicode character 32004 = U+7D04, (ucs) possibly declared in uni-125.def. (ucs) Type H to see if it is available with options. See the ucs package documentation for explanation. Type H <return> for immediate help. ... l.12 ...おける時間スケールの平均は約 $ Overfull \hbox (516.82526pt too wide) in paragraph at lines 12--13 []$[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]$ [1] (./math.aux) ) (see the transcript file for additional information) Output written on math.dvi (1 page, 572 bytes). Transcript written on math.log.

l ^ 2 / 6D $です。

上記の場合、各サブボリュームのサイズは0.25であり、拡散係数は1です。 したがって、時間スケールは約0.01秒になります。

Species A と B の分子が同じサブボリュームにある場合、反応に約1.5秒かかるので、ほとんどの場合、拡散は反応よりも速く、同じサブボリュームで解離した場合でも分子は他のサブボリュームに移動します。

小さな では、より小さなサブボリュームの体積 となるので、解離後の反応速度は速くなり、サブボリューム間の拡散時間と遷移時間も小さくなります。

9.4. 分子の局在

我々は add_molecules 関数を使って、ode や gillespie と同じように World に分子を追加しました。 一方、 MesoscopicWorld では、空間表現に従って分子を入れることができます。

from ecell4 import *

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

MesoscopicWorld では、第3引数 Integer3 を add_molecules に与えることでサブボリュームと分子の位置を設定できます。

上記の例では、Speciesタイプ A は空間全体に広がっていますが、Speciesタイプ B は、ボリュームの中心のサブボリュームにのみ位置するようになっています。

またこれらの一情報をチェックするには、座標と共に num_molecules関数を使います。

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

120
0
120

さらに、Jupyter Notebook 環境があれば、 ecell4.viz モジュールを使って分子の位置を可視化することができます。

# viz.plot_world(w, radius=0.01)
viz.plot_world(w, interactive=False)

viz.plot_world 関数は、 World を与えることによって、Jupyter Notebookのセルで分子の位置を可視化します。

分子のサイズは radius で設定できます。

それでは 分子の局在を World に設定し、これをシミュレートしてみましょう。

上記の例では、拡散係数 1 とWorldサイズ 1 を設定しているので、かき混ぜる時間は 10 秒で十分です。

シミュレーションの後で viz.plot_world を再度呼び出して結果を確認してください。

9.5. 分子の初期位置とその反応

これは、分子の局在が反応にどのように影響するかを調べる極端な例です。

%matplotlib inline
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)

このモデルは単純な結合反応のみからなります。 World はx軸に長い直方体であり、分子は中心から離れた位置にあります。

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)
viz.plot_world(w, interactive=False)

別の注意点として、Integer3(0, 0, 0) または Integer3(9, 0, 0) を設定しない理由があります。

E-Cell4では、基本的にすべてに周期的境界条件を採用しています。

したがって、前述の2つのサブボリュームは実際には隣接しています。

期待される位置を実現した後、 MesoscopicSimulator でシミュレートしてみましょう。

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)

# viz.plot_world(w, radius=0.025)
viz.plot_world(w, interactive=False)

初期座標の効果を確認するには、分子を meso で均一に配置するか、 gillespie でシミュレートすることをお勧めします。

w = meso.MesoscopicWorld(Real3(10, 1, 1), Integer3(10, 1, 1))
w.bind_to(m)
w.add_molecules(Species('A'), 1200)
w.add_molecules(Species('B'), 1200)

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

実線はバイアスをかけたケースで、破線はかけていないものになります。

バイアスがかかっている場合の反応は明らかに遅いです。

また、時系列の形状も実線と破線の間で異なっていることがわかります。

これは、最初の分離のために分子AおよびBが衝突するのにある程度の時間がかかるからです。

実際には、A と B の間の最初の距離（約4）を移動するには 秒かかります。

10. Spatiocyteを用いた一分子粒度のシミュレーション

これまでのチュートリアルで E-Cell4 の空間表現の例を示しました。 次に、「一分子粒度」と呼ばれる、より詳細な空間表現を用いてモデルをシミュレートしてみましょう。

!pip install ecell

Collecting ecell
  Downloading ecell-4.1.2-cp36-cp36m-manylinux1_x86_64.whl (45.5MB)
    100% |################################| 45.5MB 31kB/s  eta 0:00:01
Installing collected packages: ecell
Successfully installed ecell-4.1.2

%matplotlib inline
from ecell4 import *

10.1. 格子法のSpatiocyte

空間的なGillespie法では、Space をより小さな Space に分割し、サブボリューム内の分子を拡散させます。 しかし、我々は各サブボリュームにおける分子を分子の数として扱い、分子の位置は定められていませんでした。

言い換えると、空間的Gillespie法の空間分解能はサブボリューム の面に等しくなります。 この分解能を改善するには、のサイズを小さくする必要があります。 しかし、この方法では、は分子の直径の(少なくとも)10倍より大きくなければなりません。

どうすれば分子のサイズに対して空間分解能を向上させることができるでしょうか？ 答えは、1分子粒度のシミュレーションになります。 この方法は、分子の数ではなく各分子の空間的な反応拡散を用いて分子をシミュレートします。

E-Cell4には複数の一分子粒度の手法があり、ここでは格子法のSpatiocyteについて説明します。 Spatiocyteは、各分子を剛体球として扱い、分子を六方最密格子に拡散させます。

Spatiocyteは、各分子のIDとその分子の一分子粒度での位置を保持しています。 より高い空間分解能のために、SpatiocyteはSpatical Gillespieより100倍小さい時間ステップを持っています。 なぜならその拡散の時間スケールが距離の2乗とともに増加するためです。

それではSpatiocyteを試してみましょう。

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

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

m = get_model()

w = spatiocyte.SpatiocyteWorld(Real3(1, 1, 1), 0.005)  # The second argument is 'voxel_radius'.
w.bind_to(m)
w.add_molecules(Species('C'), 60)

sim = spatiocyte.SpatiocyteSimulator(w)
obs = FixedIntervalNumberObserver(0.1, ('A', 'B', 'C'))
sim.run(10, obs)

SpatiocyteWorldの第二引数には明確な違いがあります。 これをvoxel radiusといいます。 Spatiocyteでは、空間を分子の大きさで割って分子の位置を定義し、この空間の最小単位を Voxelと呼んでいます。

ほとんどの場合、分子のサイズはvoxel radiusに適しているでしょう。 この例では、一辺が の空間内の分子の半径として5 を設定しました。 シミュレーションに時間がかかりますが、結果はODEまたはGillespieと同じです。

viz.plot_number_observer(obs)

10.2. 一分子の拡散運動

次にその分解能を確認するために一分子粒度の拡散をシミュレートしてみましょう。

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

m = get_model()

w = spatiocyte.SpatiocyteWorld(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メソッドはSpatiocyteWorldの座標に粒子を置こうとします。 またnew_particleメソッドは粒子のParticleID(pid)、Particleオブジェクト(p)、そしてブール値(suc)を返します。 new_particleが粒子を置くことに成功した場合、sucはTrueになります。

粒子がすでに座標に配置されている場合は、粒子を置くことはできず、 sucはFalseになり失敗します。

pは粒子の位置、species type、半径、拡散係数を含んでいます。 pの情報はそのID pid を使って調べることが可能です。

それではまず pをチェックしてみましょう。

pid, p = w.get_particle(pid)
print(p.species().serial())  # must print: A
print(p.radius(), p.D())  # must print: (0.005, 1.0)
print(tuple(p.position()))  # must print: (0.49806291436591293, 0.49652123150307814, 0.5)

A
(0.005, 1.0)
(0.49806291436591293, 0.49652123150307814, 0.5)

get_particleメソッドは、粒子のIDを受け取り、IDと粒子を返します（もちろんそのIDは与えられるIDと同じです）。 またposition()メソッドで粒子の座標を Real3として調べることができます。 座標を直接読み取ることは難しいですが、ここではタプルに変換してからプリントしました。

上記の結果からわかるように、タプルのコーディネイトは Real3として与えられた元の位置とわずかに異なります。 これは、Spatiocyteが分子を格子上にのみ置くことができるからです。 SpatiocyteWorldは、分子を引数Real3の最も近い格子の中心に置きます。

viz.plot_worldメソッドを使って分子の座標を可視化することができます。この例ではWorldの中心にある分子が確認できます。

viz.plot_world(w, interactive=False)
# viz.plot_world(w)

また、 FixedIntervalTrajectoryObserverを使って、分子拡散プロセスの軌道を追跡することができます。

sim = spatiocyte.SpatiocyteSimulator(w)
obs = FixedIntervalTrajectoryObserver(0.002, (pid,))
sim.run(1, obs)
viz.plot_trajectory(obs, interactive=False)
# viz.plot_trajectory(obs)

ここでは viz.plot_trajectoryメソッドで軌跡を可視化しましたが、dataメソッドでReal3のリストとして結果を取得することも可能です。

print(len(obs.data()))  # => 1
print(len(obs.data()[0]))  # => 501

1
501

dataメソッドはネストされたリストを返します。 最初のインデックスは、粒子のインデックスを意味します。 2番目のインデックスは、 Real3のインデックスを意味します。 この場合、我々は1つだけパーティクルを投げたので、最初の結果は「1」であり、次の「501」は、1つのパーティクル（初期座標と1 / 0.002 = 500タイムポイントの座標）の時系列座標を意味します。

list_particles_exactメソッドとSpeciesメソッドを使って、一括してパーティクルを取得することもできます。

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

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

(u'A', (0.7103520254071217, 0.8256108849411649, 0.88))
(u'A', (0.5062278801751902, 0.034641016151377546, 0.77))
(u'A', (0.4082482904638631, 0.7188010851410841, 0.985))
(u'A', (0.30210373494325865, 0.30599564267050167, 0.68))
(u'A', (0.28577380332470415, 0.26269437248127975, 1.0050000000000001))
(u'A', (0.563382640840131, 0.0202072594216369, 0.545))

list_particles_exactメソッドは覚えておいてください。 このメソッドは、他のWorldに対してもadd_moleculesメソッド同様に使えます。

別の注意点として、Spatiocyteで一つの分子を調査する適切な方法は list_voxels_exactであり、座標はボクセルのインデックス（Real3ではありません）で記述されます。

10.3 拡散係数と二次反応

私たちが扱っているモデルには二次反応が含まれています。 この二次反応とSpatiocyteの拡散係数との関係を見てみましょう。

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

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

m = get_model()

w = spatiocyte.SpatiocyteWorld(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 = spatiocyte.SpatiocyteSimulator(w)
sim.run(1.0, obs)

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(m, odew)
odesim.run(1.0, odeobs)

viz.plot_number_observer(obs, "-", odeobs, "--")

以前よりもより速い速度定数を使用しましたが、結果は同じです。 しかし、ODEシミュレーションと比べると、それらの違いを見つけることができます（実線はspatiocyte、破線はode）。 Spatiocyteが間違っているのでしょうか？ （いいえ，間違っているわけではありません。） 実はSpatiocyteの反応速度は速くなりません。一方ODEの反応速度は無限に速くなり得ます。

これは、ODEソルバーと一分子シミュレーション法の反応速度定数の定義の違いが原因です。 前者は」macroscopic」なまたは」effective」な反応速度定数と呼ばれ、後者は」microscopic」なまたは」intrinsic」な反応速度定数と呼ばれます。

「macroscopic」な速度は分子が混ざっている状態における反応速度を表し、一方、」microscopic」な速度は分子衝突における反応性を表します。 そのため 「microscopic」な視点では、最初に分子が反応するためにまず必要なのは衝突です。 しかし、Spatiocyteでは、この」microscopic」な速度を速くすると、実際の反応速度を衝突速度より速くすることはできません。 これは」拡散律速」条件と呼ばれます。 これは、偏って配置された分子が反応するのに時間を要することに似ています。

この」macroscopic」な速度定数 と3次元空間における」microscopic」な速度定数 との間には下記の関係があることが知られている。

ここで は衝突における2分子の半径の合計に、 は拡散係数の合計になります。

上記の Jupyter Notebook のセルの場合, はおおよそ 0.25 で 「microscopic」 な速度定数は 1.0 になります。 そのため 「macroscopic」 な速度定数は約 0.2 になります。

(しかし、Spatiocyteの設定を指定しない限り、二次反応速度は より遅くなければならず、解離定数 も となります。)

一分子シミュレーション法は、よく混ざっている系(すなわち拡散係数が無限大)が想定されたODEやGillespie法とは違って、分子の」拡散」と」反応」を正確に分離することができます。

しかし、」microscopic」な速度定数 が十分に小さい場合、」macroscopic」な速度定数は」microscopic」な速度定数とほとんど同じになります(これは」reaction late-limit」と呼ばれています)。

10.4. Spatiocyte における 「構造(Structure)」

次に、細胞膜のような構造(Structure)を作り出す方法について説明します。

E-Cell4全体としてのStructure機能はまだ開発中ですが、(E-Cell4が持つシミュレータの一つである)SpatiocyteにおいてはStructureがある程度サポートできています。

その一例として、球のStructureを見てみましょう。

球の内部の分子の拡散を(その内側に)制限するためには、まずその球を作成します。

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

m = get_model()

w = spatiocyte.SpatiocyteWorld(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

539805

ここで World の状態を可視化してみましょう。

viz.plot_world(w, interactive=False)
# viz.plot_world(w)

Sphere クラスの第1引数は球の中心であり、第2引数は半径です。 次に、 C という名前の Species を作成して球の内部に追加しました。

SpatiocyteのStructureは、空間を Voxel で満たすことによって記述されます。

上記の例では、球体の Voxel は C という名前の Species で占められています。

上記のように viz.plot_world を使ってそれらの分布を見ることができます。 (しかし、Speciesの数はそのすべてを可視化するには大きすぎます。 そのため、私たちはその一部だけをプロットしています。しかし実際には球体は完全にSpeciesで占められています。)

次に、この球の内部を移動するSpeciesを作成します。

そのためにSpeciesにlocationの属性を与えます。 その後は、add_molecules関数を使ってWorldに分子を投げ込むだけです。

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

viz.plot_world(w, species_list=('A',), interactive=False)  # visualize A-molecules only
# viz.plot_world(w, species_list=('A',))  # visualize A-molecules only

ここで、私達は Species A の軌道を Species CのStructureに制限したので、add_moleculesはそのように働きます。 注意として、 add_moleculesの前にStructureを作成する必要があることがあります。

そしてFixedIntervalTrajectoryObserverを使って拡散領域の制限をチェックすることができます。

pid_list = [pid for pid, p in w.list_particles(Species('A'))[: 10]]
obs = FixedIntervalTrajectoryObserver(1e-3, pid_list)
sim = spatiocyte.SpatiocyteSimulator(w)
sim.run(1, obs)
viz.plot_trajectory(obs, interactive=False)
# viz.plot_trajectory(obs)

pid_listはA分子の ParticleIDの最初の10個のリストです。 軌道はこの10Speciesによって着色されています。 確かに、軌道は球体内に制限されています。

10.5 構造(Structure)と反応

最後に、Structure間の分子のトランスロケーションについて説明します。

location 属性のない Species は、どのStructureのメンバーでもありません。 上記の例では、Species A で location 属性を記述しないと、球の外側に A が置かれます。

次に、平面構造を作成しましょう。 面を作成するには、3つの Real3 を使用する必要があります。 それは原点(origin)と，軸ベクトルの2点(unit0、unit1)です。 プログラムでは ps = PlanarSurface(origin、unit0、unit1) となります。

面 ps に Species A があり、通常の Species B がある状態を想像してみてください。

with species_attributes():
    A | {'D': '0.1', 'location': 'M'}
    B | {'D': '1'}

m  = get_model()

w = spatiocyte.SpatiocyteWorld(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=('B', 'M'), interactive=False)
# viz.plot_world(w, species_list=('B', 'M'))

見るのは難しいかもしれませんが、実際にはSpecies B は面上でないところにしか置かれていません。

それではどうすれば、このSpecies B を面 M に吸収させ、Species A を合成することができるでしょうか?

with species_attributes():
    A | {'D': '0.1', 'location': 'M'}
    B | {'D': '1'}

with reaction_rules():
    B + M == A | (1.0, 1.5)

m = get_model()

これは、 Species B が Structure M に衝突したとき、 Aになることを意味します。

一方、Species A はStructureから解離し、逆反応の方向としては M そして B になります。

これでStructureを用いたモデルをシミュレーションできるようになりました。

w.bind_to(m)

sim = spatiocyte.SpatiocyteSimulator(w)
obs = NumberObserver(('A', 'B'))
sim.run(2, obs)

viz.plot_number_observer(obs)
viz.plot_world(w, species_list=('A', 'B'), interactive=False)
# viz.plot_world(w, species_list=('A', 'B'))

Structureからの解離反応では、そのStructureを書くことが省略できます。 したがって、A > B は、上の A > B + M と同じ意味になります。

しかし, 結合反応においては これと同様の省略はできません。

周りに M がない状態で B から A を作ることは不可能だからです。 対照的に、Species A は、球 M 上のどこでも Species B を生成することができます。

一次反応は、Structureの存るところまたは無いところのいずれかで起こります。

しかし、結合反応の場合には、二次反応は一次反応に変わり、左辺の M を無視すれば速度定数の意味も変わります。

API