diff --git a/ecell4/core/VoxelSpaceBase.cpp b/ecell4/core/VoxelSpaceBase.cpp
index 941d4b55a..a799e3b64 100644
--- a/ecell4/core/VoxelSpaceBase.cpp
+++ b/ecell4/core/VoxelSpaceBase.cpp
@@ -21,6 +21,14 @@ Integer VoxelSpaceBase::num_molecules(const Species &sp) const
     Integer count(0);
     SpeciesExpressionMatcher sexp(sp);
 
+    {
+        const auto cnt = sexp.count(vacant_->species());
+        if (cnt > 0)
+        {
+            count += vacant_->size() * cnt;
+        }
+    }
+
     for (voxel_pool_map_type::const_iterator itr(voxel_pools_.begin());
          itr != voxel_pools_.end(); ++itr)
     {
@@ -62,6 +70,11 @@ bool VoxelSpaceBase::has_voxel(const ParticleID &pid) const
 
 Integer VoxelSpaceBase::num_voxels_exact(const Species &sp) const
 {
+    if (sp == vacant_->species())
+    {
+        return vacant_->size();
+    }
+
     {
         voxel_pool_map_type::const_iterator itr(voxel_pools_.find(sp));
         if (itr != voxel_pools_.end())
@@ -87,6 +100,11 @@ Integer VoxelSpaceBase::num_voxels(const Species &sp) const
     Integer count(0);
     SpeciesExpressionMatcher sexp(sp);
 
+    if (sexp.match(vacant_->species()))
+    {
+        count += vacant_->size();
+    }
+
     for (voxel_pool_map_type::const_iterator itr(voxel_pools_.begin());
          itr != voxel_pools_.end(); ++itr)
     {
@@ -112,6 +130,11 @@ Integer VoxelSpaceBase::num_voxels() const
 {
     Integer count(0);
 
+    if (vacant_->species().serial() != "")
+    {
+        count += vacant_->size();
+    }
+
     for (voxel_pool_map_type::const_iterator itr(voxel_pools_.begin());
          itr != voxel_pools_.end(); ++itr)
     {
@@ -179,7 +202,7 @@ VoxelSpaceBase::list_voxels_exact(const Species &sp) const
 
 boost::shared_ptr<VoxelPool> VoxelSpaceBase::find_voxel_pool(const Species &sp)
 {
-    if (sp.serial() == "")
+    if (sp == vacant_->species())
         return vacant_;
 
     voxel_pool_map_type::iterator itr(voxel_pools_.find(sp));
@@ -193,7 +216,7 @@ boost::shared_ptr<VoxelPool> VoxelSpaceBase::find_voxel_pool(const Species &sp)
 boost::shared_ptr<const VoxelPool>
 VoxelSpaceBase::find_voxel_pool(const Species &sp) const
 {
-    if (sp.serial() == "")
+    if (sp == vacant_->species())
         return vacant_;
 
     voxel_pool_map_type::const_iterator itr(voxel_pools_.find(sp));
diff --git a/ecell4/core/tests/LatticeSpace_test.cpp b/ecell4/core/tests/LatticeSpace_test.cpp
index 0b60b11d2..c76855f7c 100644
--- a/ecell4/core/tests/LatticeSpace_test.cpp
+++ b/ecell4/core/tests/LatticeSpace_test.cpp
@@ -40,6 +40,8 @@ BOOST_AUTO_TEST_CASE(LatticeSpace_test_constructor) { ; }
 BOOST_AUTO_TEST_CASE(CheckVacantSize)
 {
     BOOST_CHECK_EQUAL(space.actual_size(), space.vacant()->size());
+    BOOST_CHECK_EQUAL(space.num_voxels_exact(Species("")),
+                      space.vacant()->size());
 }
 
 BOOST_AUTO_TEST_CASE(GetVoxel)
diff --git a/ecell4/core/tests/OffLatticeSpace_test.cpp b/ecell4/core/tests/OffLatticeSpace_test.cpp
index 2746355fb..d580ce7f1 100644
--- a/ecell4/core/tests/OffLatticeSpace_test.cpp
+++ b/ecell4/core/tests/OffLatticeSpace_test.cpp
@@ -41,7 +41,7 @@ struct Fixture
         for (int i(1); i < 10; ++i)
             adjoining_pairs.push_back(std::make_pair(i - 1, i));
         space = OffLatticeSpace(voxel_radius, base, positions, adjoining_pairs);
-        space.make_molecular_type(species, "");
+        space.make_molecular_type(species, "Base");
     }
 };
 
@@ -70,7 +70,7 @@ BOOST_AUTO_TEST_CASE(OffLatticeSpace_test_voxelspacebase)
     BOOST_CHECK_EQUAL(space.list_species().size(), 1);
     BOOST_CHECK_EQUAL(space.num_voxels_exact(species), 1);
     BOOST_CHECK_EQUAL(space.num_voxels(species), 1);
-    BOOST_CHECK_EQUAL(space.num_voxels(), 1);
+    BOOST_CHECK_EQUAL(space.num_voxels(), 10);
 
     BOOST_CHECK(space.has_voxel(pid));
     BOOST_CHECK(!space.has_voxel(sidgen()));
diff --git a/ecell4/python_api/spatiocyte.cpp b/ecell4/python_api/spatiocyte.cpp
index df0ff028f..4adc77b50 100644
--- a/ecell4/python_api/spatiocyte.cpp
+++ b/ecell4/python_api/spatiocyte.cpp
@@ -1,12 +1,12 @@
 #include "python_api.hpp"
 
 #include <ecell4/core/OffLatticeSpace.hpp>
-#include <ecell4/spatiocyte/SpatiocyteReactions.hpp>
+#include <ecell4/spatiocyte/OffLattice.hpp>
 #include <ecell4/spatiocyte/SpatiocyteFactory.hpp>
+#include <ecell4/spatiocyte/SpatiocyteReactions.hpp>
 #include <ecell4/spatiocyte/SpatiocyteSimulator.hpp>
 #include <ecell4/spatiocyte/SpatiocyteWorld.hpp>
 #include <ecell4/spatiocyte/Voxel.hpp>
-#include <ecell4/spatiocyte/OffLattice.hpp>
 
 #include "simulator.hpp"
 #include "simulator_factory.hpp"
@@ -21,32 +21,28 @@ namespace ecell4
 namespace python_api
 {
 
-static inline
-void define_reaction_info(py::module& m)
+static inline void define_reaction_info(py::module &m)
 {
     using container_type = ReactionInfo::container_type;
     py::class_<ReactionInfo>(m, "ReactionInfo")
-        .def(py::init<const Real, const container_type&, const container_type&>(),
-                py::arg("t"), py::arg("reactants"), py::arg("products"))
+        .def(py::init<const Real, const container_type &,
+                      const container_type &>(),
+             py::arg("t"), py::arg("reactants"), py::arg("products"))
         .def("t", &ReactionInfo::t)
         .def("reactants", &ReactionInfo::reactants)
         .def("products", &ReactionInfo::products)
         .def(py::pickle(
-            [](const ReactionInfo& self)
-            {
-                return py::make_tuple(self.t(), self.reactants(), self.products());
+            [](const ReactionInfo &self) {
+                return py::make_tuple(self.t(), self.reactants(),
+                                      self.products());
             },
-            [](py::tuple t)
-            {
+            [](py::tuple t) {
                 if (t.size() != 3)
                     throw std::runtime_error("Invalid state");
-                return ReactionInfo(
-                    t[0].cast<Real>(),
-                    t[1].cast<container_type>(),
-                    t[2].cast<container_type>()
-                );
-            }
-        ));
+                return ReactionInfo(t[0].cast<Real>(),
+                                    t[1].cast<container_type>(),
+                                    t[2].cast<container_type>());
+            }));
 
     py::class_<ReactionInfo::Item>(m, "ReactionInfoItem")
         .def_readonly("pid", &ReactionInfo::Item::pid)
@@ -54,28 +50,28 @@ void define_reaction_info(py::module& m)
         .def_readonly("voxel", &ReactionInfo::Item::voxel);
 }
 
-static inline
-void define_spatiocyte_factory(py::module& m)
+static inline void define_spatiocyte_factory(py::module &m)
 {
     py::class_<SpatiocyteFactory> factory(m, "SpatiocyteFactory");
     factory
         .def(py::init<const Real>(),
-                py::arg("voxel_radius") = SpatiocyteFactory::default_voxel_radius())
+             py::arg("voxel_radius") =
+                 SpatiocyteFactory::default_voxel_radius())
         .def("rng", &SpatiocyteFactory::rng);
     define_factory_functions(factory);
 
     m.attr("Factory") = factory;
 }
 
-static inline
-void define_spatiocyte_simulator(py::module& m)
+static inline void define_spatiocyte_simulator(py::module &m)
 {
     py::class_<SpatiocyteSimulator, Simulator, PySimulator<SpatiocyteSimulator>,
-        boost::shared_ptr<SpatiocyteSimulator>> simulator(m, "SpatiocyteSimulator");
-    simulator
-        .def(py::init<boost::shared_ptr<SpatiocyteWorld>>(), py::arg("w"))
-        .def(py::init<boost::shared_ptr<SpatiocyteWorld>, boost::shared_ptr<Model>>(),
-                py::arg("w"), py::arg("m"))
+               boost::shared_ptr<SpatiocyteSimulator>>
+        simulator(m, "SpatiocyteSimulator");
+    simulator.def(py::init<boost::shared_ptr<SpatiocyteWorld>>(), py::arg("w"))
+        .def(py::init<boost::shared_ptr<SpatiocyteWorld>,
+                      boost::shared_ptr<Model>>(),
+             py::arg("w"), py::arg("m"))
         .def("last_reactions", &SpatiocyteSimulator::last_reactions)
         .def("set_t", &SpatiocyteSimulator::set_t);
     define_simulator_functions(simulator);
@@ -83,34 +79,40 @@ void define_spatiocyte_simulator(py::module& m)
     m.attr("Simulator") = simulator;
 }
 
-static inline
-void define_spatiocyte_world(py::module& m)
+static inline void define_spatiocyte_world(py::module &m)
 {
-    py::class_<SpatiocyteWorld, ecell4::WorldInterface, PyWorldImpl<SpatiocyteWorld>,
-        boost::shared_ptr<SpatiocyteWorld>> world(m, "SpatiocyteWorld");
+    py::class_<SpatiocyteWorld, ecell4::WorldInterface,
+               PyWorldImpl<SpatiocyteWorld>, boost::shared_ptr<SpatiocyteWorld>>
+        world(m, "SpatiocyteWorld");
     world
-        .def(py::init<const Real3&>(),
-                py::arg("edge_lengths") = Real3(1.0, 1.0, 1.0))
-        .def(py::init<const Real3&, const Real&>(),
-                py::arg("edge_lengths"), py::arg("voxel_radius"))
-        .def(py::init<const Real3&, const Real&, const boost::shared_ptr<RandomNumberGenerator>&>(),
-                py::arg("edge_lengths"),
-                py::arg("voxel_radius"),
-                py::arg("rng"))
+        .def(py::init<const Real3 &>(),
+             py::arg("edge_lengths") = Real3(1.0, 1.0, 1.0))
+        .def(py::init<const Real3 &, const Real &>(), py::arg("edge_lengths"),
+             py::arg("voxel_radius"))
+        .def(py::init<const Real3 &, const Real &,
+                      const boost::shared_ptr<RandomNumberGenerator> &>(),
+             py::arg("edge_lengths"), py::arg("voxel_radius"), py::arg("rng"))
         .def(py::init<const std::string>(), py::arg("filename"))
 
-        .def("new_particle", (boost::optional<ParticleID> (SpatiocyteWorld::*)(const Particle&)) &SpatiocyteWorld::new_particle)
-        .def("new_particle", (boost::optional<ParticleID> (SpatiocyteWorld::*)(const Species&, const Real3&)) &SpatiocyteWorld::new_particle)
+        .def("new_particle", (boost::optional<ParticleID>(SpatiocyteWorld::*)(
+                                 const Particle &)) &
+                                 SpatiocyteWorld::new_particle)
+        .def("new_particle", (boost::optional<ParticleID>(SpatiocyteWorld::*)(
+                                 const Species &, const Real3 &)) &
+                                 SpatiocyteWorld::new_particle)
         .def("remove_particle", &SpatiocyteWorld::remove_particle)
-        .def("list_structure_particles", &SpatiocyteWorld::list_structure_particles)
-        .def("list_non_structure_particles", &SpatiocyteWorld::list_non_structure_particles)
+        .def("list_structure_particles",
+             &SpatiocyteWorld::list_structure_particles)
+        .def("list_non_structure_particles",
+             &SpatiocyteWorld::list_non_structure_particles)
         .def("update_particle", &SpatiocyteWorld::update_particle)
         .def("add_molecules",
-                (bool (SpatiocyteWorld::*)(const Species&, const Integer&))
-                &SpatiocyteWorld::add_molecules)
+             (bool (SpatiocyteWorld::*)(const Species &, const Integer &)) &
+                 SpatiocyteWorld::add_molecules)
         .def("add_molecules",
-                (bool (SpatiocyteWorld::*)(const Species&, const Integer&, const boost::shared_ptr<const Shape>))
-                &SpatiocyteWorld::add_molecules)
+             (bool (SpatiocyteWorld::*)(const Species &, const Integer &,
+                                        const boost::shared_ptr<const Shape>)) &
+                 SpatiocyteWorld::add_molecules)
         .def("remove_molecules", &SpatiocyteWorld::remove_molecules)
         .def("voxel_volume", &SpatiocyteWorld::voxel_volume)
         .def("get_volume", &SpatiocyteWorld::get_volume)
@@ -118,14 +120,17 @@ void define_spatiocyte_world(py::module& m)
         .def("get_species_at", &SpatiocyteWorld::get_species_at)
         .def("has_particle_at", &SpatiocyteWorld::has_particle_at)
         .def("set_value", &SpatiocyteWorld::set_value)
-        .def("new_particle", (boost::optional<ParticleID> (SpatiocyteWorld::*)(const Species&, const Voxel&))&SpatiocyteWorld::new_particle)
+        .def("new_particle", (boost::optional<ParticleID>(SpatiocyteWorld::*)(
+                                 const Species &, const Voxel &)) &
+                                 SpatiocyteWorld::new_particle)
         .def("new_voxel",
-            (boost::optional<ParticleID> (SpatiocyteWorld::*)(const Species&, const Voxel&))&SpatiocyteWorld::new_particle,
-            R"pbdoc(
+             (boost::optional<ParticleID>(SpatiocyteWorld::*)(const Species &,
+                                                              const Voxel &)) &
+                 SpatiocyteWorld::new_particle,
+             R"pbdoc(
                 .. deprecated::3.0
                    Use :func:`new_particle` instead.
-            )pbdoc"
-        )
+            )pbdoc")
         .def("new_voxel_structure", &SpatiocyteWorld::new_voxel_structure)
         .def("voxel_radius", &SpatiocyteWorld::voxel_radius)
         .def("size", &SpatiocyteWorld::size)
@@ -148,54 +153,60 @@ void define_spatiocyte_world(py::module& m)
         )pbdoc")
         .def("rng", &SpatiocyteWorld::rng)
         .def("add_offlattice",
-            [](SpatiocyteWorld& self, const Species& species, const OffLattice& offlattice)
-            {
-                self.add_space(offlattice.generate_space(species));
-            })
-        .def_static("calculate_voxel_volume", &SpatiocyteWorld::calculate_voxel_volume)
-        .def_static("calculate_hcp_lengths", &SpatiocyteWorld::calculate_hcp_lengths)
+             [](SpatiocyteWorld &self, const Species &species,
+                const OffLattice &offlattice) {
+                 const auto info = self.get_molecule_info(species);
+                 const auto updated = Species(species.serial(), info.radius,
+                                              info.D, info.loc, info.dimension);
+                 self.add_space(offlattice.generate_space(updated));
+             })
+        .def_static("calculate_voxel_volume",
+                    &SpatiocyteWorld::calculate_voxel_volume)
+        .def_static("calculate_hcp_lengths",
+                    &SpatiocyteWorld::calculate_hcp_lengths)
         .def_static("calculate_shape", &SpatiocyteWorld::calculate_shape)
