diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h
index 9466ffac7f..ea105a29ac 100644
--- a/include/amici/abstract_model.h
+++ b/include/amici/abstract_model.h
@@ -299,6 +299,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param tcl total abundances for conservation laws
      * @param sx current state sensitivity
@@ -307,8 +308,9 @@ class AbstractModel {
      */
     virtual void fstau(
         realtype* stau, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx,
-        realtype const* tcl, realtype const* sx, int ip, int ie
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, realtype const* tcl, realtype const* sx, int ip,
+        int ie
     );
 
     /**
@@ -542,6 +544,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param ie event index
      * @param xdot new model right hand side
@@ -552,9 +555,10 @@ class AbstractModel {
      */
     virtual void fdeltaxB(
         realtype* deltaxB, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx, int ie,
-        realtype const* xdot, realtype const* xdot_old, realtype const* x_old,
-        realtype const* xB, realtype const* tcl
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, int ie, realtype const* xdot,
+        realtype const* xdot_old, realtype const* x_old, realtype const* xB,
+        realtype const* tcl
     );
 
     /**
@@ -565,6 +569,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param ip sensitivity index
      * @param ie event index
@@ -575,9 +580,9 @@ class AbstractModel {
      */
     virtual void fdeltaqB(
         realtype* deltaqB, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx, int ip,
-        int ie, realtype const* xdot, realtype const* xdot_old,
-        realtype const* x_old, realtype const* xB
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, int ip, int ie, realtype const* xdot,
+        realtype const* xdot_old, realtype const* x_old, realtype const* xB
     );
 
     /**
@@ -1059,11 +1064,13 @@ class AbstractModel {
      * @brief Compute explicit roots of the model.
      * @param p parameter vector
      * @param k constant vector
+     * @param w vector with helper variables
      * @return A vector of length ne_solver, each containing a vector of
      * explicit roots for the corresponding event.
      */
     virtual std::vector<std::vector<realtype>> fexplicit_roots(
-        [[maybe_unused]] realtype const* p, [[maybe_unused]] realtype const* k
+        [[maybe_unused]] realtype const* p, [[maybe_unused]] realtype const* k,
+        [[maybe_unused]] realtype const* w
     ) = 0;
 };
 
