diff --git a/neurom/core/morphology.py b/neurom/core/morphology.py
index f5b2146b..d4650611 100644
--- a/neurom/core/morphology.py
+++ b/neurom/core/morphology.py
@@ -46,34 +46,26 @@ class Section:
 
     def __init__(self, morphio_section):
         """The section constructor."""
-        self.morphio_section = morphio_section
+        self._morphio_section = morphio_section
+
+    def to_morphio(self):
+        """Returns the morphio section."""
+        return self._morphio_section
 
     @property
     def id(self):
         """Returns the section ID."""
-        return self.morphio_section.id
+        return self._morphio_section.id
 
     @property
     def parent(self):
         """Returns the parent section if non root section else None."""
-        if self.morphio_section.is_root:
-            return None
-        return Section(self.morphio_section.parent)
+        return None if self.is_root() else Section(self._morphio_section.parent)
 
     @property
     def children(self):
         """Returns a list of child section."""
-        return [Section(child) for child in self.morphio_section.children]
-
-    def append_section(self, section):
-        """Appends a section to the current section object.
-
-        Args:
-            section (morphio.Section|morphio.mut.Section|Section|morphio.PointLevel): a section
-        """
-        if isinstance(section, Section):
-            return self.morphio_section.append_section(section.morphio_section)
-        return self.morphio_section.append_section(section)
+        return [Section(child) for child in self._morphio_section.children]
 
     def is_homogeneous_point(self):
         """A section is homogeneous if it has the same type with its children."""
@@ -93,7 +85,7 @@ def is_leaf(self):
 
     def is_root(self):
         """Is tree the root node?"""
-        return self.parent is None
+        return self._morphio_section.is_root
 
     def ipreorder(self):
         """Depth-first pre-order iteration of tree nodes."""
@@ -124,10 +116,10 @@ def iupstream(self, stop_node=None):
         """
         if stop_node is None:
             def stop_condition(section):
-                return section.parent is None
+                return section.is_root()
         else:
             def stop_condition(section):
-                return section == stop_node
+                return section.is_root() or section == stop_node
 
         current_section = self
         while not stop_condition(current_section):
@@ -157,35 +149,23 @@ def ibifurcation_point(self, iter_mode=ipreorder):
 
     def __eq__(self, other):
         """Equal when its morphio section is equal."""
-        return self.morphio_section == other.morphio_section
+        return self.to_morphio().has_same_shape(other.to_morphio())
 
     def __hash__(self):
         """Hash of its id."""
         return self.id
 
-    def __nonzero__(self):
-        """If has children."""
-        return self.morphio_section is not None
-
-    __bool__ = __nonzero__
-
     @property
     def points(self):
         """Returns the section list of points the NeuroM way (points + radius)."""
-        return np.concatenate((self.morphio_section.points,
-                               self.morphio_section.diameters[:, np.newaxis] / 2.),
+        return np.concatenate((self._morphio_section.points,
+                               self._morphio_section.diameters[:, np.newaxis] / 2.),
                               axis=1)
 
-    @points.setter
-    def points(self, value):
-        """Set the points."""
-        self.morphio_section.points = np.copy(value[:, COLS.XYZ])
-        self.morphio_section.diameters = np.copy(value[:, COLS.R]) * 2
-
     @property
     def type(self):
         """Returns the section type."""
-        return NeuriteType(int(self.morphio_section.type))
+        return NeuriteType(int(self._morphio_section.type))
 
     @property
     def length(self):
@@ -234,11 +214,11 @@ def _homogeneous_subtrees(neurite):
     sub-tree.
     """
     it = neurite.root_node.ipreorder()
-    homogeneous_neurites = [Neurite(next(it).morphio_section)]
+    homogeneous_neurites = [Neurite(next(it).to_morphio())]
 
     for section in it:
         if section.type != section.parent.type:
-            homogeneous_neurites.append(Neurite(section.morphio_section))
+            homogeneous_neurites.append(Neurite(section.to_morphio()))
 
     homogeneous_types = [neurite.type for neurite in homogeneous_neurites]
 
@@ -416,17 +396,10 @@ def graft_morphology(section):
     """Returns a morphology starting at section."""
     assert isinstance(section, Section)
     m = morphio.mut.Morphology()
-    m.append_root_section(section.morphio_section)
+    m.append_root_section(section.to_morphio())
     return Morphology(m)
 
 
