diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index eaef22f99..685609d4c 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -5,6 +5,9 @@ set(EXAMPLE_FILES
     "forwardDynamicsExample.cpp"
     "inverseDynamicsExample.cpp"
 )
+if (NOT BIORBD_USE_CASADI_MATH)
+    list(APPEND EXAMPLE_FILES "writeModelExample.cpp")
+endif()
 if (MODULE_MUSCLES)
     list(APPEND EXAMPLE_FILES "forwardDynamicsFromMusclesExample.cpp")
 endif()
diff --git a/examples/python3/write_model.py b/examples/python3/write_model.py
new file mode 100644
index 000000000..ac418a695
--- /dev/null
+++ b/examples/python3/write_model.py
@@ -0,0 +1,37 @@
+"""
+This examples shows how to
+    1. Load a model
+    2. Modify one of its properties (here, the mass of a segment)
+    3. Write the modified model to a new .bioMod file
+    4. Reload the written file to confirm the change was saved
+"""
+
+from pathlib import Path
+
+import biorbd
+
+
+def main():
+    # Load a predefined model
+    current_file_dir = Path(__file__).parent
+    model = biorbd.Model(f"{current_file_dir}/../pyomecaman.bioMod")
+
+    # Modify the mass of the first segment
+    segment = model.segment(0)
+    print(f"Original mass of {segment.name().to_string()}: {segment.characteristics().mass()}")
+    segment.characteristics().setMass(segment.characteristics().mass() * 2)
+
+    # Write the modified model to a new file
+    output_path = str(current_file_dir / "write_model_output.bioMod")
+    biorbd.Writer.writeModel(model, output_path)
+
+    # Reload the file to confirm the modification was properly saved
+    reloaded_model = biorbd.Model(output_path)
+    print(
+        f"Mass of {segment.name().to_string()} after write/reload: "
+        f"{reloaded_model.segment(0).characteristics().mass()}"
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/writeModelExample.cpp b/examples/writeModelExample.cpp
new file mode 100644
index 000000000..c665493f0
--- /dev/null
+++ b/examples/writeModelExample.cpp
@@ -0,0 +1,36 @@
+#include "biorbd.h"
+
+///
+/// \brief main Modify a model and write it back to a new .bioMod file
+/// \return Nothing
+///
+/// This examples shows how to
+///     1. Load a model
+///     2. Modify one of its properties (here, the mass of a segment)
+///     3. Write the modified model to a new .bioMod file
+///     4. Reload the written file to confirm the change was saved
+///
+
+using namespace BIORBD_NAMESPACE;
+
+int main() {
+  // Load a predefined model
+  Model model("pyomecaman.bioMod");
+
+  // Modify the mass of the first segment
+  rigidbody::Segment& segment = model.segment(0);
+  std::cout << "Original mass of " << segment.name() << ": "
+            << segment.characteristics().mass() << std::endl;
+  segment.characteristics().setMass(segment.characteristics().mass() * 2);
+
+  // Write the modified model to a new file
+  utils::Path outputPath("writeModelExample_output.bioMod");
+  Writer::writeModel(model, outputPath);
+
+  // Reload the file to confirm the modification was properly saved
+  Model reloadedModel(outputPath.relativePath());
+  std::cout << "Mass of " << segment.name() << " after write/reload: "
+            << reloadedModel.segment(0).characteristics().mass() << std::endl;
+
+  return 0;
+}
diff --git a/src/ModelWriter.cpp b/src/ModelWriter.cpp
index 7498b94b3..47ec11ded 100644
--- a/src/ModelWriter.cpp
+++ b/src/ModelWriter.cpp
@@ -12,6 +12,7 @@
 #include "RigidBody/SegmentCharacteristics.h"
 #include "Utils/Matrix3d.h"
 #include "Utils/Path.h"
+#include "Utils/Range.h"
 #include "Utils/String.h"
 #include "Utils/Vector.h"
 
@@ -64,6 +65,32 @@ void Writer::writeModel(Model& model, const utils::Path& pathToWrite) {
       biorbdModelFile << sep << sep << "rotations" << sep
                       << model.segment(i).seqR() << std::endl;
     }
+    const std::vector<utils::Range>& qRanges = model.segment(i).QRanges();
+    if (qRanges.size() > 0) {
+      biorbdModelFile << sep << sep << "rangesQ" << std::endl;
+      for (size_t j = 0; j < qRanges.size(); ++j) {
+        biorbdModelFile << sep << sep << sep << qRanges[j].min() << sep
+                        << qRanges[j].max() << std::endl;
+      }
+    }
+    const std::vector<utils::Range>& qdotRanges =
+        model.segment(i).QdotRanges();
+    if (qdotRanges.size() > 0) {
+      biorbdModelFile << sep << sep << "rangesQdot" << std::endl;
+      for (size_t j = 0; j < qdotRanges.size(); ++j) {
+        biorbdModelFile << sep << sep << sep << qdotRanges[j].min() << sep
+                        << qdotRanges[j].max() << std::endl;
+      }
+    }
+    const std::vector<utils::Range>& qddotRanges =
+        model.segment(i).QddotRanges();
+    if (qddotRanges.size() > 0) {
+      biorbdModelFile << sep << sep << "rangesQddot" << std::endl;
+      for (size_t j = 0; j < qddotRanges.size(); ++j) {
+        biorbdModelFile << sep << sep << sep << qddotRanges[j].min() << sep
+                        << qddotRanges[j].max() << std::endl;
+      }
+    }
     biorbdModelFile << sep << sep << "jointDampings" << sep;
     for (auto damping : model.segment(i).jointDampings()) {
       biorbdModelFile << sep << damping;
diff --git a/test/models/two_segments.bioMod b/test/models/two_segments.bioMod
index bf6f7bcb9..a9789dcdd 100644
--- a/test/models/two_segments.bioMod
+++ b/test/models/two_segments.bioMod
@@ -4,6 +4,18 @@ segment segment1
     rt 0.1 0.2 0.3 xyz 1 2 3
     translations    xy
     rotations	z
+    rangesQ
+        -1.5 2.5
+        -1.6 2.6
+        -1.7 2.7
+    rangesQdot
+        -20.1 25.1
+        -20.2 25.2
+        -20.3 25.3
+    rangesQddot
+        -300.1 350.1
+        -300.2 350.2
+        -300.3 350.3
 	jointdampings 1 2 3
 endsegment
 
@@ -11,6 +23,12 @@ segment segment2
     parent segment1
     rt 0.1 0.2 0.3 xyz 1 2 3
     rotations	z
+    rangesQ
+        -0.5 0.6
+    rangesQdot
+        -12.5 13.5
+    rangesQddot
+        -140.5 150.5
 endsegment
 
     marker mark_seg2
diff --git a/test/test_biorbd.cpp b/test/test_biorbd.cpp
index 2e7525457..71035e816 100644
--- a/test/test_biorbd.cpp
+++ b/test/test_biorbd.cpp
@@ -16,6 +16,7 @@
 #include "RigidBody/MeshFace.h"
 #include "RigidBody/NodeSegment.h"
 #include "RigidBody/Segment.h"
+#include "Utils/Range.h"
 #include "Utils/RotoTrans.h"
 #include "Utils/RotoTransNode.h"
 #include "Utils/String.h"
@@ -60,6 +61,37 @@ TEST(FileIO, WriteModel) {
   EXPECT_FLOAT_EQ(model.segment(0).jointDampings()[2], 3.0);
   EXPECT_FLOAT_EQ(model.segment(1).jointDampings()[0], 0.0);
 
+  for (size_t k = 0; k < model.nbSegment(); ++k) {
+    const std::vector<utils::Range>& QRanges = model.segment(k).QRanges();
+    const std::vector<utils::Range>& QRangesCopy =
+        modelCopy.segment(k).QRanges();
+    ASSERT_EQ(QRangesCopy.size(), QRanges.size());
+    for (size_t i = 0; i < QRanges.size(); ++i) {
+      EXPECT_FLOAT_EQ(QRangesCopy[i].min(), QRanges[i].min());
+      EXPECT_FLOAT_EQ(QRangesCopy[i].max(), QRanges[i].max());
+    }
+
+    const std::vector<utils::Range>& QdotRanges =
+        model.segment(k).QdotRanges();
+    const std::vector<utils::Range>& QdotRangesCopy =
+        modelCopy.segment(k).QdotRanges();
+    ASSERT_EQ(QdotRangesCopy.size(), QdotRanges.size());
+    for (size_t i = 0; i < QdotRanges.size(); ++i) {
+      EXPECT_FLOAT_EQ(QdotRangesCopy[i].min(), QdotRanges[i].min());
+      EXPECT_FLOAT_EQ(QdotRangesCopy[i].max(), QdotRanges[i].max());
+    }
+
+    const std::vector<utils::Range>& QddotRanges =
+        model.segment(k).QddotRanges();
+    const std::vector<utils::Range>& QddotRangesCopy =
+        modelCopy.segment(k).QddotRanges();
+    ASSERT_EQ(QddotRangesCopy.size(), QddotRanges.size());
+    for (size_t i = 0; i < QddotRanges.size(); ++i) {
+      EXPECT_FLOAT_EQ(QddotRangesCopy[i].min(), QddotRanges[i].min());
+      EXPECT_FLOAT_EQ(QddotRangesCopy[i].max(), QddotRanges[i].max());
+    }
+  }
+
   for (size_t k = 0; k < model.nbSegment(); ++k) {
     for (size_t i = 0; i < 4; ++i) {
       for (size_t j = 0; j < 4; ++j) {