-        .def_static("calculate_volume", &SpatiocyteWorld::calculate_volume);
-
-    m.def("create_spatiocyte_world_cell_list_impl", &create_spatiocyte_world_cell_list_impl);
-    m.def("create_spatiocyte_world_vector_impl", &create_spatiocyte_world_vector_impl);
-    m.def("create_spatiocyte_world_square_offlattice_impl", &allocate_spatiocyte_world_square_offlattice_impl);
+        .def_static("calculate_volume", &SpatiocyteWorld::calculate_volume)
+        .def("list_neighbors",
+             [](const SpatiocyteWorld &self, const Voxel &voxel) {
+                 std::vector<Voxel> list;
+                 for (auto i = 0; i < self.num_neighbors(voxel); ++i)
+                 {
+                     list.push_back(self.get_neighbor(voxel, i));
+                 }
+                 return list;
+             });
+
+    m.def("create_spatiocyte_world_cell_list_impl",
+          &create_spatiocyte_world_cell_list_impl);
+    m.def("create_spatiocyte_world_vector_impl",
+          &create_spatiocyte_world_vector_impl);
+    m.def("create_spatiocyte_world_square_offlattice_impl",
+          &allocate_spatiocyte_world_square_offlattice_impl);
 
     m.attr("World") = world;
 }
 