-def graft_neuron(section):
-    """Deprecated in favor of ``graft_morphology``."""
-    warn_deprecated('`neurom.core.neuron.graft_neuron` is deprecated in favor of '
-                    '`neurom.core.morphology.graft_morphology`')  # pragma: no cover
-    return graft_morphology(section)  # pragma: no cover
-
-
 class Neurite:
     """Class representing a neurite tree."""
 
@@ -436,7 +409,12 @@ def __init__(self, root_node):
         Args:
             root_node (morphio.Section): root section
         """
-        self.morphio_root_node = root_node
+        self._root_node = root_node
+
+    @property
+    def morphio_root_node(self):
+        """Returns the morphio root section."""
+        return self._root_node
 
     @property
     def root_node(self):
@@ -511,10 +489,6 @@ def iter_sections(self, order=Section.ipreorder, neurite_order=NeuriteIter.FileO
         """
         return iter_sections(self, iterator_type=order, neurite_order=neurite_order)
 
-    def __nonzero__(self):
-        """If has root node."""
-        return bool(self.morphio_root_node)
-
     def __eq__(self, other):
         """If root node ids and types are equal."""
         return self.type == other.type and self.morphio_root_node.id == other.morphio_root_node.id
@@ -523,37 +497,37 @@ def __hash__(self):
         """Hash is made of tuple of type and root_node."""
         return hash((self.type, self.root_node))
 
-    __bool__ = __nonzero__
-
     def __repr__(self):
         """Return a string representation."""
         return 'Neurite <type: %s>' % self.type
 
 
-class Morphology(morphio.mut.Morphology):
+class Morphology:
     """Class representing a simple morphology."""
 
     def __init__(self, filename, name=None):
         """Morphology constructor.
 
         Args:
-            filename (str|Path): a filename
-            name (str): a option morphology name
+            filename (str|Path): a filename or morphio.{mut}.Morphology object
+            name (str): an optional morphology name
         """
-        super().__init__(filename)
+        self._morphio_morph = morphio.mut.Morphology(filename).as_immutable()
         self.name = name if name else 'Morphology'
-        self.morphio_soma = super().soma
-        self.neurom_soma = make_soma(self.morphio_soma)
+        self.soma = make_soma(self._morphio_morph.soma)
 
-    @property
-    def soma(self):
-        """Corresponding soma."""
-        return self.neurom_soma
+    def to_morphio(self):
+        """Returns the morphio morphology object."""
+        return self._morphio_morph
 
     @property
     def neurites(self):
         """The list of neurites."""
-        return [Neurite(root_section) for root_section in self.root_sections]
+        return [Neurite(root_section) for root_section in self._morphio_morph.root_sections]
+
+    def section(self, section_id):
+        """Returns the section with the given id."""
+        return Section(self._morphio_morph.section(section_id))
 
     @property
     def sections(self):
@@ -568,21 +542,22 @@ def points(self):
 
     def transform(self, trans):
         """Return a copy of this morphology with a 3D transformation applied."""
-        obj = Morphology(self)
-        obj.morphio_soma.points = trans(obj.morphio_soma.points)
+        mut = self._morphio_morph.as_mutable()
+        mut.soma.points = trans(mut.soma.points)
+
+        for section in mut.iter():
+            section.points = trans(section.points)
 
-        for section in obj.sections:
-            section.morphio_section.points = trans(section.morphio_section.points)
-        return obj
+        return Morphology(mut)
 
     def __copy__(self):
         """Creates a deep copy of Morphology instance."""
-        return Morphology(self, self.name)
+        return Morphology(self.to_morphio(), self.name)
 
     def __deepcopy__(self, memodict={}):
         """Creates a deep copy of Morphology instance."""
         # pylint: disable=dangerous-default-value
-        return Morphology(self, self.name)
+        return Morphology(self.to_morphio(), self.name)
 
     def __repr__(self):
         """Return a string representation."""
diff --git a/neurom/core/soma.py b/neurom/core/soma.py
index 9c6efaea..b5c0b4d6 100755
--- a/neurom/core/soma.py
+++ b/neurom/core/soma.py
@@ -71,13 +71,6 @@ def points(self):
                                self._morphio_soma.diameters[:, np.newaxis] / 2.),
                               axis=1)
 
-    @points.setter
-    def points(self, values):
-        """Set the points."""
-        values = np.asarray(values)
-        self._morphio_soma.points = np.copy(values[:, COLS.XYZ])
-        self._morphio_soma.diameters = np.copy(values[:, COLS.R]) * 2
-
     @property
     def volume(self):
         """Gets soma volume assuming it is a sphere."""
