diff --git a/binds/python/CMakeLists.txt b/binds/python/CMakeLists.txt
index 35f82eb97..6b8fc3904 100644
--- a/binds/python/CMakeLists.txt
+++ b/binds/python/CMakeLists.txt
@@ -7,6 +7,7 @@ pybind11_add_module(_morphio SYSTEM
     bind_immutable.cpp
     bindings_utils.cpp
     bind_misc.cpp
+    bind_soma.cpp
     bind_mutable.cpp
     bind_vasculature.cpp
 )
diff --git a/binds/python/bind_immutable.cpp b/binds/python/bind_immutable.cpp
index 7c44ab9ee..d02463ca0 100644
--- a/binds/python/bind_immutable.cpp
+++ b/binds/python/bind_immutable.cpp
@@ -213,33 +213,6 @@ void bind_immutable_module(py::module& m) {
                                &morphio::EndoplasmicReticulum::filamentCounts,
                                "Returns the number of filaments for each neuronal section");
 
-
-    py::class_<morphio::Soma>(m, "Soma")
-        .def(py::init<const morphio::Soma&>())
-        .def_property_readonly(
-            "points",
-            [](morphio::Soma* soma) { return span_array_to_ndarray(soma->points()); },
-            "Returns the coordinates (x,y,z) of all soma point")
-        .def_property_readonly(
-            "diameters",
-            [](morphio::Soma* soma) { return span_to_ndarray(soma->diameters()); },
-            "Returns the diameters of all soma points")
-
-        .def_property_readonly(
-            "center",
-            [](morphio::Soma* soma) { return py::array(3, soma->center().data()); },
-            "Returns the center of gravity of the soma points")
-        .def_property_readonly("max_distance",
-                               &morphio::Soma::maxDistance,
-                               "Return the maximum distance between the center of gravity "
-                               "and any of the soma points")
-        .def_property_readonly("type", &morphio::Soma::type, "Returns the soma type")
-
-        .def_property_readonly("surface",
-                               &morphio::Soma::surface,
-                               "Returns the soma surface\n\n"
-                               "Note: the soma surface computation depends on the soma type");
-
     py::class_<morphio::Section>(m, "Section")
         .def("__str__",
              [](const morphio::Section& section) {
diff --git a/binds/python/bind_mutable.cpp b/binds/python/bind_mutable.cpp
index fd20097c6..8105c4e80 100644
--- a/binds/python/bind_mutable.cpp
+++ b/binds/python/bind_mutable.cpp
@@ -54,8 +54,7 @@ void bind_mutable_module(py::module& m) {
                                py::return_value_policy::reference)
         .def_property_readonly(
             "soma",
-            static_cast<std::shared_ptr<morphio::mut::Soma>& (morphio::mut::Morphology::*) ()>(
-                &morphio::mut::Morphology::soma),
+            [](const morphio::mut::Morphology& self) { return self.soma(); },
             "Returns a reference to the soma object\n\n"
             "Note: multiple morphologies can share the same Soma "
             "instance")
@@ -132,7 +131,7 @@ void bind_mutable_module(py::module& m) {
         .def_property_readonly("version", &morphio::mut::Morphology::version, "Returns the version")
 
         .def("remove_unifurcations",
-             static_cast<void (morphio::mut::Morphology::*) ()>(
+             static_cast<void (morphio::mut::Morphology::*)()>(
                  &morphio::mut::Morphology::removeUnifurcations),
              "Fixes the morphology single child sections and issues warnings"
              "if the section starts and ends are inconsistent")
@@ -455,42 +454,6 @@ void bind_mutable_module(py::module& m) {
              "point_level_properties"_a,
              "section_type"_a = morphio::SectionType::SECTION_UNDEFINED);
 
-    py::class_<morphio::mut::Soma, std::shared_ptr<morphio::mut::Soma>>(m, "Soma")
-        .def(py::init<const morphio::Property::PointLevel&>())
-        .def_property(
-            "points",
-            [](morphio::mut::Soma* soma) {
-                return py::array(static_cast<py::ssize_t>(soma->points().size()),
-                                 soma->points().data());
-            },
-            [](morphio::mut::Soma* soma, py::array_t<morphio::floatType> _points) {
-                soma->points() = array_to_points(_points);
-            },
-            "Returns the coordinates (x,y,z) of all soma point")
-        .def_property(
-            "diameters",
-            [](morphio::mut::Soma* soma) {
-                return py::array(static_cast<py::ssize_t>(soma->diameters().size()),
-                                 soma->diameters().data());
-            },
-            [](morphio::mut::Soma* soma, py::array_t<morphio::floatType> _diameters) {
-                soma->diameters() = _diameters.cast<std::vector<morphio::floatType>>();
-            },
-            "Returns the diameters of all soma points")
-        .def_property_readonly("type", &morphio::mut::Soma::type, "Returns the soma type")
-        .def_property_readonly("surface",
-                               &morphio::mut::Soma::surface,
-                               "Returns the soma surface\n\n"
-                               "Note: the soma surface computation depends on the soma type")
-        .def_property_readonly("max_distance",
-                               &morphio::mut::Soma::maxDistance,
-                               "Return the maximum distance between the center of gravity "
-                               "and any of the soma points")
-        .def_property_readonly(
-            "center",
-            [](morphio::mut::Soma* soma) { return py::array(3, soma->center().data()); },
-            "Returns the center of gravity of the soma points");
-
     py::class_<morphio::mut::EndoplasmicReticulum>(m, "EndoplasmicReticulum")
         .def(py::init<>())
         .def(py::init<const std::vector<uint32_t>&,
diff --git a/binds/python/bind_soma.cpp b/binds/python/bind_soma.cpp
new file mode 100644
index 000000000..2daa69c0c
--- /dev/null
+++ b/binds/python/bind_soma.cpp
@@ -0,0 +1,58 @@
+
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <morphio/soma.h>
+#include <morphio/types.h>
+
+#include "bind_soma.h"
+#include "bindings_utils.h"
+
+namespace py = pybind11;
+
+void bind_soma_module(py::module& m) {
+    py::class_<morphio::Soma, std::shared_ptr<morphio::Soma>>(m, "Soma")
+        .def(py::init<const morphio::Soma&>())
+
+        .def(py::init<const morphio::Property::PointLevel&>())
+
+        .def_property(
+            "points",
+            [](morphio::Soma* soma) {
+                return py::array(static_cast<py::ssize_t>(soma->points().size()),
+                                 soma->points().data());
+            },
+            [](morphio::Soma* soma, py::array_t<morphio::floatType> points) {
+                soma->points() = array_to_points(points);
+            },
+            "Returns the coordinates (x,y,z) of all soma point")
+
+        .def_property(
+            "diameters",
+            [](morphio::Soma* soma) {
+                return py::array(static_cast<py::ssize_t>(soma->diameters().size()),
+                                 soma->diameters().data());
+            },
+            [](morphio::Soma* soma, py::array_t<morphio::floatType> _diameters) {
+                soma->diameters() = _diameters.cast<std::vector<morphio::floatType>>();
+            },
+            "Returns the diameters of all soma points")
+
+        .def_property_readonly(
+            "center",
+            [](morphio::Soma* soma) { return py::array(3, soma->center().data()); },
+            "Returns the center of gravity of the soma points")
+
+        .def_property_readonly("max_distance",
+                               &morphio::Soma::maxDistance,
+                               "Return the maximum distance between the center of gravity "
+                               "and any of the soma points")
+
+        .def_property_readonly("type", &morphio::Soma::type, "Returns the soma type")
+
+        .def_property_readonly("surface",
+                               &morphio::Soma::surface,
+                               "Returns the soma surface\n\n"
+                               "Note: the soma surface computation depends on the soma type");
+}
diff --git a/binds/python/bind_soma.h b/binds/python/bind_soma.h
new file mode 100644
index 000000000..deed57360
--- /dev/null
+++ b/binds/python/bind_soma.h
@@ -0,0 +1,4 @@
+#pragma once
+#include <pybind11/pybind11.h>
+
+void bind_soma_module(pybind11::module& m);
diff --git a/binds/python/bindings_utils.h b/binds/python/bindings_utils.h
index 4d6ff2ebb..eb2b83f27 100644
--- a/binds/python/bindings_utils.h
+++ b/binds/python/bindings_utils.h
@@ -28,7 +28,6 @@ py::array_t<morphio::floatType> span_to_ndarray(const morphio::range<const T>& s
     return py::array(buffer_info);
 }
 
-
 /**
  * @brief "Casts" a Cpp sequence to a python array (no memory copies)
  *  Python capsule handles void pointers to objects and makes sure
diff --git a/binds/python/morphio.cpp b/binds/python/morphio.cpp
index 4afa5693c..e9a128cf1 100644
--- a/binds/python/morphio.cpp
+++ b/binds/python/morphio.cpp
@@ -2,6 +2,7 @@
 
 #include "bind_immutable.h"
 #include "bind_misc.h"
+#include "bind_soma.h"
 #include "bind_mutable.h"
 #include "bind_vasculature.h"
 
@@ -9,6 +10,7 @@ namespace py = pybind11;
 
 PYBIND11_MODULE(_morphio, m) {
     bind_misc(m);
+    bind_soma_module(m);
     bind_immutable_module(m);
 
     py::module mut_module = m.def_submodule("mut");
diff --git a/include/morphio/morphology.h b/include/morphio/morphology.h
index 9596bd108..c32e0e168 100644
--- a/include/morphio/morphology.h
+++ b/include/morphio/morphology.h
@@ -148,6 +148,7 @@ class Morphology
     const MorphologyVersion& version() const;
 
   protected:
+
     friend class mut::Morphology;
     Morphology(const Property::Properties& properties, unsigned int options);
 
diff --git a/include/morphio/mut/morphology.h b/include/morphio/mut/morphology.h
index b93492b3c..538ae2c65 100644
--- a/include/morphio/mut/morphology.h
+++ b/include/morphio/mut/morphology.h
@@ -10,8 +10,9 @@
 #include <morphio/exceptions.h>
 #include <morphio/mut/endoplasmic_reticulum.h>
 #include <morphio/mut/mitochondria.h>
+
 #include <morphio/mut/modifiers.h>
-#include <morphio/mut/soma.h>
+#include <morphio/soma.h>
 #include <morphio/properties.h>
 #include <morphio/types.h>
 
diff --git a/include/morphio/mut/soma.h b/include/morphio/mut/soma.h
deleted file mode 100644
index 016f0b533..000000000
--- a/include/morphio/mut/soma.h
+++ /dev/null
@@ -1,70 +0,0 @@
-#pragma once
-
-#include <morphio/properties.h>  // Property
-#include <morphio/soma.h>        // morphio::soma
-
-namespace morphio {
-namespace mut {
-/** Mutable(editable) morphio::Soma */
-class Soma
-{
-  public:
-    Soma() = default;
-    Soma(const Soma& soma) = default;
-
-    explicit Soma(const Property::PointLevel& pointProperties);
-    explicit Soma(const morphio::Soma& soma);
-
-    /// Return the coordinates (x,y,z) of all soma points
-    std::vector<Point>& points() noexcept {
-        return point_properties_._points;
-    }
-    const std::vector<Point>& points() const noexcept {
-        return point_properties_._points;
-    }
-
-    /// Return the diameters of all soma points
-    std::vector<morphio::floatType>& diameters() noexcept {
-        return point_properties_._diameters;
-    }
-    const std::vector<morphio::floatType>& diameters() const noexcept {
-        return point_properties_._diameters;
-    }
-
-    /// Return the soma type
-    SomaType type() const noexcept {
-        return soma_type_;
-    }
-    /**
-     * Return the center of gravity of the soma points
-     **/
-    Point center() const;
-
-    /**
-       Return the soma surface
-       Note: the soma surface computation depends on the soma type
-    **/
-    floatType surface() const;
-
-    /**
-     * Return the maximum distance between the center of gravity and any of
-     * the soma points
-     */
-    floatType maxDistance() const;
-
-    Property::PointLevel& properties() noexcept {
-        return point_properties_;
-    }
-    const Property::PointLevel& properties() const noexcept {
-        return point_properties_;
-    }
-
-  private:
-    friend class Morphology;
-
-    SomaType soma_type_ = SOMA_UNDEFINED;
-    Property::PointLevel point_properties_;
-};
-
-}  // namespace mut
-}  // namespace morphio
diff --git a/include/morphio/soma.h b/include/morphio/soma.h
index 915521567..c00f19f37 100644
--- a/include/morphio/soma.h
+++ b/include/morphio/soma.h
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <morphio/properties.h>
 #include <morphio/morphology.h>
 #include <morphio/types.h>
 
@@ -25,22 +26,33 @@ namespace morphio {
  *
  * @version unstable
  */
-class Soma
+struct Soma
 {
-  public:
+    Soma() = default;
+    Soma(const Soma& soma) = default;
+
+    explicit Soma(const Property::Properties& properties);
+    explicit Soma(const Property::PointLevel& point_properties);
+
     /// Return the  coordinates (x,y,z) of all soma points
-    range<const Point> points() const noexcept {
-        return properties_->_somaLevel._points;
+    std::vector<Point>& points() noexcept {
+        return properties_._points;
+    }
+    const std::vector<Point>& points() const noexcept {
+        return properties_._points;
     }
 
     /// Return the diameters of all soma points
-    range<const floatType> diameters() const noexcept {
-        return properties_->_somaLevel._diameters;
+    std::vector<morphio::floatType>& diameters() noexcept {
+        return properties_._diameters;
+    }
+    const std::vector<morphio::floatType>& diameters() const noexcept {
+        return properties_._diameters;
     }
 
     /// Return the soma type
     SomaType type() const noexcept {
-        return properties_->_cellLevel._somaType;
+        return soma_type_;
     }
 
     /// Return the center of gravity of the soma points
@@ -64,16 +76,18 @@ class Soma
      */
     floatType maxDistance() const;
 
+    /// Return soma properties
+    Property::PointLevel& properties() noexcept {
+        return properties_;
+    }
+    const Property::PointLevel& properties() const noexcept {
+        return properties_;
+    }
+
   private:
-    explicit Soma(const std::shared_ptr<Property::Properties>& properties);
-    // TODO: find out why the following line does not work
-    // when friend class Morphology; is removed
-    // template <typename Property>
-    // friend const morphio::Soma morphio::Morphology::soma() const;
-    friend class Morphology;
-    friend class mut::Soma;
 
-    std::shared_ptr<Property::Properties> properties_;
+    SomaType soma_type_ = SOMA_UNDEFINED;
+    Property::PointLevel properties_;
 };
 
 }  // namespace morphio
diff --git a/include/morphio/types.h b/include/morphio/types.h
index 05001f449..1f133ba25 100644
--- a/include/morphio/types.h
+++ b/include/morphio/types.h
@@ -19,7 +19,7 @@ class Section;
 template <class T>
 class SectionBase;
 
-class Soma;
+struct Soma;
 
 namespace Property {
 struct Properties;
@@ -45,7 +45,6 @@ class MitoSection;
 class Mitochondria;
 class Morphology;
 class Section;
-class Soma;
 }  // namespace mut
 
 using SectionRange = std::pair<size_t, size_t>;
diff --git a/morphio/mut/__init__.py b/morphio/mut/__init__.py
index cf0f7730a..67b73858e 100644
--- a/morphio/mut/__init__.py
+++ b/morphio/mut/__init__.py
@@ -1,6 +1,6 @@
+from .._morphio import Soma
 from .._morphio.mut import (Morphology,
                             Section,
-                            Soma,
                             MitoSection,
                             Mitochondria,
                             GlialCell,
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 8576bf922..b2b9c52b7 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -16,7 +16,6 @@ set(MORPHIO_SOURCES
     mut/modifiers.cpp
     mut/morphology.cpp
     mut/section.cpp
-    mut/soma.cpp
     mut/writers.cpp
     point_utils.cpp
     properties.cpp
diff --git a/src/morphology.cpp b/src/morphology.cpp
index f7109c77d..b820702ea 100644
--- a/src/morphology.cpp
+++ b/src/morphology.cpp
@@ -8,7 +8,6 @@
 #include <morphio/morphology.h>
 #include <morphio/section.h>
 #include <morphio/soma.h>
-
 #include <morphio/mut/morphology.h>
 
 #include "readers/morphologyASC.h"
@@ -151,7 +150,7 @@ Morphology::Morphology(const std::string& contents,
     : Morphology(loadString(contents, extension, options), options) {}
 
 Soma Morphology::soma() const {
-    return Soma(properties_);
+    return Soma(*properties_);
 }
 
 Mitochondria Morphology::mitochondria() const {
diff --git a/src/mut/morphology.cpp b/src/mut/morphology.cpp
index c1a8f2e5e..4c2e2f0e6 100644
--- a/src/mut/morphology.cpp
+++ b/src/mut/morphology.cpp
@@ -275,7 +275,7 @@ Property::Properties Morphology::buildReadOnly() const {
 
     properties._cellLevel = *_cellProperties;
     properties._cellLevel._somaType = _soma->type();
-    _appendProperties(properties._somaLevel, _soma->point_properties_);
+    _appendProperties(properties._somaLevel, _soma->properties());
 
     for (auto it = depth_begin(); it != depth_end(); ++it) {
         const std::shared_ptr<Section>& section = *it;
diff --git a/src/mut/soma.cpp b/src/mut/soma.cpp
deleted file mode 100644
index 7dbd14446..000000000
--- a/src/mut/soma.cpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#include <morphio/mut/soma.h>
-#include <morphio/soma.h>
-
-#include "../point_utils.h"     // centerOfGRavity, maxDistanceToCenterOfGravity
-#include "../shared_utils.hpp"  // _somaSurface
-
-namespace morphio {
-namespace mut {
-Soma::Soma(const Property::PointLevel& point_properties)
-    : point_properties_(point_properties) {}
-
-Soma::Soma(const morphio::Soma& soma)
-    : soma_type_(soma.type())
-    , point_properties_(soma.properties_->_somaLevel) {}
-
-Point Soma::center() const {
-    return centerOfGravity(points());
-}
-
-floatType Soma::surface() const {
-    return _somaSurface(type(), diameters(), points());
-}
-
-floatType Soma::maxDistance() const {
-    return maxDistanceToCenterOfGravity(point_properties_._points);
-}
-
-}  // end namespace mut
-}  // end namespace morphio
diff --git a/src/mut/writers.cpp b/src/mut/writers.cpp
index 7fd34d4b5..cbc7986af 100644
--- a/src/mut/writers.cpp
+++ b/src/mut/writers.cpp
@@ -62,7 +62,7 @@ bool _skipDuplicate(const std::shared_ptr<morphio::mut::Section>& section) {
     return section->diameters().front() == section->parent()->diameters().back();
 }
 
-void checkSomaHasSameNumberPointsDiameters(const morphio::mut::Soma& soma) {
+void checkSomaHasSameNumberPointsDiameters(const morphio::Soma& soma) {
     const size_t n_points = soma.points().size();
     const size_t n_diameters = soma.diameters().size();
 
diff --git a/src/readers/morphologySWC.cpp b/src/readers/morphologySWC.cpp
index 5db84a7ee..9753d7f9f 100644
--- a/src/readers/morphologySWC.cpp
+++ b/src/readers/morphologySWC.cpp
@@ -10,7 +10,6 @@
 #include <morphio/errorMessages.h>
 #include <morphio/mut/morphology.h>
 #include <morphio/mut/section.h>
-#include <morphio/mut/soma.h>
 #include <morphio/properties.h>
 
 namespace {
diff --git a/src/shared_utils.hpp b/src/shared_utils.hpp
index 7cf779a73..9289bdad9 100644
--- a/src/shared_utils.hpp
+++ b/src/shared_utils.hpp
@@ -18,41 +18,6 @@ std::string valueToString(const T a) {
     return std::to_string(a);
 }
 
-inline floatType _somaSurface(const SomaType type,
-                              const range<const floatType>& diameters,
-                              const range<const Point>& points) {
-    size_t size = points.size();
-
-    switch (type) {
-    case SOMA_SINGLE_POINT:
-    case SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS: {
-        floatType radius = diameters[0] / 2;
-        return 4 * morphio::PI * radius * radius;
-    }
-    case SOMA_CYLINDERS: {
-        // Surface is approximated as the sum of areas of the conical frustums
-        // defined by each segments. Does not include the endcaps areas
-        floatType surface = 0;
-        for (unsigned int i = 0; i < size - 1; ++i) {
-            floatType r0 = static_cast<morphio::floatType>(diameters[i]) * floatType{0.5};
-            floatType r1 = static_cast<morphio::floatType>(diameters[i + 1]) * floatType{0.5};
-            floatType h2 = euclidean_distance(points[i], points[i + 1]);
-            auto s = morphio::PI * (r0 + r1) * std::sqrt((r0 - r1) * (r0 - r1) + h2 * h2);
-            surface += s;
-        }
-        return surface;
-    }
-    case SOMA_SIMPLE_CONTOUR: {
-        throw NotImplementedError("Surface is not implemented for SOMA_SIMPLE_CONTOUR");
-    }
-    case SOMA_UNDEFINED:
-    default: {
-        morphio::readers::ErrorMessages err;
-        throw SomaError(err.ERROR_NOT_IMPLEMENTED_UNDEFINED_SOMA("Soma::surface"));
-    }
-    }
-}
-
 template <typename T>
 void _appendVector(std::vector<T>& to, const std::vector<T>& from, int offset) {
     to.insert(to.end(), from.begin() + offset, from.end());
diff --git a/src/soma.cpp b/src/soma.cpp
index e86b9b7aa..2d37037b9 100644
--- a/src/soma.cpp
+++ b/src/soma.cpp
@@ -1,20 +1,31 @@
+#include <limits>
+
+#include <morphio/errorMessages.h>
 #include <morphio/section.h>
 #include <morphio/soma.h>
 #include <morphio/vector_types.h>
 
-#include "point_utils.h"
-#include "shared_utils.hpp"
+#include "point_utils.h"  // centerOfGravity, euclideanDistance, maxDistanceToCenterOfGravity
 
 namespace morphio {
-Soma::Soma(const std::shared_ptr<Property::Properties>& properties)
-    : properties_(properties) {}
+
+Soma::Soma(const Property::Properties& properties)
+    : soma_type_(properties._cellLevel._somaType)
+    , properties_(properties._somaLevel) {}
+
+
+Soma::Soma(const Property::PointLevel& point_properties)
+    : properties_(point_properties) {}
+
 
 Point Soma::center() const {
-    return centerOfGravity(properties_->_somaLevel._points);
+    return centerOfGravity(points());
 }
 
+
 floatType Soma::volume() const {
-    switch (properties_->_cellLevel._somaType) {
+
+    switch (soma_type_) {
     case SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS: {
         floatType radius = diameters()[0] / 2;
         return 4 * morphio::PI * radius * radius;
@@ -25,16 +36,50 @@ floatType Soma::volume() const {
     case SOMA_SIMPLE_CONTOUR:
     case SOMA_UNDEFINED:
     default:
-        throw;
+        morphio::readers::ErrorMessages err;
+        throw SomaError(err.ERROR_NOT_IMPLEMENTED_UNDEFINED_SOMA("Soma::volume"));
     }
 }
 
+
 floatType Soma::surface() const {
-    return _somaSurface(type(), diameters(), points());
+    const auto& soma_points = points();
+    const auto& soma_diameters = diameters();
+
+    size_t size = soma_points.size();
+
+    switch (soma_type_) {
+    case SOMA_SINGLE_POINT:
+    case SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS: {
+        floatType radius = soma_diameters[0] / 2;
+        return 4 * morphio::PI * radius * radius;
+    }
+    case SOMA_CYLINDERS: {
+        // Surface is approximated as the sum of areas of the conical frustums
+        // defined by each segments. Does not include the endcaps areas
+        floatType surface = 0;
+        for (unsigned int i = 0; i < size - 1; ++i) {
+            floatType r0 = static_cast<morphio::floatType>(soma_diameters[i]) * floatType{0.5};
+            floatType r1 = static_cast<morphio::floatType>(soma_diameters[i + 1]) * floatType{0.5};
+            floatType h2 = euclidean_distance(soma_points[i], soma_points[i + 1]);
+            auto s = morphio::PI * (r0 + r1) * std::sqrt((r0 - r1) * (r0 - r1) + h2 * h2);
+            surface += s;
+        }
+        return surface;
+    }
+    case SOMA_SIMPLE_CONTOUR: {
+        throw NotImplementedError("Surface is not implemented for SOMA_SIMPLE_CONTOUR");
+    }
+    case SOMA_UNDEFINED:
+    default: {
+        morphio::readers::ErrorMessages err;
+        throw SomaError(err.ERROR_NOT_IMPLEMENTED_UNDEFINED_SOMA("Soma::surface"));
+    }
+    }
 }
 
+
 floatType Soma::maxDistance() const {
-    return maxDistanceToCenterOfGravity(properties_->_somaLevel._points);
+    return maxDistanceToCenterOfGravity(points());
 }
-
 }  // namespace morphio
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index df3ed349b..076263936 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,5 +1,6 @@
 set(TESTS_SRC
         main.cpp
+        test_soma.cpp
         test_immutable_morphology.cpp
         test_mitochondria.cpp
         test_morphology_readers.cpp
diff --git a/tests/test_mutable_morphology.cpp b/tests/test_mutable_morphology.cpp
index 9384155a7..d86a5fc74 100644
--- a/tests/test_mutable_morphology.cpp
+++ b/tests/test_mutable_morphology.cpp
@@ -6,6 +6,17 @@
 #include <filesystem>
 namespace fs = std::filesystem;
 
+TEST_CASE("soma", "[mutableMorphology]") {
+
+    auto morph = morphio::mut::Morphology("data/simple-heterogeneous-neurite.swc");
+
+    morph.soma();
+    morph.soma()->points();
+    morph.soma()->diameters();
+
+    morph.soma()->diameters() = {2., 3., 4.};
+}
+
 TEST_CASE("isHeterogeneous", "[mutableMorphology]") {
     auto morph = morphio::mut::Morphology("data/simple-heterogeneous-neurite.swc");
 
diff --git a/tests/test_soma.cpp b/tests/test_soma.cpp
new file mode 100644
index 000000000..a083be0ce
--- /dev/null
+++ b/tests/test_soma.cpp
@@ -0,0 +1,142 @@
+
+#include <catch2/catch.hpp>
+#include <cmath>
+
+#include <morphio/enums.h>
+#include <morphio/soma.h>
+
+#include <morphio/errorMessages.h>
+#include <morphio/morphology.h>
+#include <morphio/mut/morphology.h>
+
+
+TEST_CASE("soma-default-constructor", "[soma]") {
+    const morphio::Soma soma;
+
+    REQUIRE(soma.type() == morphio::enums::SomaType::SOMA_UNDEFINED);
+    REQUIRE(soma.points().empty());
+    REQUIRE(soma.diameters().empty());
+
+    REQUIRE(soma.properties()._points.empty());
+    REQUIRE(soma.properties()._diameters.empty());
+    REQUIRE(soma.properties()._perimeters.empty());
+
+    morphio::Point center = soma.center();
+
+    REQUIRE(std::isnan(center[0]));
+    REQUIRE(std::isnan(center[1]));
+    REQUIRE(std::isnan(center[2]));
+
+    REQUIRE_THROWS_AS(soma.surface(), morphio::SomaError);
+    REQUIRE_THROWS_AS(soma.volume(), morphio::SomaError);
+
+    REQUIRE(soma.maxDistance() == 0);
+}
+
+
+TEST_CASE("soma-point-properties-constructor", "[soma]") {
+    auto properties = morphio::Property::PointLevel();
+
+    properties._points = {{0., 1., 2.}, {3., 4., 5.}};
+    properties._diameters = {0.2, 0.3};
+
+    morphio::Soma soma(properties);
+
+    REQUIRE(soma.points() == properties._points);
+    REQUIRE(soma.diameters() == properties._diameters);
+    REQUIRE(soma.type() == morphio::enums::SomaType::SOMA_UNDEFINED);
+
+    REQUIRE(soma.center() == morphio::Point({1.5, 2.5, 3.5}));
+    REQUIRE_THAT(soma.maxDistance(), Catch::WithinAbs(2.598076, 1e-6));
+}
+
+
+TEST_CASE("soma-properties-constructor", "[soma]") {
+    auto properties = morphio::Property::Properties();
+
+    std::vector<morphio::Point> expected_points = {{0., 1., 2.}, {3., 4., 5.}, {6., 7., 8.}};
+    std::vector<morphio::floatType> expected_diameters = {0.2, 0.3, 0.4};
+    morphio::enums::SomaType expected_soma_type = morphio::enums::SomaType::SOMA_SIMPLE_CONTOUR;
+
+    properties._somaLevel._points = expected_points;
+    properties._somaLevel._diameters = expected_diameters;
+    properties._cellLevel._somaType = expected_soma_type;
+
+    morphio::Soma soma(properties);
+
+    REQUIRE(soma.points() == expected_points);
+    REQUIRE(soma.diameters() == expected_diameters);
+    REQUIRE(soma.type() == expected_soma_type);
+}
+
+
+TEST_CASE("soma-copy-constructor", "[soma]") {
+    auto properties = morphio::Property::PointLevel();
+
+    properties._points = {{0., 1., 2.}, {3., 4., 5.}};
+    properties._diameters = {0.2, 0.3};
+
+    // allocate a soma in the heap to delete it afterwards
+    // and check if the copy maintains its data
+    morphio::Soma* soma = new morphio::Soma(properties);
+    morphio::Soma soma_copy(*soma);
+
+    auto expected_points = soma->points();
+    auto expected_diameters = soma->diameters();
+    auto expected_type = soma->type();
+
+    REQUIRE(expected_points == soma_copy.points());
+    REQUIRE(expected_diameters == soma_copy.diameters());
+    REQUIRE(expected_type == soma_copy.type());
+
+    delete soma;
+
+    REQUIRE(expected_points == soma_copy.points());
+    REQUIRE(expected_diameters == soma_copy.diameters());
+    REQUIRE(expected_type == soma_copy.type());
+}
+
+
+TEST_CASE("soma-immutable-morphology-constructor", "[soma]") {
+    const auto soma = morphio::Morphology("data/h5/v1/Neuron.h5").soma();
+
+    REQUIRE(soma.type() == morphio::SomaType::SOMA_SIMPLE_CONTOUR);
+
+    std::vector<morphio::Point> expected_soma_points = {{0.0, 0.0, 0.0},
+                                                        {0.0, 0.2, 0.0},
+                                                        {0.1, 0.1, 0.0}};
+    REQUIRE(soma.points() == expected_soma_points);
+
+    std::vector<morphio::floatType> expected_soma_diameters = {0.2, 0.2, 0.2};
+    REQUIRE(soma.diameters() == expected_soma_diameters);
+}
+
+
+TEST_CASE("soma-mutable-morphology-constructor", "[soma]") {
+    auto morph = morphio::mut::Morphology("data/h5/v1/Neuron.h5");
+
+    const auto soma = *morph.soma();
+
+    REQUIRE(soma.type() == morphio::SomaType::SOMA_SIMPLE_CONTOUR);
+
+    std::vector<morphio::Point> expected_soma_points = {{0.0, 0.0, 0.0},
+                                                        {0.0, 0.2, 0.0},
+                                                        {0.1, 0.1, 0.0}};
+    REQUIRE(soma.points() == expected_soma_points);
+
+    std::vector<morphio::floatType> expected_soma_diameters = {0.2, 0.2, 0.2};
+    REQUIRE(soma.diameters() == expected_soma_diameters);
+}
+
+
+TEST_CASE("soma-mutable-morphology-mutation", "[soma]") {
+    auto morph = morphio::mut::Morphology("data/h5/v1/Neuron.h5");
+
+    std::vector<morphio::Point> expected_soma_points = {{0.1, 0.1, 0.0}};
+    morph.soma()->points() = expected_soma_points;
+    REQUIRE(morph.soma()->points() == expected_soma_points);
+
+    std::vector<morphio::floatType> expected_soma_diameters = {3.0};
+    morph.soma()->diameters() = expected_soma_diameters;
+    REQUIRE(morph.soma()->diameters() == expected_soma_diameters);
+}