diff --git a/models/model_calvetti_py/dwdw.cpp b/models/model_calvetti_py/dwdw.cpp
index 4fa96d39f4..70f1cba23b 100644
--- a/models/model_calvetti_py/dwdw.cpp
+++ b/models/model_calvetti_py/dwdw.cpp
@@ -7,8 +7,8 @@
 namespace amici {
 namespace model_model_calvetti_py {
 
-static constexpr std::array<sunindextype, 18> dwdw_colptrs_model_calvetti_py_ = {
-    0, 1, 2, 3, 4, 5, 6, 6, 6, 7, 9, 12, 16, 16, 16, 16, 16, 16
+static constexpr std::array<sunindextype, 17> dwdw_colptrs_model_calvetti_py_ = {
+    0, 1, 2, 3, 4, 5, 6, 6, 6, 7, 9, 12, 16, 16, 16, 16, 16
 };
 
 void dwdw_colptrs_model_calvetti_py(SUNMatrixWrapper &dwdw){
diff --git a/models/model_calvetti_py/dxdotdw.cpp b/models/model_calvetti_py/dxdotdw.cpp
index ac66441f99..07c21529b0 100644
--- a/models/model_calvetti_py/dxdotdw.cpp
+++ b/models/model_calvetti_py/dxdotdw.cpp
@@ -7,8 +7,8 @@
 namespace amici {
 namespace model_model_calvetti_py {
 
-static constexpr std::array<sunindextype, 18> dxdotdw_colptrs_model_calvetti_py_ = {
-    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 5, 7, 7
+static constexpr std::array<sunindextype, 17> dxdotdw_colptrs_model_calvetti_py_ = {
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 5, 7
 };
 
 void dxdotdw_colptrs_model_calvetti_py(SUNMatrixWrapper &dxdotdw){
diff --git a/models/model_calvetti_py/model_calvetti_py.cpp b/models/model_calvetti_py/model_calvetti_py.cpp
index 79758b2a17..7ddbedab18 100644
--- a/models/model_calvetti_py/model_calvetti_py.cpp
+++ b/models/model_calvetti_py/model_calvetti_py.cpp
@@ -45,7 +45,7 @@ ObservableScaling::lin, // y[4]
 ObservableScaling::lin, // y[5]
 };
 
-std::array<const char*, 17> expression_names = {
+std::array<const char*, 16> expression_names = {
     "C1ss", // w[0]
 "C2ss", // w[1]
 "C3ss", // w[2]
@@ -62,7 +62,6 @@ std::array<const char*, 17> expression_names = {
 "rate_of_V1", // w[13]
 "rate_of_V2", // w[14]
 "rate_of_V3", // w[15]
-"flux_r0", // w[16]
 };
 
 std::array<const char*, 0> free_parameter_ids = {
@@ -96,7 +95,7 @@ std::array<const char*, 6> observable_ids = {
 "obs_f2", // y[5]
 };
 
-std::array<const char*, 17> expression_ids = {
+std::array<const char*, 16> expression_ids = {
     "C1ss", // w[0]
 "C2ss", // w[1]
 "C3ss", // w[2]
@@ -113,7 +112,6 @@ std::array<const char*, 17> expression_ids = {
 "rate_of_V1", // w[13]
 "rate_of_V2", // w[14]
 "rate_of_V3", // w[15]
-"flux_r0", // w[16]
 };
 
 std::array<int, 6> state_idxs_solver = {
diff --git a/models/model_calvetti_py/model_calvetti_py.h b/models/model_calvetti_py/model_calvetti_py.h
index 70943feeb3..fcb5b0451c 100644
--- a/models/model_calvetti_py/model_calvetti_py.h
+++ b/models/model_calvetti_py/model_calvetti_py.h
@@ -19,12 +19,12 @@ extern std::array<const char*, 6> fixed_parameter_names;
 extern std::array<const char*, 6> state_names;
 extern std::array<const char*, 6> observable_names;
 extern std::array<const ObservableScaling, 6> observable_scalings;
-extern std::array<const char*, 17> expression_names;
+extern std::array<const char*, 16> expression_names;
 extern std::array<const char*, 0> free_parameter_ids;
 extern std::array<const char*, 6> fixed_parameter_ids;
 extern std::array<const char*, 6> state_ids;
 extern std::array<const char*, 6> observable_ids;
-extern std::array<const char*, 17> expression_ids;
+extern std::array<const char*, 16> expression_ids;
 extern std::array<int, 6> state_idxs_solver;
 
 extern void Jy_model_calvetti_py(realtype *Jy, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my);
@@ -125,7 +125,7 @@ class Model_model_calvetti_py : public amici::Model_DAE {
                   .ne = 4,
                   .ne_solver = 0,
                   .nspl = 0,
-                  .nw = 17,
+                  .nw = 16,
                   .ndwdx = 15,
                   .ndwdp = 0,
                   .ndwdw = 16,
@@ -557,7 +557,7 @@ class Model_model_calvetti_py : public amici::Model_DAE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "63f2e78a186ee722a79b0da2797e0c116ec81065";
+        return "c4a395b6462dc1f26ad1472e6b2fa71fe0d54382";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_calvetti_py/w.h b/models/model_calvetti_py/w.h
index 796ffea41b..1da13469c4 100644
--- a/models/model_calvetti_py/w.h
+++ b/models/model_calvetti_py/w.h
@@ -14,4 +14,3 @@
 #define rate_of_V1 w[13]
 #define rate_of_V2 w[14]
 #define rate_of_V3 w[15]
-#define flux_r0 w[16]
diff --git a/models/model_dirac_py/CMakeLists.txt b/models/model_dirac_py/CMakeLists.txt
index 1273db6cb8..dadf82ddf1 100644
--- a/models/model_dirac_py/CMakeLists.txt
+++ b/models/model_dirac_py/CMakeLists.txt
@@ -88,7 +88,6 @@ sigmay.h
 stau.cpp 
 stau.h 
 sx.h 
-w.h 
 wrapfunctions.cpp 
 wrapfunctions.h 
 x.h 
diff --git a/models/model_dirac_py/deltaqB.cpp b/models/model_dirac_py/deltaqB.cpp
index b95bd81ee1..3319c79f39 100644
--- a/models/model_dirac_py/deltaqB.cpp
+++ b/models/model_dirac_py/deltaqB.cpp
@@ -5,7 +5,6 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
diff --git a/models/model_dirac_py/dxdotdp_explicit.cpp b/models/model_dirac_py/dxdotdp_explicit.cpp
index 4211e8beda..300acbc0e6 100644
--- a/models/model_dirac_py/dxdotdp_explicit.cpp
+++ b/models/model_dirac_py/dxdotdp_explicit.cpp
@@ -48,7 +48,6 @@ void dxdotdp_explicit_rowvals_model_dirac_py(SUNMatrixWrapper &dxdotdp_explicit)
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "dxdotdp_explicit.h"
 
 namespace amici {
diff --git a/models/model_dirac_py/dxdotdx_explicit.cpp b/models/model_dirac_py/dxdotdx_explicit.cpp
index 245f539fd3..ddba8ddadf 100644
--- a/models/model_dirac_py/dxdotdx_explicit.cpp
+++ b/models/model_dirac_py/dxdotdx_explicit.cpp
@@ -48,7 +48,6 @@ void dxdotdx_explicit_rowvals_model_dirac_py(SUNMatrixWrapper &dxdotdx_explicit)
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "dxdotdx_explicit.h"
 
 namespace amici {
diff --git a/models/model_dirac_py/model_dirac_py.cpp b/models/model_dirac_py/model_dirac_py.cpp
index c173bf7726..e52ab31800 100644
--- a/models/model_dirac_py/model_dirac_py.cpp
+++ b/models/model_dirac_py/model_dirac_py.cpp
@@ -29,8 +29,8 @@ std::array<const ObservableScaling, 1> observable_scalings = {
     ObservableScaling::lin, // y[0]
 };
 
-std::array<const char*, 1> expression_names = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_names = {
+    
 };
 
 std::array<const char*, 4> free_parameter_ids = {
@@ -53,8 +53,8 @@ std::array<const char*, 1> observable_ids = {
     "obs_x2", // y[0]
 };
 
-std::array<const char*, 1> expression_ids = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_ids = {
+    
 };
 
 std::array<int, 2> state_idxs_solver = {
diff --git a/models/model_dirac_py/model_dirac_py.h b/models/model_dirac_py/model_dirac_py.h
index a03fcf9d77..093d78c2f9 100644
--- a/models/model_dirac_py/model_dirac_py.h
+++ b/models/model_dirac_py/model_dirac_py.h
@@ -19,12 +19,12 @@ extern std::array<const char*, 0> fixed_parameter_names;
 extern std::array<const char*, 2> state_names;
 extern std::array<const char*, 1> observable_names;
 extern std::array<const ObservableScaling, 1> observable_scalings;
-extern std::array<const char*, 1> expression_names;
+extern std::array<const char*, 0> expression_names;
 extern std::array<const char*, 4> free_parameter_ids;
 extern std::array<const char*, 0> fixed_parameter_ids;
 extern std::array<const char*, 2> state_ids;
 extern std::array<const char*, 1> observable_ids;
-extern std::array<const char*, 1> expression_ids;
+extern std::array<const char*, 0> expression_ids;
 extern std::array<int, 2> state_idxs_solver;
 
 extern void Jy_model_dirac_py(realtype *Jy, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my);
@@ -125,7 +125,7 @@ class Model_model_dirac_py : public amici::Model_ODE {
                   .ne = 1,
                   .ne_solver = 0,
                   .nspl = 0,
-                  .nw = 1,
+                  .nw = 0,
                   .ndwdx = 0,
                   .ndwdp = 0,
                   .ndwdw = 0,
@@ -544,7 +544,7 @@ class Model_model_dirac_py : public amici::Model_ODE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "63f2e78a186ee722a79b0da2797e0c116ec81065";
+        return "c4a395b6462dc1f26ad1472e6b2fa71fe0d54382";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_dirac_py/w.h b/models/model_dirac_py/w.h
deleted file mode 100644
index c399413a8c..0000000000
--- a/models/model_dirac_py/w.h
+++ /dev/null
@@ -1 +0,0 @@
-#define flux_r0 w[0]
diff --git a/models/model_dirac_py/xdot.cpp b/models/model_dirac_py/xdot.cpp
index e077d5eb8e..c7f6f6e46f 100644
--- a/models/model_dirac_py/xdot.cpp
+++ b/models/model_dirac_py/xdot.cpp
@@ -5,7 +5,6 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 
 namespace amici {
diff --git a/models/model_events_py/CMakeLists.txt b/models/model_events_py/CMakeLists.txt
index c7706edf3c..3c56163a8f 100644
--- a/models/model_events_py/CMakeLists.txt
+++ b/models/model_events_py/CMakeLists.txt
@@ -104,7 +104,6 @@ stau.cpp
 stau.h 
 sx.h 
 sx0_fixedParameters.cpp 
-w.h 
 wrapfunctions.cpp 
 wrapfunctions.h 
 x.h 
diff --git a/models/model_events_py/deltaqB.cpp b/models/model_events_py/deltaqB.cpp
index 4200430baf..060f96710d 100644
--- a/models/model_events_py/deltaqB.cpp
+++ b/models/model_events_py/deltaqB.cpp
@@ -6,7 +6,6 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
diff --git a/models/model_events_py/deltaxB.cpp b/models/model_events_py/deltaxB.cpp
index 0409fb5553..b377b92f8c 100644
--- a/models/model_events_py/deltaxB.cpp
+++ b/models/model_events_py/deltaxB.cpp
@@ -6,7 +6,6 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
diff --git a/models/model_events_py/dxdotdp_explicit.cpp b/models/model_events_py/dxdotdp_explicit.cpp
index e055edd533..ab288ad363 100644
--- a/models/model_events_py/dxdotdp_explicit.cpp
+++ b/models/model_events_py/dxdotdp_explicit.cpp
@@ -49,7 +49,6 @@ void dxdotdp_explicit_rowvals_model_events_py(SUNMatrixWrapper &dxdotdp_explicit
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "dxdotdp_explicit.h"
 
 namespace amici {
diff --git a/models/model_events_py/dxdotdx_explicit.cpp b/models/model_events_py/dxdotdx_explicit.cpp
index a7f128a5d0..a23dd73cd0 100644
--- a/models/model_events_py/dxdotdx_explicit.cpp
+++ b/models/model_events_py/dxdotdx_explicit.cpp
@@ -49,7 +49,6 @@ void dxdotdx_explicit_rowvals_model_events_py(SUNMatrixWrapper &dxdotdx_explicit
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "dxdotdx_explicit.h"
 
 namespace amici {
diff --git a/models/model_events_py/model_events_py.cpp b/models/model_events_py/model_events_py.cpp
index 056e4bf79c..2f5da7d8a7 100644
--- a/models/model_events_py/model_events_py.cpp
+++ b/models/model_events_py/model_events_py.cpp
@@ -33,8 +33,8 @@ std::array<const ObservableScaling, 1> observable_scalings = {
     ObservableScaling::lin, // y[0]
 };
 
-std::array<const char*, 1> expression_names = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_names = {
+    
 };
 
 std::array<const char*, 4> free_parameter_ids = {
@@ -61,8 +61,8 @@ std::array<const char*, 1> observable_ids = {
     "y1", // y[0]
 };
 
-std::array<const char*, 1> expression_ids = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_ids = {
+    
 };
 
 std::array<int, 3> state_idxs_solver = {
diff --git a/models/model_events_py/model_events_py.h b/models/model_events_py/model_events_py.h
index 16fde1877b..abc18ac9e5 100644
--- a/models/model_events_py/model_events_py.h
+++ b/models/model_events_py/model_events_py.h
@@ -19,12 +19,12 @@ extern std::array<const char*, 4> fixed_parameter_names;
 extern std::array<const char*, 3> state_names;
 extern std::array<const char*, 1> observable_names;
 extern std::array<const ObservableScaling, 1> observable_scalings;
-extern std::array<const char*, 1> expression_names;
+extern std::array<const char*, 0> expression_names;
 extern std::array<const char*, 4> free_parameter_ids;
 extern std::array<const char*, 4> fixed_parameter_ids;
 extern std::array<const char*, 3> state_ids;
 extern std::array<const char*, 1> observable_ids;
-extern std::array<const char*, 1> expression_ids;
+extern std::array<const char*, 0> expression_ids;
 extern std::array<int, 3> state_idxs_solver;
 
 extern void Jy_model_events_py(realtype *Jy, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my);
@@ -125,7 +125,7 @@ class Model_model_events_py : public amici::Model_ODE {
                   .ne = 6,
                   .ne_solver = 2,
                   .nspl = 0,
-                  .nw = 1,
+                  .nw = 0,
                   .ndwdx = 0,
                   .ndwdp = 0,
                   .ndwdw = 0,
@@ -579,7 +579,7 @@ class Model_model_events_py : public amici::Model_ODE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "63f2e78a186ee722a79b0da2797e0c116ec81065";
+        return "c4a395b6462dc1f26ad1472e6b2fa71fe0d54382";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_events_py/w.h b/models/model_events_py/w.h
deleted file mode 100644
index c399413a8c..0000000000
--- a/models/model_events_py/w.h
+++ /dev/null
@@ -1 +0,0 @@
-#define flux_r0 w[0]
diff --git a/models/model_events_py/xdot.cpp b/models/model_events_py/xdot.cpp
index cc32dc0c59..52dd31b05a 100644
--- a/models/model_events_py/xdot.cpp
+++ b/models/model_events_py/xdot.cpp
@@ -6,7 +6,6 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 
 namespace amici {
diff --git a/models/model_jakstat_adjoint_py/dxdotdw.cpp b/models/model_jakstat_adjoint_py/dxdotdw.cpp
index 6d8582496d..034db3a925 100644
--- a/models/model_jakstat_adjoint_py/dxdotdw.cpp
+++ b/models/model_jakstat_adjoint_py/dxdotdw.cpp
@@ -7,8 +7,8 @@
 namespace amici {
 namespace model_model_jakstat_adjoint_py {
 
-static constexpr std::array<sunindextype, 3> dxdotdw_colptrs_model_jakstat_adjoint_py_ = {
-    0, 2, 2
+static constexpr std::array<sunindextype, 2> dxdotdw_colptrs_model_jakstat_adjoint_py_ = {
+    0, 2
 };
 
 void dxdotdw_colptrs_model_jakstat_adjoint_py(SUNMatrixWrapper &dxdotdw){
diff --git a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.cpp b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.cpp
index 1d933731aa..85d326e72d 100644
--- a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.cpp
+++ b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.cpp
@@ -54,9 +54,8 @@ ObservableScaling::lin, // y[1]
 ObservableScaling::lin, // y[2]
 };
 
-std::array<const char*, 2> expression_names = {
+std::array<const char*, 1> expression_names = {
     "u", // w[0]
-"flux_r0", // w[1]
 };
 
 std::array<const char*, 17> free_parameter_ids = {
@@ -102,9 +101,8 @@ std::array<const char*, 3> observable_ids = {
 "obs_spline", // y[2]
 };
 
-std::array<const char*, 2> expression_ids = {
+std::array<const char*, 1> expression_ids = {
     "u", // w[0]
-"flux_r0", // w[1]
 };
 
 std::array<int, 9> state_idxs_solver = {
diff --git a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h
index bdad8fe252..81571c4e9a 100644
--- a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h
+++ b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h
@@ -19,12 +19,12 @@ extern std::array<const char*, 2> fixed_parameter_names;
 extern std::array<const char*, 9> state_names;
 extern std::array<const char*, 3> observable_names;
 extern std::array<const ObservableScaling, 3> observable_scalings;
-extern std::array<const char*, 2> expression_names;
+extern std::array<const char*, 1> expression_names;
 extern std::array<const char*, 17> free_parameter_ids;
 extern std::array<const char*, 2> fixed_parameter_ids;
 extern std::array<const char*, 9> state_ids;
 extern std::array<const char*, 3> observable_ids;
-extern std::array<const char*, 2> expression_ids;
+extern std::array<const char*, 1> expression_ids;
 extern std::array<int, 9> state_idxs_solver;
 
 extern void Jy_model_jakstat_adjoint_py(realtype *Jy, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my);
@@ -125,7 +125,7 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE {
                   .ne = 0,
                   .ne_solver = 0,
                   .nspl = 1,
-                  .nw = 2,
+                  .nw = 1,
                   .ndwdx = 0,
                   .ndwdp = 5,
                   .ndwdw = 0,
@@ -552,7 +552,7 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "63f2e78a186ee722a79b0da2797e0c116ec81065";
+        return "c4a395b6462dc1f26ad1472e6b2fa71fe0d54382";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_jakstat_adjoint_py/w.h b/models/model_jakstat_adjoint_py/w.h
index c3be8d685c..bcc7e79109 100644
--- a/models/model_jakstat_adjoint_py/w.h
+++ b/models/model_jakstat_adjoint_py/w.h
@@ -1,2 +1 @@
 #define u w[0]
-#define flux_r0 w[1]
diff --git a/models/model_nested_events_py/CMakeLists.txt b/models/model_nested_events_py/CMakeLists.txt
index 52f7dd707d..9db2b46bc6 100644
--- a/models/model_nested_events_py/CMakeLists.txt
+++ b/models/model_nested_events_py/CMakeLists.txt
@@ -90,7 +90,6 @@ stau.cpp
 stau.h 
 sx.h 
 sx0.cpp 
-w.h 
 wrapfunctions.cpp 
 wrapfunctions.h 
 x.h 
diff --git a/models/model_nested_events_py/deltaqB.cpp b/models/model_nested_events_py/deltaqB.cpp
index b30e9fcfcb..8057560e4b 100644
--- a/models/model_nested_events_py/deltaqB.cpp
+++ b/models/model_nested_events_py/deltaqB.cpp
@@ -5,7 +5,6 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
diff --git a/models/model_nested_events_py/deltaxB.cpp b/models/model_nested_events_py/deltaxB.cpp
index 1cdf5b5eb6..2c1aafb290 100644
--- a/models/model_nested_events_py/deltaxB.cpp
+++ b/models/model_nested_events_py/deltaxB.cpp
@@ -5,7 +5,6 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
diff --git a/models/model_nested_events_py/dxdotdp_explicit.cpp b/models/model_nested_events_py/dxdotdp_explicit.cpp
index f265afc06f..b0c70d89f7 100644
--- a/models/model_nested_events_py/dxdotdp_explicit.cpp
+++ b/models/model_nested_events_py/dxdotdp_explicit.cpp
@@ -48,7 +48,6 @@ void dxdotdp_explicit_rowvals_model_nested_events_py(SUNMatrixWrapper &dxdotdp_e
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "dxdotdp_explicit.h"
 
 namespace amici {
diff --git a/models/model_nested_events_py/dxdotdx_explicit.cpp b/models/model_nested_events_py/dxdotdx_explicit.cpp
index c613245501..15604f3416 100644
--- a/models/model_nested_events_py/dxdotdx_explicit.cpp
+++ b/models/model_nested_events_py/dxdotdx_explicit.cpp
@@ -48,7 +48,6 @@ void dxdotdx_explicit_rowvals_model_nested_events_py(SUNMatrixWrapper &dxdotdx_e
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "dxdotdx_explicit.h"
 
 namespace amici {
diff --git a/models/model_nested_events_py/model_nested_events_py.cpp b/models/model_nested_events_py/model_nested_events_py.cpp
index 1237c430bf..15df7aa91a 100644
--- a/models/model_nested_events_py/model_nested_events_py.cpp
+++ b/models/model_nested_events_py/model_nested_events_py.cpp
@@ -29,8 +29,8 @@ std::array<const ObservableScaling, 1> observable_scalings = {
     ObservableScaling::lin, // y[0]
 };
 
-std::array<const char*, 1> expression_names = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_names = {
+    
 };
 
 std::array<const char*, 5> free_parameter_ids = {
@@ -53,8 +53,8 @@ std::array<const char*, 1> observable_ids = {
     "obs_Virus", // y[0]
 };
 
-std::array<const char*, 1> expression_ids = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_ids = {
+    
 };
 
 std::array<int, 1> state_idxs_solver = {
diff --git a/models/model_nested_events_py/model_nested_events_py.h b/models/model_nested_events_py/model_nested_events_py.h
index 8f19547b87..cb8c654a72 100644
--- a/models/model_nested_events_py/model_nested_events_py.h
+++ b/models/model_nested_events_py/model_nested_events_py.h
@@ -19,12 +19,12 @@ extern std::array<const char*, 0> fixed_parameter_names;
 extern std::array<const char*, 1> state_names;
 extern std::array<const char*, 1> observable_names;
 extern std::array<const ObservableScaling, 1> observable_scalings;
-extern std::array<const char*, 1> expression_names;
+extern std::array<const char*, 0> expression_names;
 extern std::array<const char*, 5> free_parameter_ids;
 extern std::array<const char*, 0> fixed_parameter_ids;
 extern std::array<const char*, 1> state_ids;
 extern std::array<const char*, 1> observable_ids;
-extern std::array<const char*, 1> expression_ids;
+extern std::array<const char*, 0> expression_ids;
 extern std::array<int, 1> state_idxs_solver;
 
 extern void Jy_model_nested_events_py(realtype *Jy, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my);
@@ -125,7 +125,7 @@ class Model_model_nested_events_py : public amici::Model_ODE {
                   .ne = 3,
                   .ne_solver = 2,
                   .nspl = 0,
-                  .nw = 1,
+                  .nw = 0,
                   .ndwdx = 0,
                   .ndwdp = 0,
                   .ndwdw = 0,
@@ -552,7 +552,7 @@ class Model_model_nested_events_py : public amici::Model_ODE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "63f2e78a186ee722a79b0da2797e0c116ec81065";
+        return "c4a395b6462dc1f26ad1472e6b2fa71fe0d54382";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_nested_events_py/w.h b/models/model_nested_events_py/w.h
deleted file mode 100644
index c399413a8c..0000000000
--- a/models/model_nested_events_py/w.h
+++ /dev/null
@@ -1 +0,0 @@
-#define flux_r0 w[0]
diff --git a/models/model_nested_events_py/xdot.cpp b/models/model_nested_events_py/xdot.cpp
index ae6ea44dd4..77c6a2a687 100644
--- a/models/model_nested_events_py/xdot.cpp
+++ b/models/model_nested_events_py/xdot.cpp
@@ -5,7 +5,6 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 
 namespace amici {
diff --git a/models/model_neuron_py/CMakeLists.txt b/models/model_neuron_py/CMakeLists.txt
index 52e9275e60..12295b4bc7 100644
--- a/models/model_neuron_py/CMakeLists.txt
+++ b/models/model_neuron_py/CMakeLists.txt
@@ -104,7 +104,6 @@ stau.h
 sx.h 
 sx0.cpp 
 sx0_fixedParameters.cpp 
-w.h 
 wrapfunctions.cpp 
 wrapfunctions.h 
 x.h 
diff --git a/models/model_neuron_py/deltaqB.cpp b/models/model_neuron_py/deltaqB.cpp
index 60d5fb86b0..8b7e128d05 100644
--- a/models/model_neuron_py/deltaqB.cpp
+++ b/models/model_neuron_py/deltaqB.cpp
@@ -6,7 +6,6 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
diff --git a/models/model_neuron_py/deltaxB.cpp b/models/model_neuron_py/deltaxB.cpp
index f12934bed8..321a2cc1ca 100644
--- a/models/model_neuron_py/deltaxB.cpp
+++ b/models/model_neuron_py/deltaxB.cpp
@@ -6,7 +6,6 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
diff --git a/models/model_neuron_py/dxdotdp_explicit.cpp b/models/model_neuron_py/dxdotdp_explicit.cpp
index 8ffbe4d46e..39f2601ca6 100644
--- a/models/model_neuron_py/dxdotdp_explicit.cpp
+++ b/models/model_neuron_py/dxdotdp_explicit.cpp
@@ -49,7 +49,6 @@ void dxdotdp_explicit_rowvals_model_neuron_py(SUNMatrixWrapper &dxdotdp_explicit
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "dxdotdp_explicit.h"
 
 namespace amici {
diff --git a/models/model_neuron_py/dxdotdx_explicit.cpp b/models/model_neuron_py/dxdotdx_explicit.cpp
index 8b632d03ba..8f0df86361 100644
--- a/models/model_neuron_py/dxdotdx_explicit.cpp
+++ b/models/model_neuron_py/dxdotdx_explicit.cpp
@@ -49,7 +49,6 @@ void dxdotdx_explicit_rowvals_model_neuron_py(SUNMatrixWrapper &dxdotdx_explicit
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "dxdotdx_explicit.h"
 
 namespace amici {
diff --git a/models/model_neuron_py/model_neuron_py.cpp b/models/model_neuron_py/model_neuron_py.cpp
index 70995ad2e2..c9ba737562 100644
--- a/models/model_neuron_py/model_neuron_py.cpp
+++ b/models/model_neuron_py/model_neuron_py.cpp
@@ -30,8 +30,8 @@ std::array<const ObservableScaling, 1> observable_scalings = {
     ObservableScaling::lin, // y[0]
 };
 
-std::array<const char*, 1> expression_names = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_names = {
+    
 };
 
 std::array<const char*, 4> free_parameter_ids = {
@@ -55,8 +55,8 @@ std::array<const char*, 1> observable_ids = {
     "y1", // y[0]
 };
 
-std::array<const char*, 1> expression_ids = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_ids = {
+    
 };
 
 std::array<int, 2> state_idxs_solver = {
diff --git a/models/model_neuron_py/model_neuron_py.h b/models/model_neuron_py/model_neuron_py.h
index 8f27330f47..6327f49d2c 100644
--- a/models/model_neuron_py/model_neuron_py.h
+++ b/models/model_neuron_py/model_neuron_py.h
@@ -19,12 +19,12 @@ extern std::array<const char*, 2> fixed_parameter_names;
 extern std::array<const char*, 2> state_names;
 extern std::array<const char*, 1> observable_names;
 extern std::array<const ObservableScaling, 1> observable_scalings;
-extern std::array<const char*, 1> expression_names;
+extern std::array<const char*, 0> expression_names;
 extern std::array<const char*, 4> free_parameter_ids;
 extern std::array<const char*, 2> fixed_parameter_ids;
 extern std::array<const char*, 2> state_ids;
 extern std::array<const char*, 1> observable_ids;
-extern std::array<const char*, 1> expression_ids;
+extern std::array<const char*, 0> expression_ids;
 extern std::array<int, 2> state_idxs_solver;
 
 extern void Jy_model_neuron_py(realtype *Jy, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my);
@@ -125,7 +125,7 @@ class Model_model_neuron_py : public amici::Model_ODE {
                   .ne = 1,
                   .ne_solver = 1,
                   .nspl = 0,
-                  .nw = 1,
+                  .nw = 0,
                   .ndwdx = 0,
                   .ndwdp = 0,
                   .ndwdw = 0,
@@ -574,7 +574,7 @@ class Model_model_neuron_py : public amici::Model_ODE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "63f2e78a186ee722a79b0da2797e0c116ec81065";
+        return "c4a395b6462dc1f26ad1472e6b2fa71fe0d54382";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_neuron_py/w.h b/models/model_neuron_py/w.h
deleted file mode 100644
index c399413a8c..0000000000
--- a/models/model_neuron_py/w.h
+++ /dev/null
@@ -1 +0,0 @@
-#define flux_r0 w[0]
diff --git a/models/model_neuron_py/xdot.cpp b/models/model_neuron_py/xdot.cpp
index cd1a50688e..f188d8ff09 100644
--- a/models/model_neuron_py/xdot.cpp
+++ b/models/model_neuron_py/xdot.cpp
@@ -6,7 +6,6 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
-#include "w.h"
 #include "xdot.h"
 
 namespace amici {
diff --git a/models/model_robertson_py/CMakeLists.txt b/models/model_robertson_py/CMakeLists.txt
index 3e7685365c..fbab6a9f6b 100644
--- a/models/model_robertson_py/CMakeLists.txt
+++ b/models/model_robertson_py/CMakeLists.txt
@@ -82,7 +82,6 @@ p.h
 sigmay.cpp 
 sigmay.h 
 sx0_fixedParameters.cpp 
-w.h 
 wrapfunctions.cpp 
 wrapfunctions.h 
 x.h 
diff --git a/models/model_robertson_py/dxdotdp_explicit.cpp b/models/model_robertson_py/dxdotdp_explicit.cpp
index ce679cb43b..d8193178ac 100644
--- a/models/model_robertson_py/dxdotdp_explicit.cpp
+++ b/models/model_robertson_py/dxdotdp_explicit.cpp
@@ -49,7 +49,6 @@ void dxdotdp_explicit_rowvals_model_robertson_py(SUNMatrixWrapper &dxdotdp_expli
 #include "p.h"
 #include "k.h"
 #include "dx.h"
-#include "w.h"
 #include "dxdotdp_explicit.h"
 
 namespace amici {
diff --git a/models/model_robertson_py/dxdotdx_explicit.cpp b/models/model_robertson_py/dxdotdx_explicit.cpp
index 52605f51f3..f803addb5c 100644
--- a/models/model_robertson_py/dxdotdx_explicit.cpp
+++ b/models/model_robertson_py/dxdotdx_explicit.cpp
@@ -49,7 +49,6 @@ void dxdotdx_explicit_rowvals_model_robertson_py(SUNMatrixWrapper &dxdotdx_expli
 #include "p.h"
 #include "k.h"
 #include "dx.h"
-#include "w.h"
 #include "dxdotdx_explicit.h"
 
 namespace amici {
diff --git a/models/model_robertson_py/model_robertson_py.cpp b/models/model_robertson_py/model_robertson_py.cpp
index 145282b8b9..fa7e441dfb 100644
--- a/models/model_robertson_py/model_robertson_py.cpp
+++ b/models/model_robertson_py/model_robertson_py.cpp
@@ -33,8 +33,8 @@ ObservableScaling::lin, // y[1]
 ObservableScaling::lin, // y[2]
 };
 
-std::array<const char*, 1> expression_names = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_names = {
+    
 };
 
 std::array<const char*, 3> free_parameter_ids = {
@@ -59,8 +59,8 @@ std::array<const char*, 3> observable_ids = {
 "obs_x3", // y[2]
 };
 
-std::array<const char*, 1> expression_ids = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_ids = {
+    
 };
 
 std::array<int, 3> state_idxs_solver = {
diff --git a/models/model_robertson_py/model_robertson_py.h b/models/model_robertson_py/model_robertson_py.h
index 30b1078026..26640ca340 100644
--- a/models/model_robertson_py/model_robertson_py.h
+++ b/models/model_robertson_py/model_robertson_py.h
@@ -19,12 +19,12 @@ extern std::array<const char*, 1> fixed_parameter_names;
 extern std::array<const char*, 3> state_names;
 extern std::array<const char*, 3> observable_names;
 extern std::array<const ObservableScaling, 3> observable_scalings;
-extern std::array<const char*, 1> expression_names;
+extern std::array<const char*, 0> expression_names;
 extern std::array<const char*, 3> free_parameter_ids;
 extern std::array<const char*, 1> fixed_parameter_ids;
 extern std::array<const char*, 3> state_ids;
 extern std::array<const char*, 3> observable_ids;
-extern std::array<const char*, 1> expression_ids;
+extern std::array<const char*, 0> expression_ids;
 extern std::array<int, 3> state_idxs_solver;
 
 extern void Jy_model_robertson_py(realtype *Jy, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my);
@@ -125,7 +125,7 @@ class Model_model_robertson_py : public amici::Model_DAE {
                   .ne = 0,
                   .ne_solver = 0,
                   .nspl = 0,
-                  .nw = 1,
+                  .nw = 0,
                   .ndwdx = 0,
                   .ndwdp = 0,
                   .ndwdw = 0,
@@ -536,7 +536,7 @@ class Model_model_robertson_py : public amici::Model_DAE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "63f2e78a186ee722a79b0da2797e0c116ec81065";
+        return "c4a395b6462dc1f26ad1472e6b2fa71fe0d54382";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_robertson_py/w.h b/models/model_robertson_py/w.h
deleted file mode 100644
index c399413a8c..0000000000
--- a/models/model_robertson_py/w.h
+++ /dev/null
@@ -1 +0,0 @@
-#define flux_r0 w[0]
diff --git a/models/model_robertson_py/xdot.cpp b/models/model_robertson_py/xdot.cpp
index cc6b9808d2..73c3888787 100644
--- a/models/model_robertson_py/xdot.cpp
+++ b/models/model_robertson_py/xdot.cpp
@@ -6,7 +6,6 @@
 #include "p.h"
 #include "k.h"
 #include "dx.h"
-#include "w.h"
 #include "xdot.h"
 
 namespace amici {
diff --git a/models/model_steadystate_py/CMakeLists.txt b/models/model_steadystate_py/CMakeLists.txt
index f1c1e1492e..2e7e7e7bcd 100644
--- a/models/model_steadystate_py/CMakeLists.txt
+++ b/models/model_steadystate_py/CMakeLists.txt
@@ -81,7 +81,6 @@ p.h
 sigmay.cpp 
 sigmay.h 
 sx0_fixedParameters.cpp 
-w.h 
 wrapfunctions.cpp 
 wrapfunctions.h 
 x.h 
diff --git a/models/model_steadystate_py/dxdotdp_explicit.cpp b/models/model_steadystate_py/dxdotdp_explicit.cpp
index a9be02fda2..1051fc4b9c 100644
--- a/models/model_steadystate_py/dxdotdp_explicit.cpp
+++ b/models/model_steadystate_py/dxdotdp_explicit.cpp
@@ -48,7 +48,6 @@ void dxdotdp_explicit_rowvals_model_steadystate_py(SUNMatrixWrapper &dxdotdp_exp
 #include "x.h"
 #include "p.h"
 #include "k.h"
-#include "w.h"
 #include "dxdotdp_explicit.h"
 
 namespace amici {
diff --git a/models/model_steadystate_py/dxdotdx_explicit.cpp b/models/model_steadystate_py/dxdotdx_explicit.cpp
index 870d3bc63e..7606175cad 100644
--- a/models/model_steadystate_py/dxdotdx_explicit.cpp
+++ b/models/model_steadystate_py/dxdotdx_explicit.cpp
@@ -48,7 +48,6 @@ void dxdotdx_explicit_rowvals_model_steadystate_py(SUNMatrixWrapper &dxdotdx_exp
 #include "x.h"
 #include "p.h"
 #include "k.h"
-#include "w.h"
 #include "dxdotdx_explicit.h"
 
 namespace amici {
diff --git a/models/model_steadystate_py/model_steadystate_py.cpp b/models/model_steadystate_py/model_steadystate_py.cpp
index 58da853481..65de9fee33 100644
--- a/models/model_steadystate_py/model_steadystate_py.cpp
+++ b/models/model_steadystate_py/model_steadystate_py.cpp
@@ -38,8 +38,8 @@ ObservableScaling::lin, // y[1]
 ObservableScaling::lin, // y[2]
 };
 
-std::array<const char*, 1> expression_names = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_names = {
+    
 };
 
 std::array<const char*, 5> free_parameter_ids = {
@@ -69,8 +69,8 @@ std::array<const char*, 3> observable_ids = {
 "obs_x3", // y[2]
 };
 
-std::array<const char*, 1> expression_ids = {
-    "flux_r0", // w[0]
+std::array<const char*, 0> expression_ids = {
+    
 };
 
 std::array<int, 3> state_idxs_solver = {
diff --git a/models/model_steadystate_py/model_steadystate_py.h b/models/model_steadystate_py/model_steadystate_py.h
index c7d6c7f0af..dcbbbec161 100644
--- a/models/model_steadystate_py/model_steadystate_py.h
+++ b/models/model_steadystate_py/model_steadystate_py.h
@@ -19,12 +19,12 @@ extern std::array<const char*, 4> fixed_parameter_names;
 extern std::array<const char*, 3> state_names;
 extern std::array<const char*, 3> observable_names;
 extern std::array<const ObservableScaling, 3> observable_scalings;
-extern std::array<const char*, 1> expression_names;
+extern std::array<const char*, 0> expression_names;
 extern std::array<const char*, 5> free_parameter_ids;
 extern std::array<const char*, 4> fixed_parameter_ids;
 extern std::array<const char*, 3> state_ids;
 extern std::array<const char*, 3> observable_ids;
-extern std::array<const char*, 1> expression_ids;
+extern std::array<const char*, 0> expression_ids;
 extern std::array<int, 3> state_idxs_solver;
 
 extern void Jy_model_steadystate_py(realtype *Jy, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my);
@@ -125,7 +125,7 @@ class Model_model_steadystate_py : public amici::Model_ODE {
                   .ne = 0,
                   .ne_solver = 0,
                   .nspl = 0,
-                  .nw = 1,
+                  .nw = 0,
                   .ndwdx = 0,
                   .ndwdp = 0,
                   .ndwdw = 0,
@@ -536,7 +536,7 @@ class Model_model_steadystate_py : public amici::Model_ODE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "63f2e78a186ee722a79b0da2797e0c116ec81065";
+        return "c4a395b6462dc1f26ad1472e6b2fa71fe0d54382";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_steadystate_py/w.h b/models/model_steadystate_py/w.h
deleted file mode 100644
index c399413a8c..0000000000
--- a/models/model_steadystate_py/w.h
+++ /dev/null
@@ -1 +0,0 @@
-#define flux_r0 w[0]
diff --git a/models/model_steadystate_py/xdot.cpp b/models/model_steadystate_py/xdot.cpp
index 6b96b07288..9d25458d71 100644
--- a/models/model_steadystate_py/xdot.cpp
+++ b/models/model_steadystate_py/xdot.cpp
@@ -5,7 +5,6 @@
 #include "x.h"
 #include "p.h"
 #include "k.h"
-#include "w.h"
 #include "xdot.h"
 
 namespace amici {
diff --git a/python/sdist/amici/exporters/sundials/de_export.py b/python/sdist/amici/exporters/sundials/de_export.py
index 4fdfc43feb..0004ff1d26 100644
--- a/python/sdist/amici/exporters/sundials/de_export.py
+++ b/python/sdist/amici/exporters/sundials/de_export.py
@@ -427,12 +427,13 @@ def _write_function_file(self, function: str) -> None:
             if sym not in self.model.sym_names():
                 continue
 
+            if len(self.model.sym(sym)) == 0:
+                continue
+
             if sym in sparse_functions:
                 iszero = smart_is_zero_matrix(self.model.sparseeq(sym))
             elif sym in self.functions:
                 iszero = smart_is_zero_matrix(self.model.eq(sym))
-            elif len(self.model.sym(sym)) == 0:
-                continue
             else:
                 iszero = False
 
diff --git a/python/sdist/amici/importers/sbml/__init__.py b/python/sdist/amici/importers/sbml/__init__.py
index deb0b475d2..88f6d4208a 100644
--- a/python/sdist/amici/importers/sbml/__init__.py
+++ b/python/sdist/amici/importers/sbml/__init__.py
@@ -627,8 +627,12 @@ def _build_ode_model(
         ]
 
         # create dynamics without respecting conservation laws first
-        dxdt = smart_multiply(
-            self.stoichiometric_matrix, MutableDenseMatrix(fluxes)
+        dxdt = (
+            smart_multiply(
+                self.stoichiometric_matrix, MutableDenseMatrix(fluxes)
+            )
+            if len(self.stoichiometric_matrix)
+            else sp.zeros(self.stoichiometric_matrix.shape[0], 1)
         )
         # dxdt has algebraic states at the end
         assert dxdt.shape[0] - len(self.symbols[SymbolId.SPECIES]) == len(
@@ -1348,10 +1352,7 @@ def _process_reactions(self):
         Get reactions from SBML model.
         """
         reactions = self.sbml.getListOfReactions()
-        # nr (number of reactions) should have a minimum length of 1. This is
-        # to ensure that, if there are no reactions, the stoichiometric matrix
-        # and flux vector multiply to a zero vector with dimensions (nx, 1).
-        nr = max(1, len(reactions))
+        nr = len(reactions)
         nx = len(self.symbols[SymbolId.SPECIES])
         # stoichiometric matrix
         self.stoichiometric_matrix = sp.SparseMatrix(sp.zeros(nx, nr))
@@ -1364,7 +1365,7 @@ def _process_reactions(self):
             if reaction.isSetId()
             else f"flux_r{reaction_idx}"
             for reaction_idx, reaction in enumerate(reactions)
-        ] or ["flux_r0"]
+        ]
 
         reaction_ids = [
             reaction.getId() for reaction in reactions if reaction.isSetId()
@@ -1588,14 +1589,21 @@ def _process_rule_algebraic(self, rule: libsbml.AlgebraicRule):
                 ):
                     # placeholder, needs to be determined in IC calculation
                     symbol["init"] = sp.Float(0.0)
-                self.stoichiometric_matrix = (
-                    self.stoichiometric_matrix.row_insert(
-                        self.stoichiometric_matrix.shape[0],
-                        sp.SparseMatrix(
-                            [[0] * self.stoichiometric_matrix.shape[1]]
-                        ),
+                if len(self.stoichiometric_matrix):
+                    self.stoichiometric_matrix = (
+                        self.stoichiometric_matrix.row_insert(
+                            self.stoichiometric_matrix.shape[0],
+                            sp.SparseMatrix(
+                                1, self.stoichiometric_matrix.shape[1], {}
+                            ),
+                        )
+                    )
+                else:
+                    # we can't use Matrix.row_insert for (n, 0) matrices
+                    #  as this will ignore the previous number of rows
+                    self.stoichiometric_matrix = sp.SparseMatrix(
+                        self.stoichiometric_matrix.shape[0] + 1, 0, {}
                     )
-                )
             elif var_ix != self.stoichiometric_matrix.shape[0] - 1:
                 # if not the last col, move it to the end
                 # as we reorder state variables
@@ -2357,9 +2365,10 @@ def _process_conservation_laws(self, ode_model: DEModel) -> None:
             species_solver = list(range(ode_model.num_states_rdata()))
 
         # prune out species from stoichiometry and
-        self.stoichiometric_matrix = self.stoichiometric_matrix[
-            species_solver, :
-        ]
+        if len(self.stoichiometric_matrix):
+            self.stoichiometric_matrix = self.stoichiometric_matrix[
+                species_solver, :
+            ]
 
         # add the found CLs to the ode_model
         for cl in conservation_laws:
@@ -2379,9 +2388,11 @@ def _get_conservation_laws_demartino(
             quantity (including the eliminated one)
             (2) coefficients for the species in (1)
         """
-        from .conserved_quantities_demartino import (
-            compute_moiety_conservation_laws,
-        )
+        if not len(self.stoichiometric_matrix) or not len(
+            self.symbols[SymbolId.SPECIES]
+        ):
+            # no species or no reactions -> no conservation laws
+            return []
 
         sm = self.stoichiometric_matrix[
             : len(self.symbols[SymbolId.SPECIES]), :
@@ -2405,6 +2416,10 @@ def _get_conservation_laws_demartino(
         if not _non_const_conservation_laws_supported(self.sbml):
             return []
 
+        from .conserved_quantities_demartino import (
+            compute_moiety_conservation_laws,
+        )
+
         cls_state_idxs, cls_coefficients = compute_moiety_conservation_laws(
             stoichiometric_list,
             *sm.shape,
@@ -2453,6 +2468,12 @@ def _get_conservation_laws_rref(
             quantity (including the eliminated one)
             (2) coefficients for the species in (1)
         """
+        if not len(self.stoichiometric_matrix) or not len(
+            self.symbols[SymbolId.SPECIES]
+        ):
+            # no species or no reactions -> no conservation laws
+            return []
+
         import numpy as np
         from numpy.linalg import matrix_rank
 
diff --git a/tests/cpp/expected_results_py.h5 b/tests/cpp/expected_results_py.h5
index 1428f7e723..bd8a7f3a5d 100644
--- a/tests/cpp/expected_results_py.h5
+++ b/tests/cpp/expected_results_py.h5
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:eae94de78f59eb136564483b1b713199743593f0f9bddbab7dc0bcfa5399d4f9
-size 1847292
+oid sha256:e7c4f4ee6f9a508a2cd7b11947d61c1c53d97f6f4cfcf665e7abfbf5d46ac81b
+size 1808960