diff --git a/neurom/io/utils.py b/neurom/io/utils.py
index 393d4f64..000e6e74 100644
--- a/neurom/io/utils.py
+++ b/neurom/io/utils.py
@@ -156,7 +156,10 @@ def load_morphology(morph, reader=None):
                                                    )
                                                    )'''), reader='asc')
     """
-    if isinstance(morph, (Morphology, morphio.Morphology, morphio.mut.Morphology)):
+    if isinstance(morph, Morphology):
+        return Morphology(morph.to_morphio())
+
+    if isinstance(morph, (morphio.Morphology, morphio.mut.Morphology)):
         return Morphology(morph)
 
     if reader:
diff --git a/pyproject.toml b/pyproject.toml
index f6c16894..ca242af0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,3 +4,8 @@ requires = [
     "wheel",
 ]
 build-backend = "setuptools.build_meta"
+
+[tool.pytest.ini_options]
+testpaths = [
+    "tests",
+]
diff --git a/tests/core/test_iter.py b/tests/core/test_iter.py
index 01f3bc40..92f416be 100644
--- a/tests/core/test_iter.py
+++ b/tests/core/test_iter.py
@@ -98,10 +98,8 @@ def test_iter_population():
 
 
 def test_iter_sections_default():
-
-    ref = [s for n in POP.neurites for s in n.iter_sections()]
-    assert (ref ==
-                          [n for n in iter_sections(POP)])
+    ref = [s.id for n in POP.neurites for s in n.iter_sections()]
+    assert (ref == [n.id for n in iter_sections(POP)])
 
 def test_iter_sections_default_pop():
     ref = [s.id for n in POP.neurites for s in n.iter_sections()]
diff --git a/tests/core/test_neuron.py b/tests/core/test_neuron.py
index 4e1fbf73..ff9979ec 100644
--- a/tests/core/test_neuron.py
+++ b/tests/core/test_neuron.py
@@ -85,12 +85,6 @@ def test_for_morphio():
                         [1., 1., 1., 0.5],
                         [2., 2., 2., 0.5]])
 
-    neurom_m.soma.points = [[1, 1, 1, 1],
-                                 [2, 2, 2, 2]]
-    assert_array_equal(neurom_m.soma.points,
-                       [[1, 1, 1, 1],
-                        [2, 2, 2, 2]])
-
 
 def _check_cloned_morphology(m, m2):
     # check if two morphs are identical
@@ -116,10 +110,6 @@ def _check_cloned_morphology(m, m2):
     for neu1, neu2 in zip(m.neurites, m2.neurites):
         assert neu1 is not neu2
 
-    # check if changes are propagated between morphs
-    m2.soma.radius = 10.
-    assert m.soma.radius != m2.soma.radius
-
 
 def test_copy():
     m = nm.load_morphology(SWC_PATH / 'simple.swc')
diff --git a/tests/core/test_section.py b/tests/core/test_section.py
index 93708504..1e86df4b 100644
--- a/tests/core/test_section.py
+++ b/tests/core/test_section.py
@@ -48,10 +48,14 @@ def test_section_base_func():
     # __nonzero__
     assert section
 
+
 def test_section_tree():
     m = nm.load_morphology(str(SWC_PATH / 'simple.swc'))
 
     assert m.sections[0].parent is None
+
+    assert m.sections[0] == m.sections[0]
+
     assert m.sections[0] == m.sections[0].children[0].parent
 
     assert_array_equal([s.is_root() for s in m.sections],
@@ -70,6 +74,11 @@ def test_section_tree():
                        [0])
     assert_array_equal([s.id for s in m.sections[2].iupstream()],
                        [2, 0])
+    assert_array_equal([s.id for s in m.sections[2].iupstream(stop_node=m.sections[2])],
+                       [2])
+    # if a stop node that is not upstream is given, it should stop at root
+    assert_array_equal([s.id for s in m.sections[2].iupstream(stop_node=m.sections[1])],
+                       [2, 0])
     assert_array_equal([s.id for s in m.neurites[0].root_node.ileaf()],
                        [1, 2])
     assert_array_equal([s.id for s in m.sections[2].ileaf()],
@@ -78,31 +87,3 @@ def test_section_tree():
                        [0])
     assert_array_equal([s.id for s in m.neurites[0].root_node.ibifurcation_point()],
                        [0])
-
-
-def test_append_section():
-    n = nm.load_morphology(SWC_PATH / 'simple.swc')
-    s = n.sections[0]
-
-    s.append_section(n.sections[-1])
-    assert len(s.children) == 3
-    assert s.children[-1].id == 6
-    assert s.children[-1].type == n.sections[-1].type
-
-    s.append_section(n.sections[-1].morphio_section)
-    assert len(s.children) == 4
-    assert s.children[-1].id == 7
-    assert s.children[-1].type == n.sections[-1].type
-
-
-def test_set_points():
-    n = nm.load_morphology(SWC_PATH / 'simple.swc')
-    s = n.sections[0]
-    s.points = np.array([
-        [0, 5, 0, 2],
-        [0, 7, 0, 2],
-    ])
-    assert_array_equal(s.points, np.array([
-        [0, 5, 0, 2],
-        [0, 7, 0, 2],
-    ]))
diff --git a/tests/features/test_morphology.py b/tests/features/test_morphology.py
index 2337c839..cf755a10 100644
--- a/tests/features/test_morphology.py
+++ b/tests/features/test_morphology.py
@@ -62,12 +62,15 @@
 
 def _add_neurite_trunk(morph, elevation, azimuth, neurite_type=SectionType.basal_dendrite):
     """Add a neurite from the elevation and azimuth to a given morphology."""
+    mut = morph.to_morphio().as_mutable()
     new_pts = np.array(
         morphmath.vector_from_spherical(elevation, azimuth),
         ndmin=2
     )
+
     point_lvl = PointLevel(new_pts, [1])
-    morph.append_root_section(point_lvl, neurite_type)
+    mut.append_root_section(point_lvl, neurite_type)
+    return Morphology(mut)
 
 
 def test_soma_volume():
@@ -154,14 +157,23 @@ def test_trunk_section_lengths():
 
 
 def test_trunk_origin_radii():
-    morph = Morphology(SIMPLE)
-    morph.section(0).diameters = [2, 1]
-    morph.section(3).diameters = [2, 0.5]
-
+    morph = load_swc(
+    """
+     1 1  0  0 0 1. -1
+     2 3  0  0 0 1.0  1
+     3 3  0  5 0 0.5  2
+     4 3 -5  5 0 0.  3
+     5 3  6  5 0 0.  3
+     6 2  0  0 0 1.0  1
+     7 2  0 -4 0 0.25  6
+     8 2  6 -4 0 0.  7
+     9 2 -5 -4 0 0.  7
+    """
+    )
     ret = morphology.trunk_origin_radii(morph)
     assert ret == [1.0, 1.0]
 
-    ret = morphology.trunk_origin_radii(morph, min_length_filter=1)
+    ret = morphology.trunk_origin_radii(morph, min_length_filter=1.0)
     assert_array_almost_equal(ret, [0.5, 0.25])
 
     with pytest.warns(
@@ -259,8 +271,8 @@ def test_trunk_angles():
     morph = load_morphology(SWC_PATH / 'simple_trunk.swc')
 
     # Add two basals
-    _add_neurite_trunk(morph, np.pi / 3, np.pi / 4)
-    _add_neurite_trunk(morph, -np.pi / 3, -np.pi / 4)
+    morph = _add_neurite_trunk(morph, np.pi / 3, np.pi / 4)
+    morph = _add_neurite_trunk(morph, -np.pi / 3, -np.pi / 4)
 
     ret = morphology.trunk_angles(morph)
     assert_array_almost_equal(ret, [np.pi / 2, 0.387596, 1.183199, 1.183199, 0.387596, np.pi / 2])
@@ -302,8 +314,8 @@ def test_trunk_angles_inter_types():
     morph = load_morphology(SWC_PATH / 'simple_trunk.swc')
 
     # Add two basals
-    _add_neurite_trunk(morph, np.pi / 3, np.pi / 4)
-    _add_neurite_trunk(morph, -np.pi / 3, -np.pi / 4)
+    morph = _add_neurite_trunk(morph, np.pi / 3, np.pi / 4)
+    morph = _add_neurite_trunk(morph, -np.pi / 3, -np.pi / 4)
 
     # Test with no source
     ret = morphology.trunk_angles_inter_types(
@@ -378,8 +390,8 @@ def test_trunk_angles_from_vector():
     morph = load_morphology(SWC_PATH / 'simple_trunk.swc')
 
     # Add two basals
-    _add_neurite_trunk(morph, np.pi / 3, np.pi / 4)
-    _add_neurite_trunk(morph, -np.pi / 3, -np.pi / 4)
+    morph = _add_neurite_trunk(morph, np.pi / 3, np.pi / 4)
+    morph = _add_neurite_trunk(morph, -np.pi / 3, -np.pi / 4)
 
     # Test with no neurite selected
     ret = morphology.trunk_angles_from_vector(