!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()は全く同じことを指しています。