-static inline
-void define_voxel(py::module& m)
+static inline void define_voxel(py::module &m)
 {
     py::class_<Voxel>(m, "Voxel")
-        .def("position", &Voxel::position)
 #ifndef NDEBUG
         .def_readonly("coordinate", &Voxel::coordinate)
 #endif
-        .def("list_neighbors",
-            [](const Voxel& self)
-            {
-                std::vector<Voxel> list;
-                for (auto i = 0; i < self.num_neighbors(); ++i)
-                {
-                    list.push_back(self.get_neighbor(i));
-                }
-                return list;
-            });
+        .def("position", &Voxel::position);
 }
 
-static inline
-void define_offlattice(py::module& m)
+static inline void define_offlattice(py::module &m)
 {
     py::class_<OffLattice>(m, "OffLattice")
-        .def(py::init<const Real, const OffLattice::positions_type, const OffLattice::adjoining_pairs_type>())
+        .def(py::init<const Real, const OffLattice::positions_type,
+                      const OffLattice::adjoining_pairs_type>())
         .def(py::init<const Real, const OffLattice::positions_type>())
         .def("voxel_radius", &OffLattice::voxel_radius)
         .def("positions", &OffLattice::positions)
         .def("adjoining_pairs", &OffLattice::adjoining_pairs);
 }
 
-void setup_spatiocyte_module(py::module& m)
+void setup_spatiocyte_module(py::module &m)
 {
     define_reaction_info(m);
     define_offlattice(m);
@@ -205,7 +216,6 @@ void setup_spatiocyte_module(py::module& m)
     define_voxel(m);
 }
 
-}
-
-}
+} // namespace python_api
 