diff --git a/include/amici/model.h b/include/amici/model.h
index fe15c1066e..88181d3897 100644
--- a/include/amici/model.h
+++ b/include/amici/model.h
@@ -1591,7 +1591,8 @@ class Model : public AbstractModel, public ModelDimensions {
     }
 
     [[nodiscard]] std::vector<std::vector<realtype>> fexplicit_roots(
-        [[maybe_unused]] realtype const* p, [[maybe_unused]] realtype const* k
+        [[maybe_unused]] realtype const* p, [[maybe_unused]] realtype const* k,
+        [[maybe_unused]] realtype const* w
     ) override {
         if (ne != ne_solver) {
             throw AmiException(
diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h
index 66f6559a1d..47fe33ffdf 100644
--- a/include/amici/model_dae.h
+++ b/include/amici/model_dae.h
@@ -359,11 +359,13 @@ class Model_DAE : public Model {
      * @param p parameter vector
      * @param k constants vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx Vector with the derivative states
      **/
     virtual void froot(
         realtype* root, realtype t, realtype const* x, double const* p,
-        double const* k, realtype const* h, realtype const* dx
+        double const* k, realtype const* h, realtype const* w,
+        realtype const* dx
     );
 
     /**
@@ -444,7 +446,7 @@ class Model_DAE : public Model {
      * @param x Vector with the states
      * @param p parameter vector
      * @param k constants vector
-     * @param h heavyside vector
+     * @param h Heaviside vector
      * @param dx Vector with the derivative states
      * @param w vector with helper variables
      */
diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h
index b3cb844cd8..5dff674169 100644
--- a/include/amici/model_ode.h
+++ b/include/amici/model_ode.h
@@ -316,11 +316,13 @@ class Model_ODE : public Model {
      * @param p parameter vector
      * @param k constants vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param tcl total abundances for conservation laws
      **/
     virtual void froot(
         realtype* root, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* tcl
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* tcl
     );
 
     /**
diff --git a/models/model_calvetti_py/explicit_roots.cpp b/models/model_calvetti_py/explicit_roots.cpp
index 4598f5ec00..76e9d11849 100644
--- a/models/model_calvetti_py/explicit_roots.cpp
+++ b/models/model_calvetti_py/explicit_roots.cpp
@@ -4,11 +4,12 @@
 #include <algorithm>
 #include <vector>
 #include "k.h"
+#include "w.h"
 
 namespace amici {
 namespace model_model_calvetti_py {
 
-std::vector<std::vector<realtype>> explicit_roots_model_calvetti_py(const realtype *p, const realtype *k){
+std::vector<std::vector<realtype>> explicit_roots_model_calvetti_py(const realtype *p, const realtype *k, const realtype *w){
     return {
         {10},
         {10},
diff --git a/models/model_calvetti_py/model_calvetti_py.h b/models/model_calvetti_py/model_calvetti_py.h
index 4d49018ad2..bf5ea1242f 100644
--- a/models/model_calvetti_py/model_calvetti_py.h
+++ b/models/model_calvetti_py/model_calvetti_py.h
@@ -38,7 +38,7 @@ extern void dJydy_rowvals_model_calvetti_py(SUNMatrixWrapper &rowvals, int index
 
 
 
-extern void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
+extern void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 
 
 
@@ -99,7 +99,7 @@ extern void x_solver_model_calvetti_py(realtype *x_solver, const realtype *x_rda
 extern std::vector<HermiteSpline> create_splines_model_calvetti_py(const realtype *p, const realtype *k);
 
 
-extern std::vector<std::vector<realtype>> explicit_roots_model_calvetti_py(const realtype *p, const realtype *k);
+extern std::vector<std::vector<realtype>> explicit_roots_model_calvetti_py(const realtype *p, const realtype *k, const realtype *w);
 /**
  * @brief AMICI-generated model subclass.
  */
@@ -200,10 +200,10 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fdeltasx(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old) override {}
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
 
 
     void fdrzdp(realtype *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {}
@@ -323,8 +323,8 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {
-        root_model_calvetti_py(root, t, x, p, k, h, tcl);
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {
+        root_model_calvetti_py(root, t, x, p, k, h, w, tcl);
     }
 
 
@@ -339,7 +339,7 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
 
@@ -411,8 +411,8 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {}
 
 
-    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k) override {
-        return explicit_roots_model_calvetti_py(p, k);
+    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k, const realtype *w) override {
+        return explicit_roots_model_calvetti_py(p, k, w);
     }
 
 
@@ -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 "40190b46b1b398e321314ded4169fe910b37c484";
+        return "3fb84cd5df12639f17b179d681e8ba4b5be8a160";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_calvetti_py/root.cpp b/models/model_calvetti_py/root.cpp
index 65b8c9dc1f..07327917db 100644
--- a/models/model_calvetti_py/root.cpp
+++ b/models/model_calvetti_py/root.cpp
@@ -5,11 +5,12 @@
 #include "x.h"
 #include "k.h"
 #include "h.h"
+#include "w.h"
 
 namespace amici {
 namespace model_model_calvetti_py {
 
-void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){
+void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl){
     root[0] = t - 10;
     root[1] = 10 - t;
     root[2] = 12 - t;
diff --git a/models/model_dirac_py/deltaqB.cpp b/models/model_dirac_py/deltaqB.cpp
index aa3c488c56..b95bd81ee1 100644
--- a/models/model_dirac_py/deltaqB.cpp
+++ b/models/model_dirac_py/deltaqB.cpp
@@ -5,6 +5,7 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -13,7 +14,7 @@
 namespace amici {
 namespace model_model_dirac_py {
 
-void deltaqB_model_dirac_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
+void deltaqB_model_dirac_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
     switch(ie) {
         case 0:
             switch(ip) {
diff --git a/models/model_dirac_py/explicit_roots.cpp b/models/model_dirac_py/explicit_roots.cpp
index 5061cc9fd9..e7074251fb 100644
--- a/models/model_dirac_py/explicit_roots.cpp
+++ b/models/model_dirac_py/explicit_roots.cpp
@@ -8,7 +8,7 @@
 namespace amici {
 namespace model_model_dirac_py {
 
-std::vector<std::vector<realtype>> explicit_roots_model_dirac_py(const realtype *p, const realtype *k){
+std::vector<std::vector<realtype>> explicit_roots_model_dirac_py(const realtype *p, const realtype *k, const realtype *w){
     return {
         {p2}
     };
diff --git a/models/model_dirac_py/model_dirac_py.h b/models/model_dirac_py/model_dirac_py.h
index a078a48ff5..b4b4bc3034 100644
--- a/models/model_dirac_py/model_dirac_py.h
+++ b/models/model_dirac_py/model_dirac_py.h
@@ -38,7 +38,7 @@ extern void dJydy_rowvals_model_dirac_py(SUNMatrixWrapper &rowvals, int index);
 
 
 
-extern void root_model_dirac_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
+extern void root_model_dirac_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 
 
 
@@ -77,11 +77,11 @@ extern void xdot_model_dirac_py(realtype *xdot, const realtype t, const realtype
 extern void y_model_dirac_py(realtype *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w);
 
 
-extern void stau_model_dirac_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
+extern void stau_model_dirac_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
 extern void deltax_model_dirac_py(double *deltax, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old);
 extern void deltasx_model_dirac_py(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old);
 
-extern void deltaqB_model_dirac_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
+extern void deltaqB_model_dirac_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
 
 extern void x_solver_model_dirac_py(realtype *x_solver, const realtype *x_rdata);
 
@@ -99,7 +99,7 @@ extern void x_solver_model_dirac_py(realtype *x_solver, const realtype *x_rdata)
 extern std::vector<HermiteSpline> create_splines_model_dirac_py(const realtype *p, const realtype *k);
 
 
-extern std::vector<std::vector<realtype>> explicit_roots_model_dirac_py(const realtype *p, const realtype *k);
+extern std::vector<std::vector<realtype>> explicit_roots_model_dirac_py(const realtype *p, const realtype *k, const realtype *w);
 /**
  * @brief AMICI-generated model subclass.
  */
@@ -201,11 +201,11 @@ class Model_model_dirac_py : public amici::Model_ODE {
     }
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
-        deltaqB_model_dirac_py(deltaqB, t, x, p, k, h, dx, ip, ie, xdot, xdot_old, x_old, xB);
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
+        deltaqB_model_dirac_py(deltaqB, t, x, p, k, h, w, dx, ip, ie, xdot, xdot_old, x_old, xB);
     }
 
 
@@ -314,8 +314,8 @@ class Model_model_dirac_py : public amici::Model_ODE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {
-        root_model_dirac_py(root, t, x, p, k, h, tcl);
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {
+        root_model_dirac_py(root, t, x, p, k, h, w, tcl);
     }
 
 
@@ -330,8 +330,8 @@ class Model_model_dirac_py : public amici::Model_ODE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
-        stau_model_dirac_py(stau, t, x, p, k, h, dx, tcl, sx, ip, ie);
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
+        stau_model_dirac_py(stau, t, x, p, k, h, w, dx, tcl, sx, ip, ie);
     }
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
@@ -398,8 +398,8 @@ class Model_model_dirac_py : public amici::Model_ODE {
     void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {}
 
 
-    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k) override {
-        return explicit_roots_model_dirac_py(p, k);
+    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k, const realtype *w) override {
+        return explicit_roots_model_dirac_py(p, k, w);
     }
 
 
@@ -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 "40190b46b1b398e321314ded4169fe910b37c484";
+        return "3fb84cd5df12639f17b179d681e8ba4b5be8a160";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_dirac_py/root.cpp b/models/model_dirac_py/root.cpp
index b435bbe736..f2f2d5137a 100644
--- a/models/model_dirac_py/root.cpp
+++ b/models/model_dirac_py/root.cpp
@@ -9,7 +9,7 @@
 namespace amici {
 namespace model_model_dirac_py {
 
-void root_model_dirac_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){
+void root_model_dirac_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl){
     root[0] = -p2 + t;
 }
 
diff --git a/models/model_dirac_py/stau.cpp b/models/model_dirac_py/stau.cpp
index 399f7c43d1..c485104410 100644
--- a/models/model_dirac_py/stau.cpp
+++ b/models/model_dirac_py/stau.cpp
@@ -10,7 +10,7 @@
 namespace amici {
 namespace model_model_dirac_py {
 
-void stau_model_dirac_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
+void stau_model_dirac_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
     switch(ie) {
         case 0:
             switch(ip) {
diff --git a/models/model_events_py/deltaqB.cpp b/models/model_events_py/deltaqB.cpp
index ace0e672ae..4200430baf 100644
--- a/models/model_events_py/deltaqB.cpp
+++ b/models/model_events_py/deltaqB.cpp
@@ -6,6 +6,7 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -14,7 +15,7 @@
 namespace amici {
 namespace model_model_events_py {
 
-void deltaqB_model_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
+void deltaqB_model_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
     switch(ie) {
         case 2:
         case 3:
diff --git a/models/model_events_py/deltaxB.cpp b/models/model_events_py/deltaxB.cpp
index 3c23a63a37..0409fb5553 100644
--- a/models/model_events_py/deltaxB.cpp
+++ b/models/model_events_py/deltaxB.cpp
@@ -6,6 +6,7 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -14,7 +15,7 @@
 namespace amici {
 namespace model_model_events_py {
 
-void deltaxB_model_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl){
+void deltaxB_model_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl){
     switch(ie) {
         case 0:
             deltaxB[1] = xB0*(dx1dt - xdot_old0)/(Heaviside_4 + p2*x1*std::exp(-1.0/10.0*t) - p3*x2 + x3 - 1) + xB1*(dx2dt - xdot_old1)/(Heaviside_4 + p2*x1*std::exp(-1.0/10.0*t) - p3*x2 + x3 - 1) + xB2*(dx3dt - xdot_old2)/(Heaviside_4 + p2*x1*std::exp(-1.0/10.0*t) - p3*x2 + x3 - 1);
diff --git a/models/model_events_py/explicit_roots.cpp b/models/model_events_py/explicit_roots.cpp
index 4c6d86f0e3..bd12f13ef9 100644
--- a/models/model_events_py/explicit_roots.cpp
+++ b/models/model_events_py/explicit_roots.cpp
@@ -9,7 +9,7 @@
 namespace amici {
 namespace model_model_events_py {
 
-std::vector<std::vector<realtype>> explicit_roots_model_events_py(const realtype *p, const realtype *k){
+std::vector<std::vector<realtype>> explicit_roots_model_events_py(const realtype *p, const realtype *k, const realtype *w){
     return {
         {p4},
         {p4},
diff --git a/models/model_events_py/model_events_py.h b/models/model_events_py/model_events_py.h
index 9b9c519f0f..0bf2ec82f7 100644
--- a/models/model_events_py/model_events_py.h
+++ b/models/model_events_py/model_events_py.h
@@ -38,7 +38,7 @@ extern void dJzdz_model_events_py(realtype *dJzdz, const int iz, const realtype
 extern void Jrz_model_events_py(realtype *Jrz, const int iz, const realtype *p, const realtype *k, const realtype *rz, const realtype *sigmaz);
 extern void dJrzdsigma_model_events_py(realtype *dJrzdsigma, const int iz, const realtype *p, const realtype *k, const realtype *rz, const realtype *sigmaz);
 extern void dJrzdz_model_events_py(realtype *dJrzdz, const int iz, const realtype *p, const realtype *k, const realtype *rz, const realtype *sigmaz);
-extern void root_model_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
+extern void root_model_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 
 
 
@@ -77,11 +77,11 @@ extern void xdot_model_events_py(realtype *xdot, const realtype t, const realtyp
 extern void y_model_events_py(realtype *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w);
 extern void z_model_events_py(realtype *z, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h);
 extern void rz_model_events_py(realtype *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h);
-extern void stau_model_events_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
+extern void stau_model_events_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
 
 extern void deltasx_model_events_py(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old);
-extern void deltaxB_model_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl);
-extern void deltaqB_model_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
+extern void deltaxB_model_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl);
+extern void deltaqB_model_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
 
 extern void x_solver_model_events_py(realtype *x_solver, const realtype *x_rdata);
 
@@ -99,7 +99,7 @@ extern void x_solver_model_events_py(realtype *x_solver, const realtype *x_rdata
 extern std::vector<HermiteSpline> create_splines_model_events_py(const realtype *p, const realtype *k);
 
 
-extern std::vector<std::vector<realtype>> explicit_roots_model_events_py(const realtype *p, const realtype *k);
+extern std::vector<std::vector<realtype>> explicit_roots_model_events_py(const realtype *p, const realtype *k, const realtype *w);
 /**
  * @brief AMICI-generated model subclass.
  */
@@ -216,13 +216,13 @@ class Model_model_events_py : public amici::Model_ODE {
     }
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {
-        deltaxB_model_events_py(deltaxB, t, x, p, k, h, dx, ie, xdot, xdot_old, x_old, xB, tcl);
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {
+        deltaxB_model_events_py(deltaxB, t, x, p, k, h, w, dx, ie, xdot, xdot_old, x_old, xB, tcl);
     }
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
-        deltaqB_model_events_py(deltaqB, t, x, p, k, h, dx, ip, ie, xdot, xdot_old, x_old, xB);
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
+        deltaqB_model_events_py(deltaqB, t, x, p, k, h, w, dx, ip, ie, xdot, xdot_old, x_old, xB);
     }
 
 
@@ -337,8 +337,8 @@ class Model_model_events_py : public amici::Model_ODE {
     }
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {
-        root_model_events_py(root, t, x, p, k, h, tcl);
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {
+        root_model_events_py(root, t, x, p, k, h, w, tcl);
     }
 
 
@@ -357,8 +357,8 @@ class Model_model_events_py : public amici::Model_ODE {
     }
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
-        stau_model_events_py(stau, t, x, p, k, h, dx, tcl, sx, ip, ie);
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
+        stau_model_events_py(stau, t, x, p, k, h, w, dx, tcl, sx, ip, ie);
     }
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
@@ -433,8 +433,8 @@ class Model_model_events_py : public amici::Model_ODE {
     void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {}
 
 
-    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k) override {
-        return explicit_roots_model_events_py(p, k);
+    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k, const realtype *w) override {
+        return explicit_roots_model_events_py(p, k, w);
     }
 
 
@@ -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 "40190b46b1b398e321314ded4169fe910b37c484";
+        return "3fb84cd5df12639f17b179d681e8ba4b5be8a160";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_events_py/root.cpp b/models/model_events_py/root.cpp
index 0e0fa0fe45..30ac963273 100644
--- a/models/model_events_py/root.cpp
+++ b/models/model_events_py/root.cpp
@@ -10,7 +10,7 @@
 namespace amici {
 namespace model_model_events_py {
 
-void root_model_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){
+void root_model_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl){
     root[0] = x2 - x3;
     root[1] = x1 - x3;
     root[2] = p4 - t;
diff --git a/models/model_events_py/stau.cpp b/models/model_events_py/stau.cpp
index d2da6b4bdb..9ffc7e19fe 100644
--- a/models/model_events_py/stau.cpp
+++ b/models/model_events_py/stau.cpp
@@ -11,7 +11,7 @@
 namespace amici {
 namespace model_model_events_py {
 
-void stau_model_events_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
+void stau_model_events_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
     switch(ie) {
         case 0:
             switch(ip) {
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 993b769aeb..71c22b1b1c 100644
--- a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h
+++ b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h
@@ -195,10 +195,10 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE {
     void fdeltasx(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old) override {}
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
 
 
     void fdrzdp(realtype *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {}
@@ -324,7 +324,7 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {}
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {}
 
 
     void frz(realtype *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
@@ -338,7 +338,7 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
 
@@ -408,7 +408,7 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE {
     void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {}
 
 
-    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k) override { return {}; }
+    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k, const realtype *w) override { return {}; }
 
 
     std::string get_name() const override {
@@ -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 "40190b46b1b398e321314ded4169fe910b37c484";
+        return "3fb84cd5df12639f17b179d681e8ba4b5be8a160";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_nested_events_py/deltaqB.cpp b/models/model_nested_events_py/deltaqB.cpp
index 147390e017..b30e9fcfcb 100644
--- a/models/model_nested_events_py/deltaqB.cpp
+++ b/models/model_nested_events_py/deltaqB.cpp
@@ -5,6 +5,7 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -13,7 +14,7 @@
 namespace amici {
 namespace model_model_nested_events_py {
 
-void deltaqB_model_nested_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
+void deltaqB_model_nested_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
     switch(ie) {
         case 2:
             switch(ip) {
diff --git a/models/model_nested_events_py/deltaxB.cpp b/models/model_nested_events_py/deltaxB.cpp
index e4d155859f..1cdf5b5eb6 100644
--- a/models/model_nested_events_py/deltaxB.cpp
+++ b/models/model_nested_events_py/deltaxB.cpp
@@ -5,6 +5,7 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -13,7 +14,7 @@
 namespace amici {
 namespace model_model_nested_events_py {
 
-void deltaxB_model_nested_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl){
+void deltaxB_model_nested_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl){
     switch(ie) {
         case 0:
             deltaxB[0] = xB0*(dVirusdt - xdot_old0)/(Heaviside_1*Virus*rho_V - Virus*delta_V);
diff --git a/models/model_nested_events_py/explicit_roots.cpp b/models/model_nested_events_py/explicit_roots.cpp
index 8bbdb8c8d9..14201f899c 100644
--- a/models/model_nested_events_py/explicit_roots.cpp
+++ b/models/model_nested_events_py/explicit_roots.cpp
@@ -8,7 +8,7 @@
 namespace amici {
 namespace model_model_nested_events_py {
 
-std::vector<std::vector<realtype>> explicit_roots_model_nested_events_py(const realtype *p, const realtype *k){
+std::vector<std::vector<realtype>> explicit_roots_model_nested_events_py(const realtype *p, const realtype *k, const realtype *w){
     return {
         {t_0}
     };
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 120258bd7c..53ce2a747f 100644
--- a/models/model_nested_events_py/model_nested_events_py.h
+++ b/models/model_nested_events_py/model_nested_events_py.h
@@ -38,7 +38,7 @@ extern void dJydy_rowvals_model_nested_events_py(SUNMatrixWrapper &rowvals, int
 
 
 
-extern void root_model_nested_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
+extern void root_model_nested_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 
 
 
@@ -77,11 +77,11 @@ extern void xdot_model_nested_events_py(realtype *xdot, const realtype t, const
 extern void y_model_nested_events_py(realtype *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w);
 
 
-extern void stau_model_nested_events_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
+extern void stau_model_nested_events_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
 extern void deltax_model_nested_events_py(double *deltax, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old);
 extern void deltasx_model_nested_events_py(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old);
-extern void deltaxB_model_nested_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl);
-extern void deltaqB_model_nested_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
+extern void deltaxB_model_nested_events_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl);
+extern void deltaqB_model_nested_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
 
 extern void x_solver_model_nested_events_py(realtype *x_solver, const realtype *x_rdata);
 
@@ -99,7 +99,7 @@ extern void x_solver_model_nested_events_py(realtype *x_solver, const realtype *
 extern std::vector<HermiteSpline> create_splines_model_nested_events_py(const realtype *p, const realtype *k);
 
 
-extern std::vector<std::vector<realtype>> explicit_roots_model_nested_events_py(const realtype *p, const realtype *k);
+extern std::vector<std::vector<realtype>> explicit_roots_model_nested_events_py(const realtype *p, const realtype *k, const realtype *w);
 /**
  * @brief AMICI-generated model subclass.
  */
@@ -203,13 +203,13 @@ class Model_model_nested_events_py : public amici::Model_ODE {
     }
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {
-        deltaxB_model_nested_events_py(deltaxB, t, x, p, k, h, dx, ie, xdot, xdot_old, x_old, xB, tcl);
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {
+        deltaxB_model_nested_events_py(deltaxB, t, x, p, k, h, w, dx, ie, xdot, xdot_old, x_old, xB, tcl);
     }
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
-        deltaqB_model_nested_events_py(deltaqB, t, x, p, k, h, dx, ip, ie, xdot, xdot_old, x_old, xB);
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
+        deltaqB_model_nested_events_py(deltaqB, t, x, p, k, h, w, dx, ip, ie, xdot, xdot_old, x_old, xB);
     }
 
 
@@ -318,8 +318,8 @@ class Model_model_nested_events_py : public amici::Model_ODE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {
-        root_model_nested_events_py(root, t, x, p, k, h, tcl);
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {
+        root_model_nested_events_py(root, t, x, p, k, h, w, tcl);
     }
 
 
@@ -334,8 +334,8 @@ class Model_model_nested_events_py : public amici::Model_ODE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
-        stau_model_nested_events_py(stau, t, x, p, k, h, dx, tcl, sx, ip, ie);
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
+        stau_model_nested_events_py(stau, t, x, p, k, h, w, dx, tcl, sx, ip, ie);
     }
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {
@@ -406,8 +406,8 @@ class Model_model_nested_events_py : public amici::Model_ODE {
     void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {}
 
 
-    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k) override {
-        return explicit_roots_model_nested_events_py(p, k);
+    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k, const realtype *w) override {
+        return explicit_roots_model_nested_events_py(p, k, w);
     }
 
 
@@ -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 "40190b46b1b398e321314ded4169fe910b37c484";
+        return "3fb84cd5df12639f17b179d681e8ba4b5be8a160";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_nested_events_py/root.cpp b/models/model_nested_events_py/root.cpp
index a5121b379b..0b9186b9dc 100644
--- a/models/model_nested_events_py/root.cpp
+++ b/models/model_nested_events_py/root.cpp
@@ -9,7 +9,7 @@
 namespace amici {
 namespace model_model_nested_events_py {
 
-void root_model_nested_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){
+void root_model_nested_events_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl){
     root[0] = Virus - 1;
     root[1] = 1 - Virus;
     root[2] = t - t_0;
diff --git a/models/model_nested_events_py/stau.cpp b/models/model_nested_events_py/stau.cpp
index 0c26a548e3..da044ed1ad 100644
--- a/models/model_nested_events_py/stau.cpp
+++ b/models/model_nested_events_py/stau.cpp
@@ -10,7 +10,7 @@
 namespace amici {
 namespace model_model_nested_events_py {
 
-void stau_model_nested_events_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
+void stau_model_nested_events_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
     switch(ie) {
         case 0:
             switch(ip) {
diff --git a/models/model_neuron_py/deltaqB.cpp b/models/model_neuron_py/deltaqB.cpp
index aa6a33f49b..60d5fb86b0 100644
--- a/models/model_neuron_py/deltaqB.cpp
+++ b/models/model_neuron_py/deltaqB.cpp
@@ -6,6 +6,7 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -14,7 +15,7 @@
 namespace amici {
 namespace model_model_neuron_py {
 
-void deltaqB_model_neuron_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
+void deltaqB_model_neuron_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
     switch(ie) {
         case 0:
             switch(ip) {
diff --git a/models/model_neuron_py/deltaxB.cpp b/models/model_neuron_py/deltaxB.cpp
index 4122018dfd..f12934bed8 100644
--- a/models/model_neuron_py/deltaxB.cpp
+++ b/models/model_neuron_py/deltaxB.cpp
@@ -6,6 +6,7 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -14,7 +15,7 @@
 namespace amici {
 namespace model_model_neuron_py {
 
-void deltaxB_model_neuron_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl){
+void deltaxB_model_neuron_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl){
     switch(ie) {
         case 0:
             deltaxB[0] = xB0*(xdot_old0/(I0 - u + (1.0/25.0)*std::pow(v, 2) + 5*v + 140) + (dvdt - xdot_old0)/(I0 - u + (1.0/25.0)*std::pow(v, 2) + 5*v + 140) - 1) + xB1*(dudt - xdot_old1)/(I0 - u + (1.0/25.0)*std::pow(v, 2) + 5*v + 140);
diff --git a/models/model_neuron_py/model_neuron_py.h b/models/model_neuron_py/model_neuron_py.h
index d300ce9b03..57aa2aec2c 100644
--- a/models/model_neuron_py/model_neuron_py.h
+++ b/models/model_neuron_py/model_neuron_py.h
@@ -38,7 +38,7 @@ extern void dJzdz_model_neuron_py(realtype *dJzdz, const int iz, const realtype
 extern void Jrz_model_neuron_py(realtype *Jrz, const int iz, const realtype *p, const realtype *k, const realtype *rz, const realtype *sigmaz);
 extern void dJrzdsigma_model_neuron_py(realtype *dJrzdsigma, const int iz, const realtype *p, const realtype *k, const realtype *rz, const realtype *sigmaz);
 extern void dJrzdz_model_neuron_py(realtype *dJrzdz, const int iz, const realtype *p, const realtype *k, const realtype *rz, const realtype *sigmaz);
-extern void root_model_neuron_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
+extern void root_model_neuron_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 
 
 
@@ -77,11 +77,11 @@ extern void xdot_model_neuron_py(realtype *xdot, const realtype t, const realtyp
 extern void y_model_neuron_py(realtype *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w);
 extern void z_model_neuron_py(realtype *z, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h);
 extern void rz_model_neuron_py(realtype *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h);
-extern void stau_model_neuron_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
+extern void stau_model_neuron_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
 extern void deltax_model_neuron_py(double *deltax, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old);
 extern void deltasx_model_neuron_py(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old);
-extern void deltaxB_model_neuron_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl);
-extern void deltaqB_model_neuron_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
+extern void deltaxB_model_neuron_py(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl);
+extern void deltaqB_model_neuron_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
 
 extern void x_solver_model_neuron_py(realtype *x_solver, const realtype *x_rdata);
 
@@ -213,13 +213,13 @@ class Model_model_neuron_py : public amici::Model_ODE {
     }
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {
-        deltaxB_model_neuron_py(deltaxB, t, x, p, k, h, dx, ie, xdot, xdot_old, x_old, xB, tcl);
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {
+        deltaxB_model_neuron_py(deltaxB, t, x, p, k, h, w, dx, ie, xdot, xdot_old, x_old, xB, tcl);
     }
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
-        deltaqB_model_neuron_py(deltaqB, t, x, p, k, h, dx, ip, ie, xdot, xdot_old, x_old, xB);
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
+        deltaqB_model_neuron_py(deltaqB, t, x, p, k, h, w, dx, ip, ie, xdot, xdot_old, x_old, xB);
     }
 
 
@@ -332,8 +332,8 @@ class Model_model_neuron_py : public amici::Model_ODE {
     }
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {
-        root_model_neuron_py(root, t, x, p, k, h, tcl);
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {
+        root_model_neuron_py(root, t, x, p, k, h, w, tcl);
     }
 
 
@@ -352,8 +352,8 @@ class Model_model_neuron_py : public amici::Model_ODE {
     }
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
-        stau_model_neuron_py(stau, t, x, p, k, h, dx, tcl, sx, ip, ie);
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
+        stau_model_neuron_py(stau, t, x, p, k, h, w, dx, tcl, sx, ip, ie);
     }
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {
@@ -430,7 +430,7 @@ class Model_model_neuron_py : public amici::Model_ODE {
     void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {}
 
 
-    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k) override { return {}; }
+    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k, const realtype *w) override { return {}; }
 
 
     std::string get_name() const override {
@@ -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 "40190b46b1b398e321314ded4169fe910b37c484";
+        return "3fb84cd5df12639f17b179d681e8ba4b5be8a160";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_neuron_py/root.cpp b/models/model_neuron_py/root.cpp
index 4cf154c42d..9710c9f5a6 100644
--- a/models/model_neuron_py/root.cpp
+++ b/models/model_neuron_py/root.cpp
@@ -10,7 +10,7 @@
 namespace amici {
 namespace model_model_neuron_py {
 
-void root_model_neuron_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){
+void root_model_neuron_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl){
     root[0] = v - 30;
 }
 
diff --git a/models/model_neuron_py/stau.cpp b/models/model_neuron_py/stau.cpp
index 31b38f5e67..26ccb3a75d 100644
--- a/models/model_neuron_py/stau.cpp
+++ b/models/model_neuron_py/stau.cpp
@@ -11,7 +11,7 @@
 namespace amici {
 namespace model_model_neuron_py {
 
-void stau_model_neuron_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
+void stau_model_neuron_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
     switch(ie) {
         case 0:
             switch(ip) {
diff --git a/models/model_robertson_py/model_robertson_py.h b/models/model_robertson_py/model_robertson_py.h
index 85caa538e9..b3aef8f8d1 100644
--- a/models/model_robertson_py/model_robertson_py.h
+++ b/models/model_robertson_py/model_robertson_py.h
@@ -195,10 +195,10 @@ class Model_model_robertson_py : public amici::Model_DAE {
     void fdeltasx(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old) override {}
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
 
 
     void fdrzdp(realtype *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {}
@@ -306,7 +306,7 @@ class Model_model_robertson_py : public amici::Model_DAE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {}
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {}
 
 
     void frz(realtype *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
@@ -320,7 +320,7 @@ class Model_model_robertson_py : public amici::Model_DAE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
 
@@ -392,7 +392,7 @@ class Model_model_robertson_py : public amici::Model_DAE {
     void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {}
 
 
-    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k) override { return {}; }
+    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k, const realtype *w) override { return {}; }
 
 
     std::string get_name() const override {
@@ -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 "40190b46b1b398e321314ded4169fe910b37c484";
+        return "3fb84cd5df12639f17b179d681e8ba4b5be8a160";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/models/model_steadystate_py/model_steadystate_py.h b/models/model_steadystate_py/model_steadystate_py.h
index 8393c6ddf7..4d8238872d 100644
--- a/models/model_steadystate_py/model_steadystate_py.h
+++ b/models/model_steadystate_py/model_steadystate_py.h
@@ -195,10 +195,10 @@ class Model_model_steadystate_py : public amici::Model_ODE {
     void fdeltasx(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old) override {}
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
 
 
     void fdrzdp(realtype *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {}
@@ -306,7 +306,7 @@ class Model_model_steadystate_py : public amici::Model_ODE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {}
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {}
 
 
     void frz(realtype *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
@@ -320,7 +320,7 @@ class Model_model_steadystate_py : public amici::Model_ODE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
 
@@ -392,7 +392,7 @@ class Model_model_steadystate_py : public amici::Model_ODE {
     void fdtotal_cldx_rdata_rowvals(SUNMatrixWrapper &rowvals) override {}
 
 
-    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k) override { return {}; }
+    std::vector<std::vector<realtype>> fexplicit_roots(const realtype *p, const realtype *k, const realtype *w) override { return {}; }
 
 
     std::string get_name() const override {
@@ -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 "40190b46b1b398e321314ded4169fe910b37c484";
+        return "3fb84cd5df12639f17b179d681e8ba4b5be8a160";
     }
 
     bool has_quadratic_llh() const override {
diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py
index 3807f96552..31a9650938 100644
--- a/python/sdist/amici/de_model.py
+++ b/python/sdist/amici/de_model.py
@@ -7,7 +7,9 @@
 import logging
 import re
 from collections.abc import Callable, Sequence
+from contextlib import suppress
 from itertools import chain
+from operator import itemgetter
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -48,7 +50,6 @@
     ObservableTransformation,
     _default_simplify,
     amici_time_symbol,
-    smart_subs_dict,
     toposort_symbols,
     unique_preserve_order,
 )
@@ -646,14 +647,20 @@ def num_events(self) -> int:
 
     def num_events_solver(self) -> int:
         """
-        Number of Events.
+        Number of Events that rely on numerical root-finding.
 
         :return:
             number of event symbols (length of the root vector in AMICI)
         """
-        constant_syms = set(self.sym("k")) | set(self.sym("p"))
+        # TODO(performance): we could include constant `x` here as well
+        #   (dx/dt = 0 AND not target of event assignments)
+        #   this will require passing `x` to `fexplicit_roots`
+        # TODO(performance): cache result, and don't repeat list symbols
+        # Note that`self.num_events_solver` is currently
+        #  NOT the same as len(self.get_implicit_roots())
+        static_syms = self._static_symbols(["k", "p", "w"])
         return sum(
-            not event.has_explicit_trigger_times(constant_syms)
+            not event.has_explicit_trigger_times(static_syms)
             for event in self.events()
         )
 
@@ -938,6 +945,29 @@ def static_indices(self, name: str) -> list[int]:
 
         raise NotImplementedError(name)
 
+    def _static_symbols(self, names: list[str]) -> set[sp.Symbol]:
+        """
+        Return the static symbols among the given model entities.
+
+        E.g., `static_symbols(["p", "w"])` returns all symbols in `p` and `w`
+        that do not depend on time, neither directly nor indirectly.
+        """
+        result = set()
+
+        for name in names:
+            if name in ("k", "p"):
+                result |= set(self.sym(name))
+            elif name == "w":
+                static_indices = self.static_indices("w")
+                if len(static_indices) == 1:
+                    result.add(self.sym("w")[static_indices[0]])
+                elif len(static_indices) > 1:
+                    result |= set(itemgetter(*static_indices)(self.sym("w")))
+            else:
+                raise ValueError(name)
+
+        return result
+
     def dynamic_indices(self, name: str) -> list[int]:
         """
         Return the indices of dynamic expressions in the given model entity.
@@ -1114,6 +1144,7 @@ def generate_basic_variables(self) -> None:
         Generates the symbolic identifiers for all variables in
         ``DEModel._variable_prototype``
         """
+        self.parse_events()
         self._reorder_events()
 
         for var in self._variable_prototype:
@@ -1132,36 +1163,13 @@ def parse_events(self) -> None:
         and replaces the formulae of the found roots by identifiers of AMICI's
         Heaviside function implementation in the right-hand side
         """
-        # toposorted w_sym -> w_expr for substitution of 'w' in trigger function
-        #  do only once. `w` is not modified during this function.
-        w_toposorted = toposort_symbols(
-            dict(
-                zip(
-                    [expr.get_sym() for expr in self._expressions],
-                    [expr.get_val() for expr in self._expressions],
-                    strict=True,
-                )
-            )
-        )
-
         # Track all roots functions in the right-hand side
         roots = copy.deepcopy(self._events)
         for state in self._differential_states:
-            state.set_dt(
-                self._process_heavisides(state.get_dt(), roots, w_toposorted)
-            )
+            state.set_dt(self._process_heavisides(state.get_dt(), roots))
 
         for expr in self._expressions:
-            expr.set_val(
-                self._process_heavisides(expr.get_val(), roots, w_toposorted)
-            )
-
-        # remove all possible Heavisides from roots, which may arise from
-        # the substitution of `'w'` in `_collect_heaviside_roots`
-        for root in roots:
-            root.set_val(
-                self._process_heavisides(root.get_val(), roots, w_toposorted)
-            )
+            expr.set_val(self._process_heavisides(expr.get_val(), roots))
 
         # Now add the found roots to the model components
         for root in roots:
@@ -1171,30 +1179,65 @@ def parse_events(self) -> None:
             # add roots of heaviside functions
             self.add_component(root)
 
-        # Substitute 'w' expressions into root expressions, to avoid rewriting
-        # 'root.cpp' and 'stau.cpp' headers to include 'w.h'.
-        for event in self.events():
-            event.set_val(event.get_val().subs(w_toposorted))
-
     def _reorder_events(self) -> None:
         """
         Re-order events - first those that require root tracking,
         then the others.
         """
-        constant_syms = set(self.sym("k")) | set(self.sym("p"))
+        # Currently, the C++ simulations relies on the order of events:
+        #  those that require numerical root-finding must come first, then
+        #  those with explicit trigger times that don't depend on dynamic
+        #  variables.
+        #  TODO: This re-ordering here is a bit ugly, because we already need
+        #   to generate certain model equations to perform this ordering.
+        #   Ideally, we'd split froot into explicit and implicit parts during
+        #   code generation instead (as already done for jax models).
+        if not self.events():
+            return
+
+        static_syms = self._static_symbols(["k", "p", "w"])
+
+        # ensure that we don't have computed any root-related symbols/equations
+        #  yet, because the re-ordering might invalidate them
+        # check after `self._static_symbols` which itself generates certain
+        #  equations
+        if (
+            generated := set(self._syms)
+            | set(self._eqs)
+            | set(self._sparsesyms)
+            | set(self._sparseeqs)
+        ) and (
+            "root" in generated
+            or any(
+                name.startswith("droot") or name.endswith("droot")
+                for name in generated
+            )
+        ):
+            raise AssertionError(
+                "This function must be called before computing any "
+                "root-related symbols/equations. "
+                "The following symbols/equations are already "
+                f"generated: {generated}"
+            )
+
         self._events = list(
             chain(
                 itertools.filterfalse(
-                    lambda e: e.has_explicit_trigger_times(constant_syms),
+                    lambda e: e.has_explicit_trigger_times(static_syms),
                     self._events,
                 ),
                 filter(
-                    lambda e: e.has_explicit_trigger_times(constant_syms),
+                    lambda e: e.has_explicit_trigger_times(static_syms),
                     self._events,
                 ),
             )
         )
 
+        # regenerate after re-ordering
+        with suppress(KeyError):
+            del self._syms["h"]
+        self.sym("h")
+
     def get_appearance_counts(self, idxs: list[int]) -> list[int]:
         """
         Counts how often a state appears in the time derivative of
@@ -1503,18 +1546,26 @@ def _compute_equation(self, name: str) -> None:
             self._eqs[name] = smart_jacobian(self.eq("root"), time_symbol)
 
         elif name == "drootdt_total":
+            # root(t, x(t), w(t, x(t)))
+            # drootdt_total = drootdt + drootdx * dxdt + drootdw * dwdt_total
+            # dwdt_total = dwdt + dwdx * dxdt
             self._eqs[name] = self.eq("drootdt")
-            # backsubstitution of optimized right-hand side terms into RHS
-            # calling subs() is costly. We can skip it if we don't have any
-            # state-dependent roots.
+            xdot = self.eq("xdot")
+
+            # += drootdx * dxdt
             if self.num_states_solver() and not smart_is_zero_matrix(
-                drootdx := self.eq("drootdx")
-            ):
-                w_sorted = toposort_symbols(
-                    dict(zip(self.sym("w"), self.eq("w"), strict=True))
+                drootdx_explicit := smart_jacobian(
+                    self.eq("root"), self.sym("x")
                 )
-                tmp_xdot = smart_subs_dict(self.eq("xdot"), w_sorted)
-                self._eqs[name] += smart_multiply(drootdx, tmp_xdot)
+            ):
+                self._eqs[name] += smart_multiply(drootdx_explicit, xdot)
+
+            # += drootdw * dwdt_total
+            if not smart_is_zero_matrix(drootdw := self.eq("drootdw")):
+                dwdt = self.eq("dwdt")
+                dwdx_explicit = smart_jacobian(self.eq("w"), self.sym("x"))
+                dwdt_total = dwdt + smart_multiply(dwdx_explicit, xdot)
+                self._eqs[name] += smart_multiply(drootdw, dwdt_total)
 
         elif name == "deltax":
             # fill boluses for Heaviside functions, as empty state updates
@@ -1598,18 +1649,36 @@ def _compute_equation(self, name: str) -> None:
             ]
 
         elif name == "dtaudx":
+            drootdx_explicit = smart_jacobian(self.eq("root"), self.sym("x"))
+            drootdt_total = self.eq("drootdt_total")
+            drootdw = self.eq("drootdw")
+            dwdx = self.eq("dwdx")
+
             self._eqs[name] = [
-                self.eq("drootdx")[ie, :] / self.eq("drootdt_total")[ie]
-                if not self.eq("drootdt_total")[ie].is_zero
-                else sp.zeros(*self.eq("drootdx")[ie, :].shape)
+                (
+                    # drootdx + drootdw * dwdx
+                    drootdx_explicit[ie, :] + drootdw[ie, :] * dwdx
+                )
+                / drootdt_total[ie]
+                if not drootdt_total[ie].is_zero
+                else sp.zeros(*drootdx_explicit[ie, :].shape)
                 for ie in range(self.num_events())
             ]
 
         elif name == "dtaudp":
+            drootdp_explicit = smart_jacobian(self.eq("root"), self.sym("p"))
+            drootdt_total = self.eq("drootdt_total")
+            drootdw = self.eq("drootdw")
+            drootdp = self.eq("drootdp")
+            dwdp = self.eq("dwdp")
             self._eqs[name] = [
-                self.eq("drootdp")[ie, :] / self.eq("drootdt_total")[ie]
-                if not self.eq("drootdt_total")[ie].is_zero
-                else sp.zeros(*self.eq("drootdp")[ie, :].shape)
+                (
+                    # drootdp + drootdw * dwdp
+                    drootdp_explicit[ie, :] + drootdw[ie, :] * dwdp
+                )
+                / drootdt_total[ie]
+                if not drootdt_total[ie].is_zero
+                else sp.zeros(*drootdp[ie, :].shape)
                 for ie in range(self.num_events())
             ]
 
@@ -1757,6 +1826,9 @@ def _compute_equation(self, name: str) -> None:
                 smart_jacobian(self.eq("w")[self.num_cons_law() :, :], x)
             )
 
+        elif name == "dwdt":
+            self._eqs[name] = smart_jacobian(self.eq("w"), time_symbol)
+
         elif name == "iroot":
             self._eqs[name] = sp.Matrix(
                 [
@@ -1896,6 +1968,7 @@ def _derivative(self, eq: str, var: str, name: str = None) -> None:
             ("xdot", "x"): "w",  # has generic implementation in c++ code
             ("w", "w"): "tcl",  # dtcldw = 0
             ("w", "x"): "tcl",  # dtcldx = 0
+            ("root", "w"): "tcl",  # dtcldw = 0
         }
         # automatically detect chainrule
         chainvars = [
@@ -1945,7 +2018,9 @@ def _derivative(self, eq: str, var: str, name: str = None) -> None:
                             "attach this model."
                         )
 
-        if name == "dydw" and not smart_is_zero_matrix(derivative):
+        elif name in ("dydw", "drootdw") and not smart_is_zero_matrix(
+            derivative
+        ):
             dwdw = self.eq("dwdw")
             # h(k) = d{eq}dw*dwdw^k* (k=1)
             h = smart_multiply(derivative, dwdw)
@@ -2281,9 +2356,6 @@ def _get_unique_root(
             unique identifier for root, or ``None`` if the root is not
             time-dependent
         """
-        if not self._expr_is_time_dependent(root_found):
-            return None
-
         for root in roots:
             if (difference := (root_found - root.get_val())).is_zero or (
                 self._simplify and self._simplify(difference).is_zero
@@ -2301,6 +2373,12 @@ def _get_unique_root(
                 use_values_from_trigger_time=True,
             )
         )
+
+        if not self._expr_is_time_dependent(root_found):
+            # Not time-dependent. Return None, but we still need to create
+            # the event above
+            return None
+
         return roots[-1].get_sym()
 
     def _collect_heaviside_roots(
@@ -2331,7 +2409,6 @@ def _process_heavisides(
         self,
         dxdt: sp.Expr,
         roots: list[Event],
-        w_toposorted: dict[sp.Symbol, sp.Expr],
     ) -> sp.Expr:
         """
         Parses the RHS of a state variable, checks for Heaviside functions,
@@ -2343,8 +2420,6 @@ def _process_heavisides(
             right-hand side of state variable
         :param roots:
             list of known root functions with identifier
-        :param w_toposorted:
-            `w` symbols->expressions sorted in topological order
         :returns:
             dxdt with Heaviside functions replaced by amici helper variables
         """
@@ -2353,15 +2428,6 @@ def _process_heavisides(
         heavisides = []
         # run through the expression tree and get the roots
         tmp_roots_old = self._collect_heaviside_roots((dxdt,))
-        # substitute 'w' symbols in the root expression by their equations,
-        #  because currently,
-        #  1) root functions must not depend on 'w'
-        #  2) the check for time-dependence currently assumes only state
-        #     variables are implicitly time-dependent
-        tmp_roots_old = [
-            (a.subs(w_toposorted), b.subs(w_toposorted))
-            for a, b in tmp_roots_old
-        ]
         for tmp_root_old, tmp_x0_old in unique_preserve_order(tmp_roots_old):
             # we want unique identifiers for the roots
             tmp_root_new = self._get_unique_root(tmp_root_old, roots)
diff --git a/python/sdist/amici/exporters/sundials/cxx_functions.py b/python/sdist/amici/exporters/sundials/cxx_functions.py
index ae9eb7013d..df266fbb18 100644
--- a/python/sdist/amici/exporters/sundials/cxx_functions.py
+++ b/python/sdist/amici/exporters/sundials/cxx_functions.py
@@ -144,7 +144,7 @@ def var_in_signature(self, varname: str, ode: bool = True) -> bool:
     "root": _FunctionInfo(
         "realtype *root, const realtype t, const realtype *x, "
         "const realtype *p, const realtype *k, const realtype *h, "
-        "const realtype *tcl"
+        "const realtype *w, const realtype *tcl"
     ),
     "dwdp": _FunctionInfo(
         "realtype *dwdp, const realtype t, const realtype *x, "
@@ -288,7 +288,7 @@ def var_in_signature(self, varname: str, ode: bool = True) -> bool:
     "stau": _FunctionInfo(
         "realtype *stau, const realtype t, const realtype *x, "
         "const realtype *p, const realtype *k, const realtype *h, "
-        "const realtype *dx, "
+        "const realtype *w, const realtype *dx, "
         "const realtype *tcl, const realtype *sx, const int ip, "
         "const int ie"
     ),
@@ -313,7 +313,7 @@ def var_in_signature(self, varname: str, ode: bool = True) -> bool:
     "deltaxB": _FunctionInfo(
         "realtype *deltaxB, const realtype t, const realtype *x, "
         "const realtype *p, const realtype *k, const realtype *h, "
-        "const realtype *dx, "
+        "const realtype *w, const realtype *dx, "
         "const int ie, const realtype *xdot, const realtype *xdot_old, "
         "const realtype *x_old, "
         "const realtype *xB, const realtype *tcl"
@@ -321,7 +321,7 @@ def var_in_signature(self, varname: str, ode: bool = True) -> bool:
     "deltaqB": _FunctionInfo(
         "realtype *deltaqB, const realtype t, const realtype *x, "
         "const realtype *p, const realtype *k, const realtype *h, "
-        "const realtype *dx, "
+        "const realtype *w, const realtype *dx, "
         "const int ip, const int ie, const realtype *xdot, "
         "const realtype *xdot_old, const realtype *x_old, const realtype *xB"
     ),
@@ -408,7 +408,7 @@ def var_in_signature(self, varname: str, ode: bool = True) -> bool:
         "const realtype *p, const realtype *k, const realtype *h"
     ),
     "explicit_roots": _FunctionInfo(
-        "const realtype *p, const realtype *k",
+        "const realtype *p, const realtype *k, const realtype *w",
         return_type="std::vector<std::vector<realtype>>",
         default_return_value="{}",
         header=[
diff --git a/python/sdist/amici/exporters/sundials/de_export.py b/python/sdist/amici/exporters/sundials/de_export.py
index e96e2fe95b..06f5ecca68 100644
--- a/python/sdist/amici/exporters/sundials/de_export.py
+++ b/python/sdist/amici/exporters/sundials/de_export.py
@@ -914,14 +914,14 @@ def _get_create_splines_body(self):
     def _get_explicit_roots_body(self) -> list[str]:
         events = self.model.events()
         lines = []
-        constant_syms = set(self.model.sym("k")) | set(self.model.sym("p"))
+        static_syms = self.model._static_symbols(["k", "p", "w"])
 
         for event_idx, event in enumerate(events):
             if not (
                 tigger_times := {
                     tt
                     for tt in event.get_trigger_times()
-                    if tt.free_symbols.issubset(constant_syms)
+                    if tt.free_symbols.issubset(static_syms)
                 }
             ):
                 continue
diff --git a/python/sdist/amici/importers/pysb/__init__.py b/python/sdist/amici/importers/pysb/__init__.py
index 7b3f9ce749..235ae8b72c 100644
--- a/python/sdist/amici/importers/pysb/__init__.py
+++ b/python/sdist/amici/importers/pysb/__init__.py
@@ -398,7 +398,6 @@ def ode_model_from_pysb_importer(
 
     _process_stoichiometric_matrix(model, ode, fixed_parameters)
 
-    ode.parse_events()
     ode.generate_basic_variables()
 
     return ode
diff --git a/python/sdist/amici/importers/sbml/__init__.py b/python/sdist/amici/importers/sbml/__init__.py
index 37e2776f8a..38c7d29021 100644
--- a/python/sdist/amici/importers/sbml/__init__.py
+++ b/python/sdist/amici/importers/sbml/__init__.py
@@ -700,9 +700,7 @@ def _build_ode_model(
         if hybridization:
             ode_model._process_hybridization(hybridization)
 
-        ode_model.parse_events()
         # substitute SBML-rateOf constructs
-        #  must be done after parse_events, but before generate_basic_variables
         self._process_sbml_rate_of(ode_model)
 
         # fill in 'self._sym' based on prototypes and components in ode_model
@@ -3142,14 +3140,6 @@ def do_subs(expr, rate_ofs) -> sp.Expr:
                 component.set_val(do_subs(component.get_val(), rate_ofs))
 
             if isinstance(component, Event):
-                if rate_ofs:
-                    # currently, `root` cannot depend on `w`.
-                    #  this could be changed, but for now,
-                    #  we just flatten out w expressions
-                    component.set_val(
-                        smart_subs_dict(component.get_val(), w_toposorted)
-                    )
-
                 # for events, also substitute in state updates
                 for target, assignment in component._assignments.items():
                     if rate_ofs := assignment.find(rate_of_func):
diff --git a/python/tests/test_events.py b/python/tests/test_events.py
index d0558b6fd1..655613d4f9 100644
--- a/python/tests/test_events.py
+++ b/python/tests/test_events.py
@@ -1294,3 +1294,59 @@ def test_gh2926(tempdir):
     rdata = amici.run_simulation(model, solver)
     assert rdata.status == amici.AMICI_SUCCESS
     assert rdata.by_id("x1").tolist() == [1.0, 1.0, 2.0, 2.0]
+
+
+@skip_on_valgrind
+def test_event_with_w_dependent_trigger(tempdir):
+    """Test sensitivities for events with trigger depending on
+    cascading expressions in `w`."""
+
+    model_name = "test_event_with_w_dependent_trigger"
+    model = antimony2amici(
+        r"""
+        one = 1
+        two = 2
+        a := two^2 - 2 # 2
+        b := a^2 - 1  # 3
+        c := b * 2 + x / 10  # 6 + time / 10
+        d := c + a + (a - 1) * time / 10 # -> d = 8 + time / 5
+
+        x = 0
+        x' = a / 2  # x' = 1
+        target = 0
+
+        E1: at x >= d:  # triggers at time = 8 + time / 5 <=> time = 10
+            target = x + one;
+        """,
+        model_name=model_name,
+        output_dir=tempdir,
+    )
+
+    model.set_timepoints([0, 5, 9, 11])
+    solver = model.create_solver()
+    solver.set_sensitivity_order(SensitivityOrder.first)
+    solver.set_sensitivity_method(SensitivityMethod.forward)
+    solver.set_relative_tolerance(1e-14)
+
+    # generate synthetic measurements
+    rdata = amici.run_simulation(model, solver)
+    assert rdata.status == amici.AMICI_SUCCESS
+    # check that event triggered correctly
+    assert np.isclose(rdata.by_id("target")[-1], 11.0)
+    edata = amici.ExpData(rdata, 1, 0)
+
+    # check sensitivities against finite differences
+
+    for sens_method in (
+        SensitivityMethod.forward,
+        SensitivityMethod.adjoint,
+    ):
+        solver.set_sensitivity_method(sens_method)
+        check_derivatives(
+            model,
+            solver,
+            edata=edata,
+            atol=1e-6,
+            rtol=1e-6,
+            epsilon=1e-8,
+        )
diff --git a/python/tests/test_pysb.py b/python/tests/test_pysb.py
index ef43983a6a..ed796ebe68 100644
--- a/python/tests/test_pysb.py
+++ b/python/tests/test_pysb.py
@@ -328,7 +328,7 @@ def test_names_and_ids(pysb_example_presimulation_module):
 
 
 @skip_on_valgrind
-def test_heavyside_and_special_symbols():
+def test_heaviside_and_special_symbols():
     pysb.SelfExporter.cleanup()  # reset pysb
     pysb.SelfExporter.do_export = True
 
diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py
index 574bfcf7b2..48d44341b8 100644
--- a/python/tests/test_sbml_import.py
+++ b/python/tests/test_sbml_import.py
@@ -1041,6 +1041,10 @@ def test_regression_2700(tempdir):
     model_name = "regression_2700"
     antimony2amici(
         """
+    # we need some differential state, otherwise import will currently fail
+    x = 0
+    x' = 0
+
     a = 1
     # condition is always true, so `pp` should be 1
     pp := piecewise(1, a >= 1 && a <= 1, 0)
diff --git a/src/abstract_model.cpp b/src/abstract_model.cpp
index 15b9b9acf1..929e9033db 100644
--- a/src/abstract_model.cpp
+++ b/src/abstract_model.cpp
@@ -59,8 +59,8 @@ void AbstractModel::fdx0(AmiVector& /*x0*/, AmiVector& /*dx0*/) {
 void AbstractModel::fstau(
     realtype* /*stau*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
-    realtype const* /*dx*/, realtype const* /*tcl*/, realtype const* /*sx*/,
-    int const /*ip*/, int const /*ie*/
+    realtype const* /*w*/, realtype const* /*dx*/, realtype const* /*tcl*/,
+    realtype const* /*sx*/, int const /*ip*/, int const /*ie*/
 ) {
     throw AmiException(
         "Requested functionality is not supported as %s is "
@@ -230,9 +230,9 @@ void AbstractModel::fdeltasx(
 void AbstractModel::fdeltaxB(
     realtype* /*deltaxB*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
-    realtype const* /*dx*/, int const /*ie*/, realtype const* /*xdot*/,
-    realtype const* /*xdot_old*/, realtype const* /*x_old*/,
-    realtype const* /*xB*/, realtype const* /*tcl*/
+    realtype const* /*w*/, realtype const* /*dx*/, int const /*ie*/,
+    realtype const* /*xdot*/, realtype const* /*xdot_old*/,
+    realtype const* /*x_old*/, realtype const* /*xB*/, realtype const* /*tcl*/
 ) {
     throw AmiException(
         "Requested functionality is not supported as %s is "
@@ -244,8 +244,8 @@ void AbstractModel::fdeltaxB(
 void AbstractModel::fdeltaqB(
     realtype* /*deltaqB*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
-    realtype const* /*dx*/, int const /*ip*/, int const /*ie*/,
-    realtype const* /*xdot*/, realtype const* /*xdot_old*/,
+    realtype const* /*w*/, realtype const* /*dx*/, int const /*ip*/,
+    int const /*ie*/, realtype const* /*xdot*/, realtype const* /*xdot_old*/,
     realtype const* /*x_old*/, realtype const* /*xB*/
 ) {
     throw AmiException(
diff --git a/src/model.cpp b/src/model.cpp
index f40d7d2987..b8844bf6fe 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -442,9 +442,12 @@ void Model::reinit_events(
 void Model::reinit_explicit_roots() {
     explicit_roots_.clear();
 
-    // evaluate timepoints
+    // Evaluate event trigger timepoints.
+    // This assumes that `fw` was called before.
+    // That happens during model initialization or pre-equilibration.
     auto const exp_roots = fexplicit_roots(
-        state_.unscaled_parameters.data(), state_.fixed_parameters.data()
+        state_.unscaled_parameters.data(), state_.fixed_parameters.data(),
+        derived_state_.w_.data()
     );
     Expects(exp_roots.size() == gsl::narrow<size_t>(ne - ne_solver));
 
@@ -1447,8 +1450,8 @@ void Model::get_event_time_sensitivity(
         fstau(
             &stau.at(ip), t, compute_x_pos(x),
             state_.unscaled_parameters.data(), state_.fixed_parameters.data(),
-            state_.h.data(), dx.data(), state_.total_cl.data(), sx.data(ip),
-            plist(ip), ie
+            state_.h.data(), derived_state_.w_.data(), dx.data(),
+            state_.total_cl.data(), sx.data(ip), plist(ip), ie
         );
     }
 }
@@ -1520,8 +1523,8 @@ void Model::add_adjoint_state_event_update(
     fdeltaxB(
         derived_state_.deltaxB_.data(), t, compute_x_pos(x),
         state_.unscaled_parameters.data(), state_.fixed_parameters.data(),
-        state_.h.data(), dx.data(), ie, xdot.data(), xdot_old.data(),
-        x_old.data(), xB.data(), state_.total_cl.data()
+        state_.h.data(), derived_state_.w_.data(), dx.data(), ie, xdot.data(),
+        xdot_old.data(), x_old.data(), xB.data(), state_.total_cl.data()
     );
 
     if (always_check_finite_) {
@@ -1546,7 +1549,8 @@ void Model::add_adjoint_quadrature_event_update(
         fdeltaqB(
             derived_state_.deltaqB_.data(), t, compute_x_pos(x),
             state_.unscaled_parameters.data(), state_.fixed_parameters.data(),
-            state_.h.data(), dx.data(), plist(ip), ie, xdot.data(),
+            state_.h.data(), derived_state_.w_.data(), dx.data(), plist(ip), ie,
+            xdot.data(),
             xdot_old.data(), x_old.data(), xB.data()
         );
 
diff --git a/src/model_dae.cpp b/src/model_dae.cpp
index 94b5a8f300..b747810a98 100644
--- a/src/model_dae.cpp
+++ b/src/model_dae.cpp
@@ -95,10 +95,13 @@ void Model_DAE::froot(
 ) {
     std::ranges::fill(root, 0.0);
     auto const x_pos = compute_x_pos(x);
+    // TODO(performance) only dynamic expressions? consider storing a flag
+    //  for whether froot actually depends on `w`.
+    fw(t, N_VGetArrayPointerConst(x_pos));
     froot(
         root.data(), t, N_VGetArrayPointerConst(x_pos),
         state_.unscaled_parameters.data(), state_.fixed_parameters.data(),
-        state_.h.data(), N_VGetArrayPointerConst(dx)
+        state_.h.data(), derived_state_.w_.data(), N_VGetArrayPointerConst(dx)
     );
 }
 
@@ -222,7 +225,7 @@ void Model_DAE::fJSparse(
 void Model_DAE::froot(
     realtype* /*root*/, realtype const /*t*/, realtype const* /*x*/,
     double const* /*p*/, double const* /*k*/, realtype const* /*h*/,
-    realtype const* /*dx*/
+    realtype const* /*w*/, realtype const* /*dx*/
 ) {
     throw AmiException(
         "Requested functionality is not supported as %s is not "
diff --git a/src/model_ode.cpp b/src/model_ode.cpp
index 46fde9a726..caa900c5eb 100644
--- a/src/model_ode.cpp
+++ b/src/model_ode.cpp
@@ -84,11 +84,14 @@ void Model_ODE::froot(
     realtype const t, const_N_Vector x, gsl::span<realtype> root
 ) {
     auto const x_pos = compute_x_pos(x);
+    // TODO(performance) only dynamic expressions? consider storing a flag
+    //  for whether froot actually depends on `w`.
+    fw(t, N_VGetArrayPointerConst(x_pos));
     std::ranges::fill(root, 0.0);
     froot(
         root.data(), t, N_VGetArrayPointerConst(x_pos),
         state_.unscaled_parameters.data(), state_.fixed_parameters.data(),
-        state_.h.data(), state_.total_cl.data()
+        state_.h.data(), derived_state_.w_.data(), state_.total_cl.data()
     );
 }
 
@@ -206,7 +209,7 @@ void Model_ODE::fJSparse_rowvals(SUNMatrixWrapper& /*JSparse*/) {
 void Model_ODE::froot(
     realtype* /*root*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
-    realtype const* /*tcl*/
+    realtype const* /*w*/, realtype const* /*tcl*/
 ) {
     throw AmiException(
         "Requested functionality is not supported as %s is not "