ここでは、World クラスの基礎を説明します。
E-Cell4では、spatiocyte.SpatiocyteWorld, egfrd.EGFRDWorld, bd.BDWorld, meso.MesoscopicWorld, gillespie.GillespieWorld, そして ode.ODEWorldの6種類のWorldクラスがサポートされています。
ほとんどのソフトウェアでは、初期状態は Modelの一部であると考えられています。
しかし、E-Cell4では、初期条件は Modelとは別にWorldとして設定する必要があります。
World には、現在時刻、分子数、分子座標、構造、乱数発生器などの状態に関する情報が格納されています。
一方、Model には分子間の相互作用のタイプと分子の共通の性質が含まれています。
In [1]:
from ecell4_base.core import *
In [2]:
from ecell4_base import *
Worldクラスは、アルゴリズムに固有のパラメータを決定するコンストラクタ内に異なるセットの引数を受け取ります。
しかし、少なくとも、すべてのWorldクラスは、edge_lengthsという名前のサイズでしかインスタンス化できません。
edge_lengthsのタイプはReal3であり、これはRealsのトリプレットを表します。
E-Cell4では、すべての3次元位置がReal3として扱われます。
In [3]:
edge_lengths = Real3(1, 2, 3)
w1 = gillespie.World(edge_lengths)
w2 = ode.World(edge_lengths)
w3 = spatiocyte.World(edge_lengths)
w4 = bd.World(edge_lengths)
w5 = meso.World(edge_lengths)
w6 = egfrd.World(edge_lengths)
Worldにはサイズと体積のゲッターメソッドがあります。
In [4]:
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())
次に、分子をWorldに加えましょう。
ここでは、分子の形状を伝えるために、「半径」と「拡散係数」の属性付きのSpeciesを与える必要があります。
下記の例では「0.0025」が半径に，「1」が拡散係数に相当します。
分子の位置は、必要に応じてWorldでランダムに決定されます。
add_molecules関数における10は追加する分子の数を表しています。
In [5]:
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で設定した半径と拡散係数を再び書く必要はありません（変更したくない限り）。
In [6]:
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も利用できます。
In [7]:
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)
In [8]:
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")))
num_molecules_exactと異なりnum_moleculesは与えられたSpeciesと一致するすべての数をルールベースの方法で返します。
WorldのすべてのSpeciesが分子間相互作用を持たない場合、num_moleculesはnum_molecules_exactと同じです。
In [9]:
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")))
World はシミュレーションの時間を保持しています。
In [10]:
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())
最後に、HDF5ファイルにWorldの状態を保存したりロードすることができます。
In [11]:
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
In [12]:
w1 = gillespie.World()
w2 = ode.World()
w3 = spatiocyte.World()
w4 = bd.World()
w5 = meso.World()
w6 = egfrd.World()
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")))
In [13]:
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
すべてのWorldクラスでは、コンストラクターの唯一の引数としてHDF5ファイルパスも使用できます。
In [14]:
print(gillespie.World("gillespie.h5").t())
print(ode.World("ode.h5").t())
print(spatiocyte.World("spatiocyte.h5").t())
print(bd.World("bd.h5").t())
print(meso.World("meso.h5").t())
print(egfrd.World("egfrd.h5").t())
World には、分子の座標にアクセスする共通の機能もあります。
In [15]:
w1 = gillespie.World()
w2 = ode.World()
w3 = spatiocyte.World()
w4 = bd.World()
w5 = meso.World()
w6 = egfrd.World()
まず、new_particleを使って特定の位置に分子を置くことができます。
In [16]:
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 = 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は作成されたパーティクルとそれが成功したかどうかを返します。
分子の表現における分解能は異なっています。
たとえば、gillespie.World には分子の座標に関する情報はほとんどありません。
それ故 gillespie.World は与えられた位置を無視し、分子の数をカウントアップすることのみ行います。
ParticleID は lot と serial と名付けられた Integer のペアになります。
In [17]:
print(pid6.lot(), pid6.serial())
print(pid6 == ParticleID((0, 1)))
パーティクルシミュレータ、すなわちspatiocyte、bdおよびegfrdは、idでパーティクルにアクセスするためのインターフェイスを提供しています。
has_particle は、指定された ParticleID に対してパーティクルが存在するかどうかを返します。
In [18]:
# 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
存在するかをチェックした後で、次のように get_particle でパーティクルを取得できます。
In [19]:
# 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(拡散係数) からなります。
In [20]:
# 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())
Spatiocyteの場合、粒子の位置は与えられた位置に最も近いボクセルの中心に自動的に丸められます。
パーティクルの位置を移動させることもできます。
update_particle は指定された ParticleID で指定されたパーティクルを指定された Particle に置き換え、Falseを返します。
ParticleID に対応するパーティクルが見つからない場合は、新しいパーティクルを作成してTrueを返します。
異なるタイプの Species でパーティクルを与えると、パーティクルの Species も変更されます。
In [21]:
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))
list_particles と list_particles_exact は、World の ParticleID と Particle のペアのリストを返します。
In [22]:
print(w1.list_particles_exact(sp1))
# print(w2.list_particles_exact(sp1)) # ode.World 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))
remove_particleを使用して特定のパーティクルを削除することもできます。
In [23]:
# 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
In [24]:
w = spatiocyte.World(Real3(1, 2, 3), voxel_radius=0.01)
w.bind_to(m)
Voxel と呼ばれる単一の格子のサイズは、voxel_radius メソッドによって取得できます。
spatiocyte.World には、行数、列数、およびレイヤ数を取得するメソッドがあります。
これらのサイズは、World構築時に与えられた edge_lengths に基づいて自動的に計算されます。
In [25]:
print(w.voxel_radius()) # => 0.01
print(tuple(w.shape())) # => (64, 152, 118)
# print(w.col_size(), w.row_size(), w.layer_size()) # => (64, 152, 118)
print(w.size()) # => 1147904 = 64 * 152 * 118
格子ベースの空間内の位置は、グローバル座標と呼ばれるInteger3 (列、行および層)として扱われます。
このようにして、spatiocyte.World は、Real3を格子ベースの座標に変換する機能を提供します。
In [26]:
# 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)
spatiocyte.Worldでは、グローバル座標は単一の整数に変換されます。 それは単に座標と呼ばれています。
グローバル座標で同じ方法で座標を扱うこともできます。
In [27]:
# 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は新しいボクセルを座標で作成する方法を提供します。
In [28]:
# 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)
Voxel は、species、座標、半径およびD からなります。
In [29]:
# print(v.species().serial(), v.coordinate(), v.radius(), v.D()) # => (u'A', 278033, 0.0025, 1.0)
もちろん、get_particleとlist_particles_exact同様にget_voxelとlist_voxels_exactを使ってVoxelとVoxelのリストを取得することができます。
In [30]:
# print(w.num_voxels_exact(Species("A")))
# print(w.list_voxels_exact(Species("A")))
# print(w.get_voxel(pid))
update_particleに対応するupdate_voxelでボクセルを移動および更新することができます。
In [31]:
# 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
最後に、remove_voxelは、remove_particleのようにVoxelを削除します。
In [32]:
# print(w.has_voxel(pid)) # => True
# w.remove_voxel(pid)
# print(w.has_voxel(pid)) # => False
In [33]:
w1 = gillespie.World()
w2 = ode.World()
w3 = spatiocyte.World()
w4 = bd.World()
w5 = meso.World()
w6 = egfrd.World()
Shapeオブジェクトを使用すると、分子の初期位置をWorldの一部に限定できます。
以下の場合、60個の分子が与えられた球の内部に配置されます。
ここに置かれた分子の拡散はShapeに制限されていません。
このShapeは、初期化専用です。
In [34]:
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を定義します。
In [35]:
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にバインドできます。
In [36]:
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とはサブボリュームのリストを意味します。 したがって、構造は他の粒子の侵入を避けることはできません。
Random Number GeneratorもWorldの一部です。
ode.World以外のすべてのWorldでは、Random Number Generatorが格納されており、シミュレーションにランダムな値が必要なときに更新されます。
E-Cell4では、Random Number GeneratorとしてGSLRandomNumberGeneratorクラスが1つだけ実装されています。
In [37]:
rng1 = GSLRandomNumberGenerator()
print([rng1.uniform_int(1, 6) for _ in range(20)])
引数を指定しない場合、Random Number Generatorは常に0のシードで初期化されます。
In [38]:
rng2 = GSLRandomNumberGenerator()
print([rng2.uniform_int(1, 6) for _ in range(20)]) # => same as above
次のように、シードを整数で初期化することができます。
In [39]:
rng2 = GSLRandomNumberGenerator()
rng2.seed(15)
print([rng2.uniform_int(1, 6) for _ in range(20)])
入力なしでシード関数を呼び出すと、シードは現在時刻の情報から生み出されます。
In [40]:
rng2 = GSLRandomNumberGenerator()
rng2.seed()
print([rng2.uniform_int(1, 6) for _ in range(20)])
GSLRandomNumberGeneratorは、乱数を得るためのいくつかの方法を提供します。
In [41]:
print(rng1.uniform(0.0, 1.0))
print(rng1.uniform_int(0, 100))
print(rng1.gaussian(1.0))
Worldは、構築時にRandom Number Generatorを受け入れます。
デフォルトでは、GSLRandomNumberGenerator()が使用されます。
そのため、Generatorを与えないと、シミュレーションの振る舞いは常に同じです（これをdeterminisiticと呼びます）。
In [42]:
rng = GSLRandomNumberGenerator()
rng.seed()
w1 = gillespie.World(Real3(1, 1, 1), rng=rng)
World中のGSLRandomNumberGeneratorへはrng関数を使ってアクセスできます。
In [43]:
print(w1.rng().uniform(0.0, 1.0))
rng()は、GSLRandomNumberGeneratorへの共有ポインタを返します。
したがって、上記の例では、rngとw1.rng()は全く同じことを指しています。