+} // namespace ecell4
diff --git a/ecell4/spatiocyte/ReactionEvent.cpp b/ecell4/spatiocyte/ReactionEvent.cpp
index e89119934..4836fcb5c 100644
--- a/ecell4/spatiocyte/ReactionEvent.cpp
+++ b/ecell4/spatiocyte/ReactionEvent.cpp
@@ -24,9 +24,12 @@ void ZerothOrderReactionEvent::fire_()
     {
         const MoleculeInfo info(world_->get_molecule_info(sp));
 
-        if (boost::shared_ptr<VoxelPool> location =
-                world_->find_voxel_pool(Species(info.loc)))
+        if (const auto space_and_location =
+                world_->find_space_and_voxel_pool(Species(info.loc)))
         {
+            const auto space = space_and_location->first;
+            const auto location = space_and_location->second;
+
             if (location->size() == 0)
             {
                 time_ += draw_dt();
@@ -35,8 +38,8 @@ void ZerothOrderReactionEvent::fire_()
 
             while (true)
             {
-                const Voxel voxel(world_->coordinate2voxel(
-                    world_->rng()->uniform_int(0, world_->size() - 1)));
+                const Voxel voxel(
+                    space, world_->rng()->uniform_int(0, space->size() - 1));
 
                 if (voxel.get_voxel_pool() != location)
                 {
diff --git a/ecell4/spatiocyte/SpatiocyteEvent.hpp b/ecell4/spatiocyte/SpatiocyteEvent.hpp
index c51b402f4..5ecd0e0ee 100644
--- a/ecell4/spatiocyte/SpatiocyteEvent.hpp
+++ b/ecell4/spatiocyte/SpatiocyteEvent.hpp
@@ -66,6 +66,7 @@ struct StepEvent : SpatiocyteEvent
 protected:
     boost::shared_ptr<Model> model_;
     boost::shared_ptr<SpatiocyteWorld> world_;
+    boost::weak_ptr<VoxelSpaceBase> space_;
     boost::shared_ptr<MoleculePool> mpool_;
 
     const Real alpha_;
@@ -87,9 +88,6 @@ struct StepEvent2D : StepEvent
                 const Species &species, const Real &t, const Real alpha = 1.0);
 
     void walk(const Real &alpha);
-
-protected:
-    std::vector<unsigned int> nids_; // neighbor indexes
 };
 
 struct ZerothOrderReactionEvent : SpatiocyteEvent
@@ -123,15 +121,20 @@ struct FirstOrderReactionEvent : SpatiocyteEvent
     ReactionInfo::Item choice()
     {
         const Species &species(rule_.reactants().at(0));
-        boost::shared_ptr<const MoleculePool> mt(
-            world_->find_molecule_pool(species));
-
-        const Integer i(rng_.lock()->uniform_int(0, mt->size() - 1));
-        const SpatiocyteWorld::coordinate_id_pair_type &info(mt->at(i));
-
-        // TODO: Calling coordinate2voxel() is invalid
-        return ReactionInfo::Item(info.pid, species,
-                                  world_->coordinate2voxel(info.coordinate));
+        if (const auto space_and_molecule_pool =
+                world_->find_space_and_molecule_pool(species))
+        {
+            const auto space = space_and_molecule_pool->first;
+            const auto molecule_pool = space_and_molecule_pool->second;
+
+            const auto i =
+                rng_.lock()->uniform_int(0, molecule_pool->size() - 1);
+            const auto &info = molecule_pool->at(i);
+
+            return ReactionInfo::Item(info.pid, species,
+                                      Voxel(space, info.coordinate));
+        }
+        throw "MoleculePool is not found";
     }
 
     boost::shared_ptr<SpatiocyteWorld> world_;
diff --git a/ecell4/spatiocyte/SpatiocyteReactions.cpp b/ecell4/spatiocyte/SpatiocyteReactions.cpp
index 472c55388..9395aed9d 100644
--- a/ecell4/spatiocyte/SpatiocyteReactions.cpp
+++ b/ecell4/spatiocyte/SpatiocyteReactions.cpp
@@ -12,8 +12,7 @@ namespace spatiocyte
 inline const std::string get_serial(boost::shared_ptr<SpatiocyteWorld> world,
                                     const Voxel &voxel)
 {
-    boost::shared_ptr<const VoxelPool> mtype(voxel.get_voxel_pool());
-    return mtype->is_vacant() ? "" : mtype->species().serial();
+    return voxel.get_voxel_pool()->species().serial();
 }
 
 inline const std::string get_location(boost::shared_ptr<SpatiocyteWorld> world,
@@ -22,8 +21,7 @@ inline const std::string get_location(boost::shared_ptr<SpatiocyteWorld> world,
     boost::shared_ptr<const VoxelPool> mtype(voxel.get_voxel_pool());
     if (mtype->is_vacant())
         return "";
-    boost::shared_ptr<const VoxelPool> ltype(mtype->location());
-    return ltype->is_vacant() ? "" : ltype->species().serial();
+    return mtype->location()->species().serial();
 }
 
 static inline void make_product(boost::shared_ptr<SpatiocyteWorld> world,
@@ -444,7 +442,7 @@ apply_second_order_reaction(boost::shared_ptr<SpatiocyteWorld> world,
         return apply_vanishment(world, reactant_item0, reactant_item1);
     case 1:
         return apply_ab2c(world, reactant_item0, reactant_item1,
-                          *(products.begin()));
+                          products.at(0));
     case 2:
         return apply_ab2cd(world, reactant_item0, reactant_item1,
                            products.at(0), products.at(1));
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index a29b7305a..af30ceccd 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -9,33 +9,40 @@ namespace ecell4
 namespace spatiocyte
 {
 
-void SpatiocyteWorld::add_space(std::unique_ptr<VoxelSpaceBase> space)
+void SpatiocyteWorld::add_space(std::unique_ptr<VoxelSpaceBase> uniq_space)
 {
-    for (coordinate_type i(0); i < space->size(); ++i)
+    boost::shared_ptr<VoxelSpaceBase> space(uniq_space.release());
+
+    for (auto i(0); i < space->size(); ++i)
     {
-        const Real3 position(space->coordinate2position(i));
-        const coordinate_type nearest(
-            get_root()->position2coordinate(position));
+        const Voxel voxel(space, i);
+        const auto position(voxel.position());
+        const Voxel nearest(get_root(),
+                            get_root()->position2coordinate(position));
 
-        for (Integer j(0); j < get_root()->num_neighbors(nearest); ++j)
+        for (Integer j(0); j < num_neighbors(nearest); ++j)
         {
-            const coordinate_type neighbor(
-                get_root()->get_neighbor(nearest, j));
-            if (length(get_root()->coordinate2position(neighbor) - position) <
-                voxel_radius() * 2)
-                interfaces_.add(neighbor, i + size_);
+            const auto neighbor(get_neighbor(nearest, j));
+            if (length(neighbor.position() - position) < voxel_radius() * 2)
+            {
+                interfaces_.add(neighbor, voxel);
+            }
         }
     }
 
     for (const auto &interface : interfaces_)
     {
-        std::vector<coordinate_type> neighbors;
-        for (Integer i(0); i < get_root()->num_neighbors(interface.first); ++i)
+        const Voxel voxel(interface.first);
+
+        std::vector<Voxel> neighbors;
+        for (auto i(0); i < num_neighbors(voxel); ++i)
         {
-            const coordinate_type neighbor(
-                get_root()->get_neighbor(interface.first, i));
+            const auto neighbor(get_neighbor(voxel, i));
+
             if (!interfaces_.find(neighbor))
+            {
                 neighbors.push_back(neighbor);
+            }
         }
 
         for (const auto &adjoining : interface.second)
@@ -45,7 +52,7 @@ void SpatiocyteWorld::add_space(std::unique_ptr<VoxelSpaceBase> space)
     }
 
     size_ += space->size();
-    spaces_.push_back(space_type(space.release()));
+    spaces_.push_back(space);
 }
 
 void SpatiocyteWorld::set_value(const Species &sp, const Real value)
@@ -149,21 +156,29 @@ bool SpatiocyteWorld::add_molecules(const Species &sp, const Integer &num)
 
     const MoleculeInfo info(get_molecule_info(sp));
 
-    Integer count(0);
-    while (count < num)
+    if (const auto space_and_location =
+            find_space_and_voxel_pool(Species(info.loc)))
     {
-        const Voxel voxel(coordinate2voxel(rng()->uniform_int(0, size() - 1)));
+        const auto space = space_and_location->first;
+        const auto location = space_and_location->second;
 
-        if (voxel.get_voxel_pool()->species().serial() != info.loc)
-        {
-            continue;
-        }
-        else if (new_particle(sp, voxel))
+        if (location->size() < num)
+            return false;
+
+        auto count(0);
+        while (count < num)
         {
-            ++count;
+            const Voxel voxel(space, rng()->uniform_int(0, space->size() - 1));
+
+            if (voxel.get_voxel_pool() != location)
+                continue;
+
+            if (new_particle(sp, voxel))
+                ++count;
         }
+        return true;
     }
-    return true;
+    return false;
 }
 
 bool SpatiocyteWorld::add_molecules(const Species &sp, const Integer &num,
@@ -227,26 +242,29 @@ SpatiocyteWorld::add_structure3(const Species &sp, const std::string &location,
                                 const boost::shared_ptr<const Shape> shape)
 {
     Integer count(0);
-    for (coordinate_type coord(0); coord < size(); ++coord)
+    for (const auto &space : spaces_)
     {
-        const Voxel voxel(coordinate2voxel(coord));
-
-        if (!this->is_inside(coord) || shape->is_inside(voxel.position()) > 0)
-        {
-            continue;
-        }
-
-        if (voxel.get_voxel_pool()->species().serial() != location)
+        for (auto coord(0); coord < space->size(); ++coord)
         {
-            throw NotSupported("Mismatch in the location. Failed to place '" +
-                               sp.serial() + "' to '" +
-                               voxel.get_voxel_pool()->species().serial() +
-                               "'. " + "'" + location + "' is expected.");
-            continue;
+            const Voxel voxel(space, coord);
+            // should check if coord doesn't point at the boundary.
+            if (!space->is_inside(coord) ||
+                shape->is_inside(voxel.position()) > 0)
+                continue;
+
+            if (voxel.get_voxel_pool()->species().serial() != location)
+            {
+                throw NotSupported(
+                    "Mismatch in the location. Failed to place '" +
+                    sp.serial() + "' to '" +
+                    voxel.get_voxel_pool()->species().serial() + "'. " + "'" +
+                    location + "' is expected.");
+                continue;
+            }
+
+            if (new_voxel_structure(sp, voxel))
+                ++count;
         }
-
-        if (new_voxel_structure(sp, voxel))
-            ++count;
     }
     return count;
 }
@@ -256,25 +274,27 @@ SpatiocyteWorld::add_structure2(const Species &sp, const std::string &location,
                                 const boost::shared_ptr<const Shape> shape)
 {
     Integer count(0);
-    for (coordinate_type coord(0); coord < size(); ++coord)
+    for (const auto &space : spaces_)
     {
-        const Voxel voxel(coordinate2voxel(coord));
-        if (!this->is_inside(coord) || !is_surface_voxel(voxel, shape))
-        {
-            continue;
-        }
-
-        if (voxel.get_voxel_pool()->species().serial() != location)
+        for (auto coord(0); coord < space->size(); ++coord)
         {
-            throw NotSupported("Mismatch in the location. Failed to place '" +
-                               sp.serial() + "' to '" +
-                               voxel.get_voxel_pool()->species().serial() +
-                               "'. " + "'" + location + "' is expected.");
-            continue;
+            const Voxel voxel(space, coord);
+            // should check if coord doesn't point at the boundary.
+            if (!space->is_inside(coord) || !is_surface_voxel(voxel, shape))
+                continue;
+            if (voxel.get_voxel_pool()->species().serial() != location)
+            {
+                throw NotSupported(
+                    "Mismatch in the location. Failed to place '" +
+                    sp.serial() + "' to '" +
+                    voxel.get_voxel_pool()->species().serial() + "'. " + "'" +
+                    location + "' is expected.");
+                continue;
+            }
+
+            if (new_voxel_structure(sp, voxel))
+                ++count;
         }
-
-        if (new_voxel_structure(sp, voxel))
-            ++count;
     }
     return count;
 }
@@ -286,8 +306,8 @@ bool SpatiocyteWorld::is_surface_voxel(
     if (L > 0 || L < -2 * voxel_radius())
         return false;
 
-    for (Integer i(0); i < voxel.num_neighbors(); ++i)
-        if (shape->is_inside(voxel.get_neighbor(i).position()) > 0)
+    for (Integer i(0); i < num_neighbors(voxel); ++i)
+        if (shape->is_inside(get_neighbor(voxel, i).position()) > 0)
             return true;
 
     return false;
@@ -301,20 +321,24 @@ void SpatiocyteWorld::remove_molecules(const Species &sp, const Integer &num)
             "The number of molecules must be positive.");
     }
 
-    boost::shared_ptr<const MoleculePool> mtype(find_molecule_pool(sp));
-    if (mtype->size() < num)
+    if (const auto space_and_molecule_pool = find_space_and_molecule_pool(sp))
     {
-        throw std::invalid_argument(
-            "The number of molecules cannot be negative.");
-    }
+        const auto space = space_and_molecule_pool->first;
+        const auto molecule_pool = space_and_molecule_pool->second;
 
-    Integer count(0);
-    while (count < num)
-    {
-        const Integer idx(rng_->uniform_int(0, mtype->size() - 1));
-        if (coordinate2voxel(mtype->at(idx).coordinate).clear())
+        if (molecule_pool->size() < num)
+            throw std::invalid_argument(
+                "The number of molecules cannot be negative.");
+
+        auto count(0);
+        while (count < num)
         {
-            ++count;
+            const auto idx(rng_->uniform_int(0, molecule_pool->size() - 1));
+            const Voxel voxel(space, molecule_pool->at(idx).coordinate);
+            if (voxel.clear())
+            {
+                ++count;
+            }
         }
     }
 }
@@ -322,22 +346,32 @@ void SpatiocyteWorld::remove_molecules(const Species &sp, const Integer &num)
 boost::optional<Voxel> SpatiocyteWorld::check_neighbor(const Voxel &voxel,
                                                        const std::string &loc)
 {
-    const std::size_t num_neighbors(voxel.num_neighbors());
+    const std::size_t num(num_neighbors(voxel));
 
     std::vector<Voxel> tmp;
-    tmp.reserve(num_neighbors);
+    tmp.reserve(num);
 
-    for (unsigned int rnd(0); rnd < num_neighbors; ++rnd)
+    for (unsigned int rnd(0); rnd < num; ++rnd)
     {
-        const Voxel neighbor(voxel.get_neighbor(rnd));
+        const Voxel neighbor(get_neighbor(voxel, rnd));
         boost::shared_ptr<const VoxelPool> mt(neighbor.get_voxel_pool());
-        const std::string serial(mt->is_vacant() ? "" : mt->species().serial());
-        if (serial == loc)
+        if (mt->species().serial() == loc)
         {
             tmp.push_back(neighbor);
         }
     }
 
+    if (const auto neighbors = neighbors_.find(voxel))
+    {
+        for (const auto &neighbor : *neighbors)
+        {
+            if (neighbor.get_voxel_pool()->species().serial() == loc)
+            {
+                tmp.push_back(neighbor);
+            }
+        }
+    }
+
     if (tmp.size() == 0)
     {
         return boost::none;
@@ -346,6 +380,33 @@ boost::optional<Voxel> SpatiocyteWorld::check_neighbor(const Voxel &voxel,
     return tmp[rng()->uniform_int(0, tmp.size() - 1)];
 }
 
+const Voxel
+SpatiocyteWorld::get_neighbor_randomly(const Voxel &voxel,
+                                       Shape::dimension_kind dimension)
+{
+    std::vector<Voxel> neighbors;
+    for (Integer idx = 0; idx < num_neighbors(voxel); ++idx)
+    {
+        const Voxel neighbor = get_neighbor(voxel, idx);
+        if (get_dimension(neighbor.get_voxel_pool()->species()) > dimension)
+        {
+            continue;
+        }
+        neighbors.push_back(neighbor);
+    }
+
+    const Integer idx(rng()->uniform_int(0, neighbors.size() - 1));
+    const auto neighbor = neighbors.at(idx);
+
+    if (const auto neighbors = interfaces_.find(neighbor))
+    {
+        const auto idx(rng()->uniform_int(0, neighbors->size() - 1));
+        return neighbors->at(idx);
+    }
+
+    return neighbor;
+}
+
 } // namespace spatiocyte
 
 } // namespace ecell4
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index 10b42626c..ea99e2d8d 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -30,18 +30,18 @@ namespace spatiocyte
 
 struct MoleculeInfo
 {
-    const Real radius;
-    const Real D;
-    const std::string loc;
-    const Shape::dimension_kind dimension;
+    Real radius;
+    Real D;
+    std::string loc;
+    Shape::dimension_kind dimension;
 };
+
 class SpatiocyteWorld : public WorldInterface
 {
 public:
     typedef LatticeSpaceVectorImpl default_root_type;
 
     typedef VoxelSpaceBase::coordinate_id_pair_type coordinate_id_pair_type;
-    typedef VoxelSpaceBase::coordinate_type coordinate_type;
 
     typedef boost::shared_ptr<VoxelSpaceBase> space_type;
     typedef std::vector<space_type> space_container_type;
@@ -254,8 +254,7 @@ class SpatiocyteWorld : public WorldInterface
         {
             if (const auto &view = space->find_voxel(pid))
             {
-                return std::make_pair(
-                    pid, gen_particle_from(space, view->species, view->voxel));
+                return std::make_pair(pid, gen_particle_from(space, *view));
             }
         }
         throw "No particle corresponding to a given ParticleID is found.";
@@ -407,6 +406,21 @@ class SpatiocyteWorld : public WorldInterface
         throw "No VoxelPool corresponding to a given Species is found";
     }
 
+    boost::optional<std::pair<space_type, boost::shared_ptr<VoxelPool>>>
+    find_space_and_voxel_pool(const Species &species) const
+    {
+        for (const auto &space : spaces_)
+        {
+            if (space->has_species(species))
+            {
+                const auto voxel_pool = space->find_voxel_pool(species);
+                return std::pair<space_type, boost::shared_ptr<VoxelPool>>(
+                    space, voxel_pool);
+            }
+        }
+        return boost::none;
+    }
+
     bool has_molecule_pool(const Species &species) const
     {
         for (const auto &space : spaces_)
@@ -438,6 +452,21 @@ class SpatiocyteWorld : public WorldInterface
         throw "No MoleculePool corresponding to a given Species is found";
     }
 
+    boost::optional<std::pair<space_type, boost::shared_ptr<MoleculePool>>>
+    find_space_and_molecule_pool(const Species &species) const
+    {
+        for (const auto &space : spaces_)
+        {
+            if (space->has_molecule_pool(species))
+            {
+                const auto molecule_pool(space->find_molecule_pool(species));
+                return std::pair<space_type, boost::shared_ptr<MoleculePool>>(
+                    space, molecule_pool);
+            }
+        }
+        return boost::none;
+    }
+
     /*
      * Coordinate Transformation
      */
@@ -513,6 +542,16 @@ class SpatiocyteWorld : public WorldInterface
      * SpatiocyteWorld API
      */
 
+    const MoleculeInfo get_molecule_info(const Species &sp) const
+    {
+        const auto itr = molecule_info_cache_.find(sp);
+        if (itr != molecule_info_cache_.end())
+        {
+            return itr->second;
+        }
+        throw NotFound("MoleculeInfo not found");
+    }
+
     /**
      * draw attributes of species and return it as a molecule info.
      * @param sp a species
@@ -520,11 +559,13 @@ class SpatiocyteWorld : public WorldInterface
      */
     const MoleculeInfo get_molecule_info(const Species &sp)
     {
-        boost::optional<Real> diffusion_coef;
-        boost::optional<Real> radius;
-        boost::optional<std::string> location;
-        Shape::dimension_kind dimension =
-            Shape::THREE; // The default dimension is 3
+        // Default
+        MoleculeInfo info = {
+            /* radius = */ voxel_radius(),
+            /* D = */ 0.0,
+            /* loc = */ "",
+            /* dimension = */ Shape::THREE,
+        };
 
         const auto itr = molecule_info_cache_.find(sp);
         if (itr != molecule_info_cache_.end())
@@ -532,10 +573,7 @@ class SpatiocyteWorld : public WorldInterface
             // return itr->second;
             // TODO: the below code is only for warning.
             //       In the future, the value should be returned immediately.
-            diffusion_coef = itr->second.D;
-            radius = itr->second.radius;
-            location = itr->second.loc;
-            dimension = itr->second.dimension;
+            info = itr->second;
         }
         else if (const auto model = lock_model())
         {
@@ -543,70 +581,54 @@ class SpatiocyteWorld : public WorldInterface
 
             if (species_from_model.has_attribute("D"))
             {
-                diffusion_coef = species_from_model.get_attribute_as<Real>("D");
+                info.D = species_from_model.get_attribute_as<Real>("D");
             }
 
             if (species_from_model.has_attribute("radius"))
             {
-                radius = species_from_model.get_attribute_as<Real>("radius");
+                info.radius =
+                    species_from_model.get_attribute_as<Real>("radius");
             }
 
             if (species_from_model.has_attribute("location"))
             {
-                location = species_from_model.get_attribute_as<std::string>(
+                info.loc = species_from_model.get_attribute_as<std::string>(
                     "location");
             }
 
-            dimension = extras::get_dimension_from_model(sp, model);
+            info.dimension = extras::get_dimension_from_model(sp, model);
         }
 
         if (sp.has_attribute("D"))
         {
             const auto new_value = sp.get_attribute_as<Real>("D");
-            if (diffusion_coef && *diffusion_coef != new_value)
+            if (info.D != new_value)
             {
                 warning("D");
-                diffusion_coef = new_value;
+                info.D = new_value;
             }
         }
 
         if (sp.has_attribute("radius"))
         {
             const auto new_value = sp.get_attribute_as<Real>("radius");
-            if (radius && *radius != new_value)
+            if (info.radius != new_value)
             {
                 warning("radius");
-                radius = new_value;
+                info.radius = new_value;
             }
         }
 
         if (sp.has_attribute("location"))
         {
             const auto new_value = sp.get_attribute_as<std::string>("location");
-            if (location && *location != new_value)
+            if (info.loc != new_value)
             {
                 warning("location");
-                location = new_value;
+                info.loc = new_value;
             }
         }
 
-        if (diffusion_coef == boost::none)
-        {
-            diffusion_coef = 0.0;
-        }
-
-        if (radius == boost::none)
-        {
-            radius = voxel_radius();
-        }
-
-        if (location == boost::none)
-        {
-            location = "";
-        }
-
-        const MoleculeInfo info = {*radius, *diffusion_coef, *location,
-                                   dimension};
         molecule_info_cache_.insert(
             molecule_info_cache_t::value_type(sp, info));
         return info;
@@ -685,7 +707,6 @@ class SpatiocyteWorld : public WorldInterface
 
     std::vector<Species> list_non_structure_species() const;
     std::vector<Species> list_structure_species() const;
-    // std::vector<coordinate_type> list_coords(const Species& sp) const;
 
     boost::optional<ParticleID> new_particle(const Species &sp,
                                              const Voxel &voxel)
@@ -735,11 +756,23 @@ class SpatiocyteWorld : public WorldInterface
     boost::optional<Voxel> check_neighbor(const Voxel &voxel,
                                           const std::string &loc);
 
-    // bool update_molecule(coordinate_type at, Species species);
+    const Integer num_neighbors(const Voxel &voxel) const
+    {
+        return voxel.space.lock()->num_neighbors(voxel.coordinate);
+    }
+
+    const Voxel get_neighbor(const Voxel &voxel, Integer nrand) const
+    {
+        return Voxel(voxel.space,
+                     voxel.space.lock()->get_neighbor(voxel.coordinate, nrand));
+    }
+
+    const Voxel get_neighbor_randomly(const Voxel &voxel,
+                                      Shape::dimension_kind dimension);
 
     const Species &draw_species(const Species &pttrn) const;
 
-    boost::shared_ptr<RandomNumberGenerator> rng() { return rng_; }
+    boost::shared_ptr<RandomNumberGenerator> rng() const { return rng_; }
 
     void bind_to(boost::shared_ptr<Model> model)
     {
@@ -775,6 +808,11 @@ class SpatiocyteWorld : public WorldInterface
         return get_molecule_info(species).dimension;
     }
 
+    Shape::dimension_kind get_dimension(const Species &species) const
+    {
+        return get_molecule_info(species).dimension;
+    }
+
     /**
      * static members
      */
@@ -803,25 +841,8 @@ class SpatiocyteWorld : public WorldInterface
     }
 
 protected:
-    bool is_inside(const coordinate_type &coord) const
-    {
-        return get_root()->is_inside(coord);
-    }
-
     space_type get_root() const { return spaces_.at(0); }
 
-    space_type get_space(coordinate_type &coordinate)
-    {
-        for (const auto &space : spaces_)
-        {
-            if (coordinate < space->size())
-                return space;
-
-            coordinate -= space->size();
-        }
-        return space_type();
-    }
-
     Integer add_structure2(const Species &sp, const std::string &location,
                            const boost::shared_ptr<const Shape> shape);
     Integer add_structure3(const Species &sp, const std::string &location,
@@ -829,10 +850,12 @@ class SpatiocyteWorld : public WorldInterface
     bool is_surface_voxel(const Voxel &voxel,
                           const boost::shared_ptr<const Shape> shape) const;
 
-    Particle gen_particle_from(const space_type &space, const Species &species,
-                               const coordinate_type coordinate) const
+    Particle gen_particle_from(const space_type &space,
+                               const VoxelView &view) const
     {
-        const auto position = space->coordinate2position(coordinate);
+        const auto species = view.species;
+
+        const auto position = space->coordinate2position(view.voxel);
         const auto minfo_iter = molecule_info_cache_.find(species);
         if (minfo_iter != molecule_info_cache_.end())
         {
@@ -869,9 +892,7 @@ class SpatiocyteWorld : public WorldInterface
     {
         return map_voxels<std::pair<ParticleID, Particle>>(
             list_fn, [this](const space_type &space, const VoxelView &view) {
-                return std::make_pair(
-                    view.pid,
-                    gen_particle_from(space, view.species, view.voxel));
+                return std::make_pair(view.pid, gen_particle_from(space, view));
             });
     }
 
@@ -879,20 +900,12 @@ class SpatiocyteWorld : public WorldInterface
     std::vector<ParticleBase<Voxel>> list_voxels_private(ListFn list_fn) const
     {
         return map_voxels<ParticleBase<Voxel>>(
-            list_fn, [](const space_type &space, const VoxelView &view) {
+            list_fn, [this](const space_type &space, const VoxelView &view) {
                 return ParticleBase<Voxel>(view.pid, view.species,
                                            Voxel(space, view.voxel));
             });
     }
 
-public:
-    // TODO: Calling this function is invalid, and this should be removed.
-    Voxel coordinate2voxel(const coordinate_type &coordinate)
-    {
-        coordinate_type coord(coordinate);
-        return Voxel(get_space(coord), coord);
-    }
-
 protected:
     typedef utils::get_mapper_mf<Species, MoleculeInfo>::type
         molecule_info_cache_t;
@@ -900,15 +913,15 @@ class SpatiocyteWorld : public WorldInterface
     std::size_t size_;
     space_container_type spaces_;
 
-    OneToManyMap<coordinate_type> interfaces_;
-    OneToManyMap<coordinate_type> neighbors_;
+    OneToManyMap<Voxel> interfaces_;
+    OneToManyMap<Voxel> neighbors_;
 
     boost::shared_ptr<RandomNumberGenerator> rng_;
     SerialIDGenerator<ParticleID> sidgen_;
 
     boost::weak_ptr<Model> model_;
     molecule_info_cache_t molecule_info_cache_;
-};
+}; // namespace spatiocyte
 
 inline SpatiocyteWorld *create_spatiocyte_world_cell_list_impl(
     const Real3 &edge_lengths, const Real &voxel_radius,
diff --git a/ecell4/spatiocyte/StepEvent.cpp b/ecell4/spatiocyte/StepEvent.cpp
index 834020be4..585e55b2d 100644
--- a/ecell4/spatiocyte/StepEvent.cpp
+++ b/ecell4/spatiocyte/StepEvent.cpp
@@ -10,9 +10,19 @@ namespace spatiocyte
 StepEvent::StepEvent(boost::shared_ptr<Model> model,
                      boost::shared_ptr<SpatiocyteWorld> world,
                      const Species &species, const Real &t, const Real alpha)
-    : SpatiocyteEvent(t), model_(model), world_(world),
-      mpool_(world_->find_molecule_pool(species)), alpha_(alpha)
+    : SpatiocyteEvent(t), model_(model), world_(world), alpha_(alpha)
 {
+    if (const auto space_and_molecule_pool =
+            world_->find_space_and_molecule_pool(species))
+    {
+        space_ = space_and_molecule_pool->first;
+        mpool_ = space_and_molecule_pool->second;
+    }
+    else
+    {
+        throw "MoleculePool is not found";
+    }
+
     time_ = t;
 }
 
@@ -48,8 +58,7 @@ void StepEvent3D::walk(const Real &alpha)
     std::size_t idx(0);
     for (const auto &info : voxels)
     {
-        const Voxel voxel(world_->coordinate2voxel(info.coordinate));
-        const Integer rnd(rng->uniform_int(0, voxel.num_neighbors() - 1));
+        const Voxel voxel(space_, info.coordinate);
 
         if (voxel.get_voxel_pool() != mpool_)
         {
@@ -58,7 +67,8 @@ void StepEvent3D::walk(const Real &alpha)
             continue;
         }
 
-        const Voxel neighbor(voxel.get_neighbor(rnd));
+        const Voxel neighbor(
+            world_->get_neighbor_randomly(voxel, Shape::THREE));
 
         if (world_->can_move(voxel, neighbor))
         {
@@ -90,10 +100,6 @@ StepEvent2D::StepEvent2D(boost::shared_ptr<Model> model,
         dt_ = R * R / D * alpha_;
 
     time_ = t + dt_;
-
-    nids_.clear();
-    for (unsigned int i(0); i < 12; ++i)
-        nids_.push_back(i);
 }
 
 void StepEvent2D::walk(const Real &alpha)
@@ -110,8 +116,7 @@ void StepEvent2D::walk(const Real &alpha)
     std::size_t idx(0);
     for (const auto &info : voxels)
     {
-        // TODO: Calling coordinate2voxel is invalid
-        const Voxel voxel(world_->coordinate2voxel(info.coordinate));
+        const Voxel voxel(space_, info.coordinate);
 
         if (voxel.get_voxel_pool() != mpool_)
         {
@@ -120,32 +125,18 @@ void StepEvent2D::walk(const Real &alpha)
             continue;
         }
 
-        const std::size_t num_neighbors(voxel.num_neighbors());
+        const Voxel neighbor(world_->get_neighbor_randomly(voxel, Shape::TWO));
 
-        ecell4::shuffle(*(rng.get()), nids_);
-        for (const auto &index : nids_)
+        if (world_->can_move(voxel, neighbor))
         {
-            if (index >= num_neighbors)
-                continue;
-
-            const Voxel neighbor(voxel.get_neighbor(index));
-            boost::shared_ptr<const VoxelPool> target(
-                neighbor.get_voxel_pool());
-
-            if (world_->get_dimension(target->species()) > Shape::TWO)
-                continue;
-
-            if (world_->can_move(voxel, neighbor))
-            {
-                if (rng->uniform(0, 1) <= alpha)
-                    world_->move(voxel, neighbor, /*candidate=*/idx);
-            }
-            else
-            {
-                attempt_reaction_(info, neighbor, alpha);
-            }
-            break;
+            if (rng->uniform(0, 1) <= alpha)
+                world_->move(voxel, neighbor, /*candidate=*/idx);
         }
+        else
+        {
+            attempt_reaction_(info, neighbor, alpha);
+        }
+
         ++idx;
     }
 }
@@ -154,16 +145,10 @@ void StepEvent::attempt_reaction_(
     const SpatiocyteWorld::coordinate_id_pair_type &info, const Voxel &dst,
     const Real &alpha)
 {
-    // TODO: Calling coordiante2voxel is invalid
-    const Voxel voxel(world_->coordinate2voxel(info.coordinate));
+    const Voxel voxel(space_, info.coordinate);
     boost::shared_ptr<const VoxelPool> from_mt(voxel.get_voxel_pool());
     boost::shared_ptr<const VoxelPool> to_mt(dst.get_voxel_pool());
 
-    if (to_mt->is_vacant())
-    {
-        return;
-    }
-
     const Species &speciesA(from_mt->species());
     const Species &speciesB(to_mt->species());
 
diff --git a/ecell4/spatiocyte/Voxel.hpp b/ecell4/spatiocyte/Voxel.hpp
index b75f1ffd6..e4ad04602 100644
--- a/ecell4/spatiocyte/Voxel.hpp
+++ b/ecell4/spatiocyte/Voxel.hpp
@@ -3,6 +3,7 @@
 
 #include <boost/weak_ptr.hpp>
 #include <ecell4/core/VoxelSpaceBase.hpp>
+#include <ecell4/core/hash.hpp>
 #include <ecell4/core/types.hpp>
 
 namespace ecell4
@@ -11,9 +12,11 @@ namespace ecell4
 namespace spatiocyte
 {
 
+class SpatiocyteWorld;
+
 struct Voxel
 {
-    typedef Integer coordinate_type;
+    typedef VoxelSpaceBase::coordinate_type coordinate_type;
 
     Voxel(boost::weak_ptr<VoxelSpaceBase> space, coordinate_type coordinate)
         : space(space), coordinate(coordinate)
@@ -36,14 +39,9 @@ struct Voxel
         return space.lock()->get_voxel_pool_at(coordinate);
     }
 
-    Integer num_neighbors() const
+    bool operator==(const Voxel &rhs) const noexcept
     {
-        return space.lock()->num_neighbors(coordinate);
-    }
-
-    Voxel get_neighbor(Integer nrand) const
-    {
-        return Voxel(space, space.lock()->get_neighbor(coordinate, nrand));
+        return space.lock() == rhs.space.lock() && coordinate == rhs.coordinate;
     }
 };
 
@@ -51,4 +49,17 @@ struct Voxel
 
 } // namespace ecell4
 
+ECELL4_DEFINE_HASH_BEGIN()
+template <>
+struct hash<ecell4::spatiocyte::Voxel>
+{
+    std::size_t operator()(const ecell4::spatiocyte::Voxel &val) const
+    {
+        auto ptr = val.space.lock().get();
+        return hash<decltype(ptr)>()(ptr) ^
+               static_cast<std::size_t>(val.coordinate);
+    }
+};
+ECELL4_DEFINE_HASH_END()
+
 #endif
diff --git a/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp b/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
index 369aa41cd..6a8c433b6 100644
--- a/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
+++ b/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
@@ -27,7 +27,7 @@ struct Fixture
 
     Fixture()
         : edge_lengths(1e-6, 1e-6, 1e-6), voxel_radius(1e-8),
-          rng(new GSLRandomNumberGenerator()), model(new NetworkModel),
+          rng(new GSLRandomNumberGenerator()), model(new NetworkModel()),
           world(edge_lengths, voxel_radius, rng)
     {
         world.bind_to(model);
@@ -63,6 +63,26 @@ BOOST_AUTO_TEST_CASE(SpatiocyteWorld_test_list_particles)
     BOOST_CHECK_EQUAL(particles.size(), 0);
 }
 
+BOOST_AUTO_TEST_CASE(SpatiocyteWorld_test_get_molecule_info)
+{
+    const Species m("M", voxel_radius, 0.0);
+    const Species a("A", voxel_radius, 1.0, "M");
+    model->add_species_attribute(m);
+    model->add_species_attribute(a);
+
+    const auto info_m = world.get_molecule_info(m);
+    BOOST_CHECK_EQUAL(info_m.radius, voxel_radius);
+    BOOST_CHECK_EQUAL(info_m.D, 0.0);
+    BOOST_CHECK_EQUAL(info_m.loc, "");
+    BOOST_CHECK_EQUAL(info_m.dimension, Shape::THREE);
+
+    const auto info_a = world.get_molecule_info(a);
+    BOOST_CHECK_EQUAL(info_a.radius, voxel_radius);
+    BOOST_CHECK_EQUAL(info_a.D, 1.0);
+    BOOST_CHECK_EQUAL(info_a.loc, "M");
+    BOOST_CHECK_EQUAL(info_a.dimension, Shape::THREE);
+}
+
 BOOST_AUTO_TEST_CASE(SpatiocyteWorld_test_update_particles)
 {
     SerialIDGenerator<ParticleID> sidgen;
@@ -127,9 +147,9 @@ BOOST_AUTO_TEST_CASE(SpatiocyteWorld_test_neighbor)
     const Species sp("TEST", 1e-8, 1e-12);
     model->add_species_attribute(sp);
 
-    for (Integer i(0); i < voxel.num_neighbors(); ++i)
+    for (Integer i(0); i < world.num_neighbors(voxel); ++i)
     {
-        world.new_particle(sp, voxel.get_neighbor(i));
+        world.new_particle(sp, world.get_neighbor(voxel, i));
     }
     std::vector<std::pair<ParticleID, Particle>> particles(
         world.list_particles());
diff --git a/tests/spatiocyte/offlattice.py b/tests/spatiocyte/offlattice.py
index ccfa9d8b9..105aed8c1 100644
--- a/tests/spatiocyte/offlattice.py
+++ b/tests/spatiocyte/offlattice.py
@@ -15,3 +15,18 @@ def test_constructor(self):
         species = Species('Base')
         world = SpatiocyteWorld(ones(), self.voxel_radius)
         world.add_offlattice(species, self.offlattice)
+
+    def test_add_molecules(self):
+        base = Species('Base', radius=self.voxel_radius, D=0.0)
+        species = Species('A', radius=self.voxel_radius, D=1.0, location='Base')
+
+        model = NetworkModel()
+        model.add_species_attribute(base)
+        model.add_species_attribute(species)
+
+        world = SpatiocyteWorld(ones(), self.voxel_radius)
+        world.bind_to(model)
+        world.add_offlattice(base, self.offlattice)
+
+        world.add_molecules(species, 5)
+        self.assertEqual(5, world.num_molecules(species))