{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 3. How to Setup the Initial Condition\n", "\n", "Here, we explain the basics of `World` classes. In E-Cell4, six types of World classes are offically supported now: `spatiocyte.SpatiocyteWorld`, `egfrd.EGFRDWorld`, `bd.BDWorld`, `meso.MesoscopicWorld`, `gillespie.GillespieWorld`, and `ode.ODEWorld`.\n", "\n", "In the most of softwares, the initial condition is supposed to be a part of `Model`. But, in E-Cell4, the initial condition must be set up as `World` separately from `Model`. `World` stores an information about the state at a time point, such as a current time, the number of molecules, coordinate of molecules, structures, and random number generator. Meanwhile, `Model` contains the type of interactions between molecules and the common properties of molecules." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import ecell4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.1. Common APIs of World\n", "\n", "Even though `World` describes the spatial representation specific to the corresponding algorithm, it has compatible APIs. In this section, the common interfaces of six `World` classes are introduced." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ecell4.core import *\n", "from ecell4.gillespie import GillespieWorld\n", "from ecell4.ode import ODEWorld\n", "from ecell4.spatiocyte import SpatiocyteWorld\n", "from ecell4.bd import BDWorld\n", "from ecell4.meso import MesoscopicWorld\n", "from ecell4.egfrd import EGFRDWorld" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`World` classes accept different sets of arguments in the constructor, which determine the parameters specific to the algorithm. However, at least, all `World` classes can be instantiated only with their size, named `edge_lengths`. The type of `edge_lengths` is `Real3`, which represents a triplet of `Real`s. In E-Cell4, all 3-dimensional positions are treated as a `Real3`. See also [8. More about 1. Brief Tour of E-Cell4 Simulations](8. More about 1. Brief Tour of E-Cell4 Simulations.ipynb)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "edge_lengths = Real3(1, 2, 3)\n", "w1 = GillespieWorld(edge_lengths)\n", "w2 = ODEWorld(edge_lengths)\n", "w3 = SpatiocyteWorld(edge_lengths)\n", "w4 = BDWorld(edge_lengths)\n", "w5 = MesoscopicWorld(edge_lengths)\n", "w6 = EGFRDWorld(edge_lengths)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`World` has getter methods for the size and volume." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "((1.0, 2.0, 3.0), 6.0)\n", "((1.0, 2.0, 3.0), 6.0)\n", "((1.0, 2.0, 3.0), 6.0)\n", "((1.0, 2.0, 3.0), 6.0)\n", "((1.0, 2.0, 3.0), 6.0)\n", "((1.0, 2.0, 3.0), 6.0)\n" ] } ], "source": [ "print(tuple(w1.edge_lengths()), w1.volume())\n", "print(tuple(w2.edge_lengths()), w2.volume())\n", "print(tuple(w3.edge_lengths()), w3.volume())\n", "print(tuple(w4.edge_lengths()), w4.volume())\n", "print(tuple(w5.edge_lengths()), w5.volume())\n", "print(tuple(w6.edge_lengths()), w6.volume())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's add molecules into the World. Here, you must give `Species` attributed with \"radius\" and \"D\" for `EGFRDWorld`, `BDWorld` or `SpatiocyteWorld` to tell the shape of molecules. Positions of the molecules are randomly determined in the `World` if needed." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sp1 = Species(\"A\", \"0.0025\", \"1\")\n", "w1.add_molecules(sp1, 10)\n", "w2.add_molecules(sp1, 10)\n", "w3.add_molecules(sp1, 10)\n", "w4.add_molecules(sp1, 10)\n", "w5.add_molecules(sp1, 10)\n", "w6.add_molecules(sp1, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once binding a `NetworkModel` to the `World`, you don't need to give attributes explicitly. The `World` will ask attributes to the bound `NetworkModel`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m = NetworkModel()\n", "m.add_species_attribute(Species(\"A\", \"0.0025\", \"1\"))\n", "m.add_species_attribute(Species(\"B\", \"0.0025\", \"1\"))\n", "\n", "w1.bind_to(m)\n", "w2.bind_to(m)\n", "w3.bind_to(m)\n", "w4.bind_to(m)\n", "w5.bind_to(m)\n", "w6.bind_to(m)\n", "w1.add_molecules(Species(\"B\"), 20)\n", "w2.add_molecules(Species(\"B\"), 20)\n", "w3.add_molecules(Species(\"B\"), 20)\n", "w4.add_molecules(Species(\"B\"), 20)\n", "w5.add_molecules(Species(\"B\"), 20)\n", "w6.add_molecules(Species(\"B\"), 20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, `remove_molecules` and `num_molecules_exact` are also available." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "w1.remove_molecules(Species(\"B\"), 5)\n", "w2.remove_molecules(Species(\"B\"), 5)\n", "w3.remove_molecules(Species(\"B\"), 5)\n", "w4.remove_molecules(Species(\"B\"), 5)\n", "w5.remove_molecules(Species(\"B\"), 5)\n", "w6.remove_molecules(Species(\"B\"), 5)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(10, 10, 10, 10, 10, 10)\n", "(15, 15, 15, 15, 15, 15)\n" ] } ], "source": [ "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\")))\n", "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\")))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`num_molecules` also count the number of molecules, but returns all the number matched with the given `Species` in the rule-based way. When all `Species` in the `World` has no site and bond, `num_molecules` is almost same with `num_molecules_exact`." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(10, 10, 10, 10, 10, 10)\n", "(15, 15, 15, 15, 15, 15)\n" ] } ], "source": [ "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\")))\n", "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\")))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`World` also owns a simulation time." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)\n", "(1.0, 1.0, 1.0, 1.0, 1.0, 1.0)\n" ] } ], "source": [ "print(w1.t(), w2.t(), w3.t(), w4.t(), w5.t(), w6.t())\n", "w1.set_t(1.0)\n", "w2.set_t(1.0)\n", "w3.set_t(1.0)\n", "w4.set_t(1.0)\n", "w5.set_t(1.0)\n", "w6.set_t(1.0)\n", "print(w1.t(), w2.t(), w3.t(), w4.t(), w5.t(), w6.t())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, you can `save` and `load` the state of a `World` into/from a HDF5 file." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "w1.save(\"gillespie.h5\")\n", "w2.save(\"ode.h5\")\n", "w3.save(\"spatiocyte.h5\")\n", "w4.save(\"bd.h5\")\n", "w5.save(\"meso.h5\")\n", "w6.save(\"egfrd.h5\")\n", "del w1, w2, w3, w4, w5, w6" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(0.0, (1.0, 1.0, 1.0), 1.0, 0, 0)\n", "(0.0, (1.0, 1.0, 1.0), 1.0, 0, 0)\n", "(0.0, (1.0, 1.0, 1.0), 1.0, 0, 0)\n", "(0.0, (1.0, 1.0, 1.0), 1.0, 0, 0)\n", "(0.0, (1.0, 1.0, 1.0), 1.0, 0, 0)\n", "(0.0, (1.0, 1.0, 1.0), 1.0, 0, 0)\n" ] } ], "source": [ "w1 = GillespieWorld()\n", "w2 = ODEWorld()\n", "w3 = SpatiocyteWorld()\n", "w4 = BDWorld()\n", "w5 = MesoscopicWorld()\n", "w6 = EGFRDWorld()\n", "print(w1.t(), tuple(w1.edge_lengths()), w1.volume(), w1.num_molecules(Species(\"A\")), w1.num_molecules(Species(\"B\")))\n", "print(w2.t(), tuple(w2.edge_lengths()), w2.volume(), w2.num_molecules(Species(\"A\")), w2.num_molecules(Species(\"B\")))\n", "print(w3.t(), tuple(w3.edge_lengths()), w3.volume(), w3.num_molecules(Species(\"A\")), w3.num_molecules(Species(\"B\")))\n", "print(w4.t(), tuple(w4.edge_lengths()), w4.volume(), w4.num_molecules(Species(\"A\")), w4.num_molecules(Species(\"B\")))\n", "print(w5.t(), tuple(w5.edge_lengths()), w5.volume(), w5.num_molecules(Species(\"A\")), w5.num_molecules(Species(\"B\")))\n", "print(w6.t(), tuple(w6.edge_lengths()), w6.volume(), w6.num_molecules(Species(\"A\")), w6.num_molecules(Species(\"B\")))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1.0, (1.0, 2.0, 3.0), 6.0, 10, 15)\n", "(1.0, (1.0, 2.0, 3.0), 6.0, 10, 15)\n", "(1.0, (1.0, 2.0, 3.0), 6.0, 10, 15)\n", "(1.0, (1.0, 2.0, 3.0), 6.0, 10, 15)\n", "(1.0, (1.0, 2.0, 3.0), 6.0, 10, 15)\n", "(1.0, (1.0, 2.0, 3.0), 6.0, 10, 15)\n" ] } ], "source": [ "w1.load(\"gillespie.h5\")\n", "w2.load(\"ode.h5\")\n", "w3.load(\"spatiocyte.h5\")\n", "w4.load(\"bd.h5\")\n", "w5.load(\"meso.h5\")\n", "w6.load(\"egfrd.h5\")\n", "print(w1.t(), tuple(w1.edge_lengths()), w1.volume(), w1.num_molecules(Species(\"A\")), w1.num_molecules(Species(\"B\")))\n", "print(w2.t(), tuple(w2.edge_lengths()), w2.volume(), w2.num_molecules(Species(\"A\")), w2.num_molecules(Species(\"B\")))\n", "print(w3.t(), tuple(w3.edge_lengths()), w3.volume(), w3.num_molecules(Species(\"A\")), w3.num_molecules(Species(\"B\")))\n", "print(w4.t(), tuple(w4.edge_lengths()), w4.volume(), w4.num_molecules(Species(\"A\")), w4.num_molecules(Species(\"B\")))\n", "print(w5.t(), tuple(w5.edge_lengths()), w5.volume(), w5.num_molecules(Species(\"A\")), w5.num_molecules(Species(\"B\")))\n", "print(w6.t(), tuple(w6.edge_lengths()), w6.volume(), w6.num_molecules(Species(\"A\")), w6.num_molecules(Species(\"B\")))\n", "del w1, w2, w3, w4, w5, w6" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All the `World` classes also accept a HDF5 file name as an unique argument of the constructor." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0\n", "1.0\n", "1.0\n", "1.0\n", "1.0\n", "1.0\n" ] } ], "source": [ "print(GillespieWorld(\"gillespie.h5\").t())\n", "print(ODEWorld(\"ode.h5\").t())\n", "print(SpatiocyteWorld(\"spatiocyte.h5\").t())\n", "print(BDWorld(\"bd.h5\").t())\n", "print(MesoscopicWorld(\"meso.h5\").t())\n", "print(EGFRDWorld(\"egfrd.h5\").t())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.2. How to Get Positions of Molecules\n", "\n", "`World` has the common functions to access coordinates of molecules too." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "w1 = GillespieWorld()\n", "w2 = ODEWorld()\n", "w3 = SpatiocyteWorld()\n", "w4 = BDWorld()\n", "w5 = MesoscopicWorld()\n", "w6 = EGFRDWorld()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, you can place a molecule at the certain position with `new_particle`." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sp1 = Species(\"A\", \"0.0025\", \"1\")\n", "pos = Real3(0.5, 0.5, 0.5)\n", "(pid1, p1), suc1 = w1.new_particle(sp1, pos)\n", "(pid2, p2), suc2 = w2.new_particle(sp1, pos)\n", "(pid3, p3), suc3 = w3.new_particle(sp1, pos)\n", "(pid4, p4), suc4 = w4.new_particle(sp1, pos)\n", "(pid5, p5), suc5 = w5.new_particle(sp1, pos)\n", "(pid6, p6), suc6 = w6.new_particle(sp1, pos)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`new_particle` returns a particle created and whether it's succeeded or not. However the resolution in representation of molecules differs. For example, `GillespieWorld` has almost no information about the coordinate of molecules. Thus, it simply ignores the given position, and just counts up the number of molecules here.\n", "\n", "`ParticleID` is a pair of `Integer`s named `lot` and `serial`." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(0, 1L)\n", "True\n" ] } ], "source": [ "print(pid6.lot(), pid6.serial())\n", "print(pid6 == ParticleID((0, 1)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Particle simulators, i.e. `spatiocyte`, `bd` and `egfrd`, provide an interface to access a particle by its id. `has_particle` returns if a particles exists or not for the given `ParticleID`." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True\n", "True\n", "True\n" ] } ], "source": [ "# print(w1.has_particle(pid1))\n", "# print(w2.has_particle(pid2))\n", "print(w3.has_particle(pid3)) # => True\n", "print(w4.has_particle(pid4)) # => True\n", "# print(w5.has_particle(pid5))\n", "print(w6.has_particle(pid6)) # => True" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After checking the existency, you can get the partcle by `get_particle` as follows." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# pid1, p1 = w1.get_particle(pid1)\n", "# pid2, p2 = w2.get_particle(pid2)\n", "pid3, p3 = w3.get_particle(pid3)\n", "pid4, p4 = w4.get_particle(pid4)\n", "# pid5, p5 = w5.get_particle(pid5)\n", "pid6, p6 = w6.get_particle(pid6)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "`Particle` consists of `species`, `position`, `radius` and `D`." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(u'A', (0.5062278801751902, 0.5080682368868706, 0.5), 0.0025, 1.0)\n", "(u'A', (0.5, 0.5, 0.5), 0.0025, 1.0)\n", "(u'A', (0.5, 0.5, 0.5), 0.0025, 1.0)\n" ] } ], "source": [ "# print(p1.species().serial(), tuple(p1.position()), p1.radius(), p1.D())\n", "# print(p2.species().serial(), tuple(p2.position()), p2.radius(), p2.D())\n", "print(p3.species().serial(), tuple(p3.position()), p3.radius(), p3.D())\n", "print(p4.species().serial(), tuple(p4.position()), p4.radius(), p4.D())\n", "# print(p5.species().serial(), tuple(p5.position()), p5.radius(), p5.D())\n", "print(p6.species().serial(), tuple(p6.position()), p6.radius(), p6.D())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the case of `spatiocyte`, a particle position is automatically round to the center of the voxel nearest to the given position.\n", "\n", "You can even move the position of the particle. `update_particle` replace the particle specified with the given `ParticleID` with the given `Particle` and return `False`. If no corresponding particle is found, create new particle and return `True`. If you give a `Particle` with the different type of `Species`, the `Species` of the `Particle` will be also changed." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "False\n", "False\n", "False\n" ] } ], "source": [ "newp = Particle(sp1, Real3(0.3, 0.3, 0.3), 0.0025, 1)\n", "# print(w1.update_particle(pid1, newp))\n", "# print(w2.update_particle(pid2, newp))\n", "print(w3.update_particle(pid3, newp))\n", "print(w4.update_particle(pid4, newp))\n", "# print(w5.update_particle(pid5, newp))\n", "print(w6.update_particle(pid6, newp))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`list_particles` and `list_particles_exact` return a list of pairs of `ParticleID` and `Particle` in the `World`. `World` automatically makes up for the gap with random numbers. For example, `GillespieWorld` returns a list of positions randomly distributed in the `World` size." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[(, )]\n", "[(, )]\n", "[(, )]\n", "[(, )]\n", "[(, )]\n" ] } ], "source": [ "print(w1.list_particles_exact(sp1))\n", "# print(w2.list_particles_exact(sp1)) # ODEWorld has no member named list_particles\n", "print(w3.list_particles_exact(sp1))\n", "print(w4.list_particles_exact(sp1))\n", "print(w5.list_particles_exact(sp1))\n", "print(w6.list_particles_exact(sp1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can remove a specific particle with `remove_particle`." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "False\n", "False\n", "False\n" ] } ], "source": [ "# w1.remove_particle(pid1)\n", "# w2.remove_particle(pid2)\n", "w3.remove_particle(pid3)\n", "w4.remove_particle(pid4)\n", "# w5.remove_particle(pid5)\n", "w6.remove_particle(pid6)\n", "# print(w1.has_particle(pid1))\n", "# print(w2.has_particle(pid2))\n", "print(w3.has_particle(pid3)) # => False\n", "print(w4.has_particle(pid4)) # => False\n", "# print(w5.has_particle(pid5))\n", "print(w6.has_particle(pid6)) # => False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.3. Lattice-based Coordinate\n", "\n", "In addition to the common interface, each `World` can have their own interfaces. As an example, we explain methods to handle lattice-based coordinate here. `SpatiocyteWorld` is based on a space discretized to hexiagonal close packing lattices, `LatticeSpace`." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "w = SpatiocyteWorld(Real3(1, 2, 3), voxel_radius=0.01)\n", "w.bind_to(m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The size of a single lattice, called `Voxel`, can be obtained by `voxel_radius()`. `SpatiocyteWorld` has methods to get the numbers of rows, columns, and layers. These sizes are automatically calculated based on the given `edge_lengths` at the construction." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.01\n", "(62, 152, 116)\n", "(62, 152, 116)\n", "1093184\n" ] } ], "source": [ "print(w.voxel_radius()) # => 0.01\n", "print(tuple(w.shape())) # => (62, 152, 116)\n", "print(w.col_size(), w.row_size(), w.layer_size()) # => (62, 152, 116)\n", "print(w.size()) # => 1093184 = 62 * 152 * 116" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A position in the lattice-based space is treated as an `Integer3`, column, row and layer, called a global coordinate. Thus, `SpatiocyteWorld` provides the function to convert the `Real3` into a lattice-based coordinate." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(31, 25, 29)\n", "(0.5062278801751902, 0.5080682368868706, 0.5)\n" ] } ], "source": [ "p1 = Real3(0.5, 0.5, 0.5)\n", "g1 = w.position2global(p1)\n", "p2 = w.global2position(g1)\n", "print(tuple(g1)) # => (31, 25, 29)\n", "print(tuple(p2)) # => (0.5062278801751902, 0.5080682368868706, 0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In `SpatiocyteWorld`, the global coordinate is translated to a single integer. It is just called a coordinate. You can also treat the coordinate as in the same way with a global coordinate." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "278033\n", "(0.5062278801751902, 0.5080682368868706, 0.5)\n", "(31, 25, 29)\n" ] } ], "source": [ "p1 = Real3(0.5, 0.5, 0.5)\n", "c1 = w.position2coordinate(p1)\n", "p2 = w.coordinate2position(c1)\n", "g1 = w.coord2global(c1)\n", "print(c1) # => 278033\n", "print(tuple(p2)) # => (0.5062278801751902, 0.5080682368868706, 0.5)\n", "print(tuple(g1)) # => (31, 25, 29)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With these coordinates, you can handle a `Voxel`, which represents a `Particle` object. Instead of `new_particle`, `new_voxel` provides the way to create a new `Voxel` with a coordinate." ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(, , True)\n" ] } ], "source": [ "c1 = w.position2coordinate(Real3(0.5, 0.5, 0.5))\n", "((pid, v), is_succeeded) = w.new_voxel(Species(\"A\"), c1)\n", "print(pid, v, is_succeeded)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A `Voxel` consists of `species`, `coordinate`, `radius` and `D`." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(u'A', 278033, 0.0025, 1.0)\n" ] } ], "source": [ "print(v.species().serial(), v.coordinate(), v.radius(), v.D()) # => (u'A', 278033, 0.0025, 1.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, you can get a voxel and list voxels with `get_voxel` and `list_voxels_exact` similar to `get_particle` and `list_particles_exact`." ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n", "[(, )]\n", "(, )\n" ] } ], "source": [ "print(w.num_voxels_exact(Species(\"A\")))\n", "print(w.list_voxels_exact(Species(\"A\")))\n", "print(w.get_voxel(pid))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can move and update the voxel with `update_voxel` corresponding to `update_particle`." ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "278058\n", "(u'A', 278058, 0.0025, 1.0)\n", "1\n" ] } ], "source": [ "c2 = w.position2coordinate(Real3(0.5, 0.5, 1.0))\n", "w.update_voxel(pid, Voxel(v.species(), c2, v.radius(), v.D()))\n", "pid, newv = w.get_voxel(pid)\n", "print(c2) # => 278058\n", "print(newv.species().serial(), newv.coordinate(), newv.radius(), newv.D()) # => (u'A', 278058, 0.0025, 1.0)\n", "print(w.num_voxels_exact(Species(\"A\"))) # => 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, `remove_voxel` remove a voxel as `remove_particle` does." ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True\n", "False\n" ] } ], "source": [ "print(w.has_voxel(pid)) # => True\n", "w.remove_voxel(pid)\n", "print(w.has_voxel(pid)) # => False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.4 Structure" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [], "source": [ "w1 = GillespieWorld()\n", "w2 = ODEWorld()\n", "w3 = SpatiocyteWorld()\n", "w4 = BDWorld()\n", "w5 = MesoscopicWorld()\n", "w6 = EGFRDWorld()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By using a `Shape` object, you can confine initial positions of molecules to a part of `World`. In the case below, 60 molecules are positioned inside the given `Sphere`. Diffusion of the molecules placed here is **NOT** restricted in the `Shape`. This `Shape` is only for the initialization." ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sp1 = Species(\"A\", \"0.0025\", \"1\")\n", "sphere = Sphere(Real3(0.5, 0.5, 0.5), 0.3)\n", "w1.add_molecules(sp1, 60, sphere)\n", "w2.add_molecules(sp1, 60, sphere)\n", "w3.add_molecules(sp1, 60, sphere)\n", "w4.add_molecules(sp1, 60, sphere)\n", "w5.add_molecules(sp1, 60, sphere)\n", "w6.add_molecules(sp1, 60, sphere)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A property of `Species`, `'location'`, is available to restrict diffusion of molecules. `'location'` is not fully supported yet, but only supported in `spatiocyte` and `meso`. `add_structure` defines a new structure given as a pair of `Species` and `Shape`." ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": true }, "outputs": [], "source": [ "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()`\n", "w3.add_structure(Species(\"M\"), membrane)\n", "w5.add_structure(Species(\"M\"), membrane)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After defining a structure, you can bind molecules to the structure as follows:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": true }, "outputs": [], "source": [ "sp2 = Species(\"B\", \"0.0025\", \"0.1\", \"M\") # `'location'` is the fourth argument\n", "w3.add_molecules(sp2, 60)\n", "w5.add_molecules(sp2, 60)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "The molecules bound to a `Species` named `B` diffuse on a structure named `M`, which has a shape of `SphericalSurface` (a hollow sphere). In `spatiocyte`, a structure is represented as a set of particles with `Species` `M` occupying a voxel. It means that molecules not belonging to the structure is not able to overlap the voxel and it causes a collision. On the other hand, in `meso`, a structure means a list of subvolumes. Thus, a structure doesn't avoid an incursion of other particles." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.5. Random Number Generator\n", "\n", "A random number generator is also a part of `World`. All `World` except `ODEWorld` store a random number generator, and updates it when the simulation needs a random value. On E-Cell4, only one class `GSLRandomNumberGenerator` is implemented as a random number generator." ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[6, 1, 2, 6, 2, 3, 6, 5, 4, 5, 5, 4, 2, 5, 4, 2, 3, 3, 2, 2]\n" ] } ], "source": [ "rng1 = GSLRandomNumberGenerator()\n", "print([rng1.uniform_int(1, 6) for _ in range(20)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With no argument, the random number generator is always initialized with a seed, `0`." ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[6, 1, 2, 6, 2, 3, 6, 5, 4, 5, 5, 4, 2, 5, 4, 2, 3, 3, 2, 2]\n" ] } ], "source": [ "rng2 = GSLRandomNumberGenerator()\n", "print([rng2.uniform_int(1, 6) for _ in range(20)]) # => same as above" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can initialize the seed with an integer as follows:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[6, 5, 2, 4, 1, 1, 3, 5, 2, 6, 4, 1, 2, 5, 2, 5, 1, 2, 2, 6]\n" ] } ], "source": [ "rng2 = GSLRandomNumberGenerator()\n", "rng2.seed(15)\n", "print([rng2.uniform_int(1, 6) for _ in range(20)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When you call the `seed` function with no input, the seed is drawn from the current time." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[6, 2, 1, 5, 6, 2, 6, 2, 6, 6, 5, 6, 2, 3, 5, 3, 6, 1, 4, 3]\n" ] } ], "source": [ "rng2 = GSLRandomNumberGenerator()\n", "rng2.seed()\n", "print([rng2.uniform_int(1, 6) for _ in range(20)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`GSLRandomNumberGenerator` provides several ways to get a random number." ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0303352042101\n", "33\n", "0.893555545521\n" ] } ], "source": [ "print(rng1.uniform(0.0, 1.0))\n", "print(rng1.uniform_int(0, 100))\n", "print(rng1.gaussian(1.0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`World` accepts a random number generator at the construction. As a default, `GSLRandomNumberGenerator()` is used. Thus, when you don't give a generator, behavior of the simulation is always same (determinisitc)." ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": true }, "outputs": [], "source": [ "rng = GSLRandomNumberGenerator()\n", "rng.seed()\n", "w1 = GillespieWorld(Real3(1, 1, 1), rng=rng)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can access the `GSLRandomNumberGenerator` in a `World` through `rng` function." ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.966634188779\n" ] } ], "source": [ "print(w1.rng().uniform(0.0, 1.0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`rng()` returns a shared pointer to the `GSLRandomNumberGenerator`. Thus, in the example above, `rng` and `w1.rng()` point exactly the same thing." ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }