diff --git a/include/amici/exception.h b/include/amici/exception.h
index 7ee91d511b..431be48cb9 100644
--- a/include/amici/exception.h
+++ b/include/amici/exception.h
@@ -62,29 +62,35 @@ class AmiException : public std::exception {
 };
 
 /**
- * @brief cvode exception handler class
+ * @brief CVODE exception handler class
  */
 class CvodeException : public AmiException {
   public:
     /**
      * @brief Constructor
-     * @param error_code error code returned by cvode function
-     * @param function cvode function name
+     * @param error_code error code returned by CVODE function
+     * @param function CVODE function name
+     * @param extra Extra text to append to error message
      */
-    CvodeException(int error_code, char const* function);
+    CvodeException(
+        int error_code, char const* function, char const* extra = nullptr
+    );
 };
 
 /**
- * @brief ida exception handler class
+ * @brief IDA exception handler class
  */
 class IDAException : public AmiException {
   public:
     /**
      * @brief Constructor
-     * @param error_code error code returned by ida function
-     * @param function ida function name
+     * @param error_code error code returned by IDA function
+     * @param function IDA function name
+     * @param extra Extra text to append to error message
      */
-    IDAException(int error_code, char const* function);
+    IDAException(
+        int error_code, char const* function, char const* extra = nullptr
+    );
 };
 
 /**
diff --git a/include/amici/forwardproblem.h b/include/amici/forwardproblem.h
index 91e5d50791..8c500bf73c 100644
--- a/include/amici/forwardproblem.h
+++ b/include/amici/forwardproblem.h
@@ -7,7 +7,6 @@
 #include "amici/vector.h"
 #include <amici/amici.h>
 
-#include <memory>
 #include <sundials/sundials_direct.h>
 #include <vector>
 
@@ -197,7 +196,9 @@ class ForwardProblem {
     SimulationState const& getSimulationStateTimepoint(int it) const {
         if (model->getTimepoint(it) == initial_state_.t)
             return getInitialSimulationState();
-        return timepoint_states_.find(model->getTimepoint(it))->second;
+        auto map_iter = timepoint_states_.find(model->getTimepoint(it));
+        assert(map_iter != timepoint_states_.end());
+        return map_iter->second;
     };
 
     /**
@@ -273,9 +274,9 @@ class ForwardProblem {
     /**
      * @brief Execute everything necessary for the handling of data points
      *
-     * @param it index of data point
+     * @param t measurement timepoint
      */
-    void handleDataPoint(int it);
+    void handleDataPoint(realtype t);
 
     /**
      * @brief Applies the event bolus to the current state
@@ -368,7 +369,8 @@ class ForwardProblem {
      * @brief Array of flags indicating which root has been found.
      *
      * Array of length nr (ne) with the indices of the user functions gi found
-     * to have a root. For i = 0, . . . ,nr 1 if gi has a root, and = 0 if not.
+     * to have a root. For i = 0, . . . ,nr 1 or -1 if gi has a root, and = 0
+     * if not. See CVodeGetRootInfo for details.
      */
     std::vector<int> roots_found_;
 
diff --git a/include/amici/model.h b/include/amici/model.h
index 72f733e6cf..f5421bed3d 100644
--- a/include/amici/model.h
+++ b/include/amici/model.h
@@ -12,7 +12,6 @@
 #include "amici/vector.h"
 
 #include <map>
-#include <memory>
 #include <vector>
 
 namespace amici {
@@ -117,6 +116,8 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param ndxdotdp_explicit Number of nonzero elements in `dxdotdp_explicit`
      * @param ndxdotdx_explicit Number of nonzero elements in `dxdotdx_explicit`
      * @param w_recursion_depth Recursion depth of fw
+     * @param state_independent_events Map of events with state-independent
+     * triggers functions, mapping trigger timepoints to event indices.
      */
     Model(
         ModelDimensions const& model_dimensions,
@@ -124,7 +125,8 @@ class Model : public AbstractModel, public ModelDimensions {
         amici::SecondOrderMode o2mode, std::vector<amici::realtype> idlist,
         std::vector<int> z2event, bool pythonGenerated = false,
         int ndxdotdp_explicit = 0, int ndxdotdx_explicit = 0,
-        int w_recursion_depth = 0
+        int w_recursion_depth = 0,
+        std::map<realtype, std::vector<int>> state_independent_events = {}
     );
 
     /** Destructor. */
@@ -1449,6 +1451,15 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     SUNMatrixWrapper const& get_dxdotdp_full() const;
 
+    /**
+     * @brief Get trigger times for events that don't require root-finding.
+     *
+     * @return List of unique trigger points for events that don't require
+     * root-finding (i.e. that trigger at predetermined timepoints),
+     * in ascending order.
+     */
+    virtual std::vector<double> get_trigger_timepoints() const;
+
     /**
      * Flag indicating whether for
      * `amici::Solver::sensi_` == `amici::SensitivityOrder::second`
@@ -1462,6 +1473,12 @@ class Model : public AbstractModel, public ModelDimensions {
     /** Logger */
     Logger* logger = nullptr;
 
+    /**
+     * @brief Map of trigger timepoints to event indices for events that don't
+     * require root-finding.
+     */
+    std::map<realtype, std::vector<int>> state_independent_events_ = {};
+
   protected:
     /**
      * @brief Write part of a slice to a buffer according to indices specified
diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h
index 375ac018fd..b35cfa1d70 100644
--- a/include/amici/model_dae.h
+++ b/include/amici/model_dae.h
@@ -40,6 +40,8 @@ class Model_DAE : public Model {
      * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit
      * @param ndxdotdx_explicit number of nonzero elements dxdotdx_explicit
      * @param w_recursion_depth Recursion depth of fw
+     * @param state_independent_events Map of events with state-independent
+     * triggers functions, mapping trigger timepoints to event indices.
      */
     Model_DAE(
         ModelDimensions const& model_dimensions,
@@ -47,12 +49,13 @@ class Model_DAE : public Model {
         const SecondOrderMode o2mode, std::vector<realtype> const& idlist,
         std::vector<int> const& z2event, bool const pythonGenerated = false,
         int const ndxdotdp_explicit = 0, int const ndxdotdx_explicit = 0,
-        int const w_recursion_depth = 0
+        int const w_recursion_depth = 0,
+        std::map<realtype, std::vector<int>> state_independent_events = {}
     )
         : Model(
             model_dimensions, simulation_parameters, o2mode, idlist, z2event,
             pythonGenerated, ndxdotdp_explicit, ndxdotdx_explicit,
-            w_recursion_depth
+            w_recursion_depth, state_independent_events
         ) {
         derived_state_.M_ = SUNMatrixWrapper(nx_solver, nx_solver);
         auto M_nnz = static_cast<sunindextype>(
diff --git a/include/amici/model_dimensions.h b/include/amici/model_dimensions.h
index f0679dbe36..b5aa1ba21e 100644
--- a/include/amici/model_dimensions.h
+++ b/include/amici/model_dimensions.h
@@ -31,6 +31,7 @@ struct ModelDimensions {
      * @param nz Number of event observables
      * @param nztrue Number of event observables of the non-augmented model
      * @param ne Number of events
+     * @param ne_solver Number of events that require root-finding
      * @param nspl Number of splines
      * @param nJ Number of objective functions
      * @param nw Number of repeating elements
@@ -58,11 +59,12 @@ struct ModelDimensions {
         int const nx_rdata, int const nxtrue_rdata, int const nx_solver,
         int const nxtrue_solver, int const nx_solver_reinit, int const np,
         int const nk, int const ny, int const nytrue, int const nz,
-        int const nztrue, int const ne, int const nspl, int const nJ,
-        int const nw, int const ndwdx, int const ndwdp, int const ndwdw,
-        int const ndxdotdw, std::vector<int> ndJydy, int const ndxrdatadxsolver,
-        int const ndxrdatadtcl, int const ndtotal_cldx_rdata, int const nnz,
-        int const ubw, int const lbw
+        int const nztrue, int const ne, int const ne_solver, int const nspl,
+        int const nJ, int const nw, int const ndwdx, int const ndwdp,
+        int const ndwdw, int const ndxdotdw, std::vector<int> ndJydy,
+        int const ndxrdatadxsolver, int const ndxrdatadtcl,
+        int const ndtotal_cldx_rdata, int const nnz, int const ubw,
+        int const lbw
     )
         : nx_rdata(nx_rdata)
         , nxtrue_rdata(nxtrue_rdata)
@@ -76,6 +78,7 @@ struct ModelDimensions {
         , nz(nz)
         , nztrue(nztrue)
         , ne(ne)
+        , ne_solver(ne_solver)
         , nspl(nspl)
         , nw(nw)
         , ndwdx(ndwdx)
@@ -104,6 +107,8 @@ struct ModelDimensions {
         Expects(nztrue >= 0);
         Expects(nztrue <= nz);
         Expects(ne >= 0);
+        Expects(ne_solver >= 0);
+        Expects(ne >= ne_solver);
         Expects(nspl >= 0);
         Expects(nw >= 0);
         Expects(ndwdx >= 0);
@@ -164,7 +169,10 @@ struct ModelDimensions {
     /** Number of events */
     int ne{0};
 
-    /** numer of spline functions in the model */
+    /** Number of events that require root-finding */
+    int ne_solver{0};
+
+    /** Number of spline functions in the model */
     int nspl{0};
 
     /** Number of common expressions */
diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h
index 6567ee5c95..e03e1867d1 100644
--- a/include/amici/model_ode.h
+++ b/include/amici/model_ode.h
@@ -39,6 +39,8 @@ class Model_ODE : public Model {
      * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit
      * @param ndxdotdx_explicit number of nonzero elements dxdotdx_explicit
      * @param w_recursion_depth Recursion depth of fw
+     * @param state_independent_events Map of events with state-independent
+     * triggers functions, mapping trigger timepoints to event indices.
      */
     Model_ODE(
         ModelDimensions const& model_dimensions,
@@ -46,12 +48,13 @@ class Model_ODE : public Model {
         const SecondOrderMode o2mode, std::vector<realtype> const& idlist,
         std::vector<int> const& z2event, bool const pythonGenerated = false,
         int const ndxdotdp_explicit = 0, int const ndxdotdx_explicit = 0,
-        int const w_recursion_depth = 0
+        int const w_recursion_depth = 0,
+        std::map<realtype, std::vector<int>> state_independent_events = {}
     )
         : Model(
             model_dimensions, simulation_parameters, o2mode, idlist, z2event,
             pythonGenerated, ndxdotdp_explicit, ndxdotdx_explicit,
-            w_recursion_depth
+            w_recursion_depth, state_independent_events
         ) {}
 
     void
diff --git a/include/amici/serialization.h b/include/amici/serialization.h
index 501b56618a..c2a69d4b3a 100644
--- a/include/amici/serialization.h
+++ b/include/amici/serialization.h
@@ -15,6 +15,7 @@
 #include <boost/iostreams/device/back_inserter.hpp>
 #include <boost/iostreams/stream.hpp>
 #include <boost/serialization/array.hpp>
+#include <boost/serialization/map.hpp>
 #include <boost/serialization/vector.hpp>
 
 /** @file serialization.h Helper functions and forward declarations for
@@ -143,6 +144,7 @@ void serialize(Archive& ar, amici::Model& m, unsigned int const /*version*/) {
     ar& m.sigma_res_;
     ar& m.steadystate_computation_mode_;
     ar& m.steadystate_sensitivity_mode_;
+    ar& m.state_independent_events_;
 }
 
 /**
@@ -260,6 +262,7 @@ void serialize(
     ar& m.nz;
     ar& m.nztrue;
     ar& m.ne;
+    ar& m.ne_solver;
     ar& m.nspl;
     ar& m.nw;
     ar& m.ndwdx;
diff --git a/matlab/@amimodel/generateC.m b/matlab/@amimodel/generateC.m
index 59eb5e37b4..869f386baf 100644
--- a/matlab/@amimodel/generateC.m
+++ b/matlab/@amimodel/generateC.m
@@ -163,6 +163,7 @@ function generateC(this)
 fprintf(fid,['                  ' num2str(this.nz) ',\n']);
 fprintf(fid,['                  ' num2str(this.nztrue) ',\n']);
 fprintf(fid,['                  ' num2str(this.nevent) ',\n']);
+fprintf(fid,['                  ' num2str(this.nevent) ',\n']);
 fprintf(fid,['                  0,\n']);
 fprintf(fid,['                  ' num2str(this.ng) ',\n']);
 fprintf(fid,['                  ' num2str(this.nw) ',\n']);
diff --git a/models/model_calvetti/model_calvetti.h b/models/model_calvetti/model_calvetti.h
index 828be82728..c8144bdf5e 100644
--- a/models/model_calvetti/model_calvetti.h
+++ b/models/model_calvetti/model_calvetti.h
@@ -45,6 +45,7 @@ class Model_model_calvetti : public amici::Model_DAE {
                   0,
                   0,
                   4,
+                  4,
                   0,
                   1,
                   38,
@@ -207,6 +208,6 @@ class Model_model_calvetti : public amici::Model_DAE {
 
 } // namespace model_model_calvetti
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_calvetti_h */
diff --git a/models/model_dirac/model_dirac.h b/models/model_dirac/model_dirac.h
index 7a762479b5..cfd943e456 100644
--- a/models/model_dirac/model_dirac.h
+++ b/models/model_dirac/model_dirac.h
@@ -45,6 +45,7 @@ class Model_model_dirac : public amici::Model_ODE {
                   0,
                   0,
                   2,
+                  2,
                   0,
                   1,
                   0,
@@ -204,6 +205,6 @@ class Model_model_dirac : public amici::Model_ODE {
 
 } // namespace model_model_dirac
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_dirac_h */
diff --git a/models/model_events/model_events.h b/models/model_events/model_events.h
index df4bb68ae7..ad6c976419 100644
--- a/models/model_events/model_events.h
+++ b/models/model_events/model_events.h
@@ -59,6 +59,7 @@ class Model_model_events : public amici::Model_ODE {
                   2,
                   2,
                   6,
+                  6,
                   0,
                   1,
                   0,
@@ -232,6 +233,6 @@ class Model_model_events : public amici::Model_ODE {
 
 } // namespace model_model_events
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_events_h */
diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint.h b/models/model_jakstat_adjoint/model_jakstat_adjoint.h
index fdac2a9f94..6d7601947a 100644
--- a/models/model_jakstat_adjoint/model_jakstat_adjoint.h
+++ b/models/model_jakstat_adjoint/model_jakstat_adjoint.h
@@ -49,6 +49,7 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE {
                   0,
                   0,
                   0,
+                  0,
                   1,
                   2,
                   1,
@@ -210,6 +211,6 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE {
 
 } // namespace model_model_jakstat_adjoint
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_jakstat_adjoint_h */
diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h
index 22ca276067..bfac0b3267 100644
--- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h
+++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h
@@ -49,6 +49,7 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE {
                   0,
                   0,
                   0,
+                  0,
                   18,
                   10,
                   2,
@@ -210,6 +211,6 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE {
 
 } // namespace model_model_jakstat_adjoint_o2
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_jakstat_adjoint_o2_h */
diff --git a/models/model_nested_events/model_nested_events.h b/models/model_nested_events/model_nested_events.h
index 9ff8f519fe..0ed43eedbd 100644
--- a/models/model_nested_events/model_nested_events.h
+++ b/models/model_nested_events/model_nested_events.h
@@ -48,6 +48,7 @@ class Model_model_nested_events : public amici::Model_ODE {
                   0,
                   0,
                   4,
+                  4,
                   0,
                   1,
                   0,
@@ -210,6 +211,6 @@ class Model_model_nested_events : public amici::Model_ODE {
 
 } // namespace model_model_nested_events
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_nested_events_h */
diff --git a/models/model_neuron/model_neuron.h b/models/model_neuron/model_neuron.h
index e8f6f5c21f..e744f57f65 100644
--- a/models/model_neuron/model_neuron.h
+++ b/models/model_neuron/model_neuron.h
@@ -62,6 +62,7 @@ class Model_model_neuron : public amici::Model_ODE {
                   1,
                   1,
                   1,
+                  1,
                   0,
                   1,
                   0,
@@ -238,6 +239,6 @@ class Model_model_neuron : public amici::Model_ODE {
 
 } // namespace model_model_neuron
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_neuron_h */
diff --git a/models/model_neuron_o2/model_neuron_o2.h b/models/model_neuron_o2/model_neuron_o2.h
index 23df2b9b33..a108e6284b 100644
--- a/models/model_neuron_o2/model_neuron_o2.h
+++ b/models/model_neuron_o2/model_neuron_o2.h
@@ -64,6 +64,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE {
                   5,
                   1,
                   1,
+                  1,
                   0,
                   5,
                   2,
@@ -242,6 +243,6 @@ class Model_model_neuron_o2 : public amici::Model_ODE {
 
 } // namespace model_model_neuron_o2
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_neuron_o2_h */
diff --git a/models/model_robertson/model_robertson.h b/models/model_robertson/model_robertson.h
index 7f4377d785..816dd2db32 100644
--- a/models/model_robertson/model_robertson.h
+++ b/models/model_robertson/model_robertson.h
@@ -47,6 +47,7 @@ class Model_model_robertson : public amici::Model_DAE {
                   0,
                   0,
                   0,
+                  0,
                   1,
                   1,
                   2,
@@ -209,6 +210,6 @@ class Model_model_robertson : public amici::Model_DAE {
 
 } // namespace model_model_robertson
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_robertson_h */
diff --git a/models/model_steadystate/model_steadystate.h b/models/model_steadystate/model_steadystate.h
index b61649f9c8..776b754b08 100644
--- a/models/model_steadystate/model_steadystate.h
+++ b/models/model_steadystate/model_steadystate.h
@@ -46,6 +46,7 @@ class Model_model_steadystate : public amici::Model_ODE {
                   0,
                   0,
                   0,
+                  0,
                   1,
                   2,
                   2,
@@ -204,6 +205,6 @@ class Model_model_steadystate : public amici::Model_ODE {
 
 } // namespace model_model_steadystate
 
-} // namespace amici 
+} // namespace amici
 
 #endif /* _amici_model_steadystate_h */
diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py
index b1fa02c421..ac857b366c 100644
--- a/python/sdist/amici/de_export.py
+++ b/python/sdist/amici/de_export.py
@@ -1425,13 +1425,24 @@ def num_expr(self) -> int:
         return len(self.sym("w"))
 
     def num_events(self) -> int:
+        """
+        Total number of Events (those for which root-functions are added and those without).
+
+        :return:
+            number of events
+        """
+        return len(self.sym("h"))
+
+    def num_events_solver(self) -> int:
         """
         Number of Events.
 
         :return:
             number of event symbols (length of the root vector in AMICI)
         """
-        return len(self.sym("h"))
+        return sum(
+            not event.triggers_at_fixed_timepoint() for event in self.events()
+        )
 
     def sym(self, name: str) -> sp.Matrix:
         """
@@ -1750,6 +1761,16 @@ def parse_events(self) -> None:
             # add roots of heaviside functions
             self.add_component(root)
 
+        # re-order events - first those that require root tracking, then the others
+        self._events = list(
+            chain(
+                itertools.filterfalse(
+                    Event.triggers_at_fixed_timepoint, self._events
+                ),
+                filter(Event.triggers_at_fixed_timepoint, self._events),
+            )
+        )
+
     def get_appearance_counts(self, idxs: List[int]) -> List[int]:
         """
         Counts how often a state appears in the time derivative of
@@ -3642,6 +3663,7 @@ def _write_model_header_cpp(self) -> None:
             "NZ": self.model.num_eventobs(),
             "NZTRUE": self.model.num_eventobs(),
             "NEVENT": self.model.num_events(),
+            "NEVENT_SOLVER": self.model.num_events_solver(),
             "NOBJECTIVE": "1",
             "NSPL": len(self.model.splines),
             "NW": len(self.model.sym("w")),
@@ -3736,6 +3758,7 @@ def _write_model_header_cpp(self) -> None:
                 )
             ),
             "Z2EVENT": ", ".join(map(str, self.model._z2event)),
+            "STATE_INDEPENDENT_EVENTS": self._get_state_independent_event_intializer(),
             "ID": ", ".join(
                 (
                     str(float(isinstance(s, DifferentialState)))
@@ -3871,6 +3894,27 @@ def _get_symbol_id_initializer_list(self, name: str) -> str:
             for idx, symbol in enumerate(self.model.sym(name))
         )
 
+    def _get_state_independent_event_intializer(self) -> str:
+        """Get initializer list for state independent events in amici::Model."""
+        map_time_to_event_idx = {}
+        for event_idx, event in enumerate(self.model.events()):
+            if not event.triggers_at_fixed_timepoint():
+                continue
+            trigger_time = float(event.get_trigger_time())
+            try:
+                map_time_to_event_idx[trigger_time].append(event_idx)
+            except KeyError:
+                map_time_to_event_idx[trigger_time] = [event_idx]
+
+        def vector_initializer(v):
+            """std::vector initializer list with elements from `v`"""
+            return f"{{{', '.join(map(str, v))}}}"
+
+        return ", ".join(
+            f"{{{trigger_time}, {vector_initializer(event_idxs)}}}"
+            for trigger_time, event_idxs in map_time_to_event_idx.items()
+        )
+
     def _write_c_make_file(self):
         """Write CMake ``CMakeLists.txt`` file for this model."""
         sources = "\n".join(
diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py
index 77d9013ad2..72ce4a95ba 100644
--- a/python/sdist/amici/de_model.py
+++ b/python/sdist/amici/de_model.py
@@ -8,6 +8,7 @@
 from .import_utils import (
     RESERVED_SYMBOLS,
     ObservableTransformation,
+    amici_time_symbol,
     cast_to_sym,
     generate_measurement_symbol,
     generate_regularization_symbol,
@@ -696,6 +697,9 @@ def __init__(
         self._state_update = state_update
         self._initial_value = initial_value
 
+        # expression(s) for the timepoint(s) at which the event triggers
+        self._t_root = sp.solve(self.get_val(), amici_time_symbol)
+
     def get_initial_value(self) -> bool:
         """
         Return the initial value for the root function.
@@ -713,3 +717,20 @@ def __eq__(self, other):
         return self.get_val() == other.get_val() and (
             self.get_initial_value() == other.get_initial_value()
         )
+
+    def triggers_at_fixed_timepoint(self) -> bool:
+        """Check whether the event triggers at a (single) fixed time-point."""
+        if len(self._t_root) != 1:
+            return False
+        return self._t_root[0].is_Number
+
+    def get_trigger_time(self) -> sp.Float:
+        """Get the time at which the event triggers.
+
+        Only for events that trigger at a single fixed time-point.
+        """
+        if not self.triggers_at_fixed_timepoint():
+            raise NotImplementedError(
+                "This event does not trigger at a fixed timepoint."
+            )
+        return self._t_root[0]
diff --git a/python/tests/test_de_model.py b/python/tests/test_de_model.py
new file mode 100644
index 0000000000..7dec534da9
--- /dev/null
+++ b/python/tests/test_de_model.py
@@ -0,0 +1,35 @@
+import sympy as sp
+from amici.de_model import Event
+from amici.import_utils import amici_time_symbol
+
+
+def test_event_trigger_time():
+    e = Event(
+        sp.Symbol("event1"), "event name", amici_time_symbol - 10, sp.Float(0)
+    )
+    assert e.triggers_at_fixed_timepoint() is True
+    assert e.get_trigger_time() == 10
+
+    # fixed, but multiple timepoints - not (yet) supported
+    e = Event(
+        sp.Symbol("event1"),
+        "event name",
+        sp.sin(amici_time_symbol),
+        sp.Float(0),
+    )
+    assert e.triggers_at_fixed_timepoint() is False
+
+    e = Event(
+        sp.Symbol("event1"), "event name", amici_time_symbol / 2, sp.Float(0)
+    )
+    assert e.triggers_at_fixed_timepoint() is True
+    assert e.get_trigger_time() == 0
+
+    # parameter-dependent triggers - not (yet) supported
+    e = Event(
+        sp.Symbol("event1"),
+        "event name",
+        amici_time_symbol - sp.Symbol("delay"),
+        sp.Float(0),
+    )
+    assert e.triggers_at_fixed_timepoint() is False
diff --git a/python/tests/test_events.py b/python/tests/test_events.py
index d2a177bded..065cdeb126 100644
--- a/python/tests/test_events.py
+++ b/python/tests/test_events.py
@@ -1,8 +1,12 @@
 """Tests for SBML events, including piecewise expressions."""
 from copy import deepcopy
 
+import amici
 import numpy as np
 import pytest
+from amici.antimony_import import antimony2amici
+from amici.gradient_check import check_derivatives
+from amici.testing import TemporaryDirectoryWinSafe as TemporaryDirectory
 from amici.testing import skip_on_valgrind
 from util import (
     check_trajectories_with_forward_sensitivities,
@@ -704,3 +708,40 @@ def expm(x):
     from mpmath import expm
 
     return np.array(expm(x).tolist()).astype(float)
+
+
+def test_handling_of_fixed_time_point_event_triggers():
+    """Test handling of events without solver-tracked root functions."""
+    ant_model = """
+    model test_events_time_based
+        event_target = 0
+        bolus = 1
+        at (time > 1): event_target = 1
+        at (time > 2): event_target = event_target + bolus
+        at (time > 3): event_target = 3
+    end
+    """
+    module_name = "test_events_time_based"
+    with TemporaryDirectory(prefix=module_name, delete=False) as outdir:
+        antimony2amici(
+            ant_model,
+            model_name=module_name,
+            output_dir=outdir,
+            verbose=True,
+        )
+        model_module = amici.import_model_module(
+            module_name=module_name, module_path=outdir
+        )
+        amici_model = model_module.getModel()
+        assert amici_model.ne == 3
+        assert amici_model.ne_solver == 0
+        amici_model.setTimepoints(np.linspace(0, 4, 200))
+        amici_solver = amici_model.getSolver()
+        rdata = amici.runAmiciSimulation(amici_model, amici_solver)
+        assert rdata.status == amici.AMICI_SUCCESS
+        assert (rdata.x[rdata.ts < 1] == 0).all()
+        assert (rdata.x[(rdata.ts >= 1) & (rdata.ts < 2)] == 1).all()
+        assert (rdata.x[(rdata.ts >= 2) & (rdata.ts < 3)] == 2).all()
+        assert (rdata.x[(rdata.ts >= 3)] == 3).all()
+
+        check_derivatives(amici_model, amici_solver, edata=None)
diff --git a/python/tests/test_swig_interface.py b/python/tests/test_swig_interface.py
index a746552b55..6a7cfec855 100644
--- a/python/tests/test_swig_interface.py
+++ b/python/tests/test_swig_interface.py
@@ -315,6 +315,7 @@ def test_unhandled_settings(pysb_example_presimulation_module):
         "setParametersByIdRegex",
         "setParametersByNameRegex",
         "setInitialStateSensitivities",
+        "get_trigger_timepoints",
     ]
     from amici.swig_wrappers import model_instance_settings
 
diff --git a/src/exception.cpp b/src/exception.cpp
index 0e031ccca0..81ec938d29 100644
--- a/src/exception.cpp
+++ b/src/exception.cpp
@@ -33,14 +33,20 @@ void AmiException::storeMessage(char const* fmt, va_list argptr) {
     vsnprintf(msg_.data(), msg_.size(), fmt, argptr);
 }
 
-CvodeException::CvodeException(int const error_code, char const* function)
+CvodeException::CvodeException(
+    int const error_code, char const* function, char const* extra
+)
     : AmiException(
-        "Cvode routine %s failed with error code %i", function, error_code
+        "CVODE routine %s failed with error code %i. %s", function, error_code,
+        extra ? extra : ""
     ) {}
 
-IDAException::IDAException(int const error_code, char const* function)
+IDAException::IDAException(
+    int const error_code, char const* function, char const* extra
+)
     : AmiException(
-        "IDA routine %s failed with error code %i", function, error_code
+        "IDA routine %s failed with error code %i. %s", function, error_code,
+        extra ? extra : ""
     ) {}
 
 IntegrationFailure::IntegrationFailure(int code, realtype t)
diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp
index 13946547ef..c090f51280 100644
--- a/src/forwardproblem.cpp
+++ b/src/forwardproblem.cpp
@@ -12,6 +12,27 @@
 
 namespace amici {
 
+/**
+ * @brief Check if the next timepoint is too close to the current timepoint.
+ *
+ * Based on CVODES' `cvHin`.
+ * @param cur_t Current time.
+ * @param t_next Next stop time.
+ * @return True if too close, false otherwise.
+ */
+bool is_next_t_too_close(realtype cur_t, realtype t_next) {
+    auto tdiff = t_next - cur_t;
+    if(tdiff == 0.0)
+        return true;
+
+    auto tdist = std::fabs(tdiff);
+    auto tround = std::numeric_limits<realtype>::epsilon() * std::max(std::fabs(cur_t), std::fabs(t_next));
+    if (tdist < 2.0 * tround)
+        return true;
+
+    return false;
+}
+
 ForwardProblem::ForwardProblem(
     ExpData const* edata, Model* model, Solver* solver,
     SteadystateProblem const* preeq
@@ -109,30 +130,69 @@ void ForwardProblem::workForwardProblem() {
 
     /* store initial state and sensitivity*/
     initial_state_ = getSimulationState();
+    // store root information at t0
+    model->froot(t_, x_, dx_, rootvals_);
+
+    // get list of trigger timepoints for fixed-time triggered events
+    auto trigger_timepoints = model->get_trigger_timepoints();
+    auto it_trigger_timepoints = std::find_if(
+        trigger_timepoints.begin(), trigger_timepoints.end(),
+        [this](auto t) { return t > this->t_; }
+    );
 
     /* loop over timepoints */
     for (it_ = 0; it_ < model->nt(); it_++) {
-        auto nextTimepoint = model->getTimepoint(it_);
+        // next output time-point
+        auto next_t_out = model->getTimepoint(it_);
 
-        if (std::isinf(nextTimepoint))
+        if (std::isinf(next_t_out))
             break;
 
-        if (nextTimepoint > model->t0()) {
-            // Solve for nextTimepoint
-            while (t_ < nextTimepoint) {
-                int status = solver->run(nextTimepoint);
-                solver->writeSolution(&t_, x_, dx_, sx_, dx_);
+        if (next_t_out > model->t0()) {
+            // Solve for next output timepoint
+            while (t_ < next_t_out) {
+                if (is_next_t_too_close(t_, next_t_out)) {
+                    // next timepoint is too close to current timepoint.
+                    // we use the state of the current timepoint.
+                    break;
+                }
+
+                // next stop time is next output timepoint or next
+                // time-triggered event
+                auto next_t_event
+                    = it_trigger_timepoints != trigger_timepoints.end()
+                          ? *it_trigger_timepoints
+                          : std::numeric_limits<realtype>::infinity();
+                auto next_t_stop = std::min(next_t_out, next_t_event);
+
+                int status = solver->run(next_t_stop);
                 /* sx will be copied from solver on demand if sensitivities
                  are computed */
+                solver->writeSolution(&t_, x_, dx_, sx_, dx_);
+
                 if (status == AMICI_ILL_INPUT) {
-                    /* clustering of roots => turn off rootfinding */
+                    /* clustering of roots => turn off root-finding */
                     solver->turnOffRootFinding();
-                } else if (status == AMICI_ROOT_RETURN) {
+                } else if (status == AMICI_ROOT_RETURN || t_ == next_t_event) {
+                    // solver-tracked or time-triggered event
+                    solver->getRootInfo(roots_found_.data());
+
+                    // check if we are at a trigger timepoint.
+                    // if so, set the root-found flag
+                    if (t_ == next_t_event) {
+                        for (auto ie : model->state_independent_events_[t_]) {
+                            // determine direction of root crossing from
+                            // root function value at the previous event
+                            roots_found_[ie] = std::copysign(1, -rootvals_[ie]);
+                        }
+                        ++it_trigger_timepoints;
+                    }
+
                     handleEvent(&tlastroot_, false, false);
                 }
             }
         }
-        handleDataPoint(it_);
+        handleDataPoint(next_t_out);
     }
 
     /* fill events */
@@ -156,13 +216,9 @@ void ForwardProblem::handleEvent(
     /* store Heaviside information at event occurrence */
     model->froot(t_, x_, dx_, rootvals_);
 
-    /* store timepoint at which the event occurred*/
+    /* store timepoint at which the event occurred */
     discs_.push_back(t_);
 
-    /* extract and store which events occurred */
-    if (!seflag && !initial_event) {
-        solver->getRootInfo(roots_found_.data());
-    }
     root_idx_.push_back(roots_found_);
 
     rval_tmp_ = rootvals_;
@@ -324,11 +380,11 @@ void ForwardProblem::handle_secondary_event(realtype* tlastroot) {
     }
 }
 
-void ForwardProblem::handleDataPoint(int /*it*/) {
+void ForwardProblem::handleDataPoint(realtype t) {
     /* We only store the simulation state if it's not the initial state, as the
        initial state is stored anyway and we want to avoid storing it twice */
-    if (t_ != model->t0() && timepoint_states_.count(t_) == 0)
-        timepoint_states_[t_] = getSimulationState();
+    if (t != model->t0() && timepoint_states_.count(t) == 0)
+        timepoint_states_[t] = getSimulationState();
     /* store diagnosis information for debugging */
     solver->storeDiagnosis();
 }
diff --git a/src/model.cpp b/src/model.cpp
index c50fcc60bf..218870cafe 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -178,12 +178,14 @@ Model::Model(
     SimulationParameters simulation_parameters, SecondOrderMode o2mode,
     std::vector<realtype> idlist, std::vector<int> z2event,
     bool const pythonGenerated, int const ndxdotdp_explicit,
-    int const ndxdotdx_explicit, int const w_recursion_depth
+    int const ndxdotdx_explicit, int const w_recursion_depth,
+    std::map<realtype, std::vector<int>> state_independent_events
 )
     : ModelDimensions(model_dimensions)
     , pythonGenerated(pythonGenerated)
     , o2mode(o2mode)
     , idlist(std::move(idlist))
+    , state_independent_events_(std::move(state_independent_events))
     , derived_state_(model_dimensions)
     , z2event_(std::move(z2event))
     , state_is_non_negative_(nx_solver, false)
@@ -297,6 +299,7 @@ bool operator==(ModelDimensions const& a, ModelDimensions const& b) {
            && (a.nx_solver_reinit == b.nx_solver_reinit) && (a.np == b.np)
            && (a.nk == b.nk) && (a.ny == b.ny) && (a.nytrue == b.nytrue)
            && (a.nz == b.nz) && (a.nztrue == b.nztrue) && (a.ne == b.ne)
+           && (a.ne_solver == b.ne_solver) && (a.nspl == b.nspl)
            && (a.nw == b.nw) && (a.ndwdx == b.ndwdx) && (a.ndwdp == b.ndwdp)
            && (a.ndwdw == b.ndwdw) && (a.ndxdotdw == b.ndxdotdw)
            && (a.ndJydy == b.ndJydy) && (a.nnz == b.nnz) && (a.nJ == b.nJ)
@@ -3071,6 +3074,20 @@ void Model::fstotal_cl(
     );
 }
 
+std::vector<double> Model::get_trigger_timepoints() const {
+    std::vector<double> trigger_timepoints(
+        state_independent_events_.size(), 0.0
+    );
+    // collect keys from state_independent_events_ which are the trigger
+    // timepoints
+    auto it = trigger_timepoints.begin();
+    for (auto const& kv : state_independent_events_) {
+        *(it++) = kv.first;
+    }
+    std::sort(trigger_timepoints.begin(), trigger_timepoints.end());
+    return trigger_timepoints;
+}
+
 const_N_Vector Model::computeX_pos(const_N_Vector x) {
     if (any_state_non_negative_) {
         for (int ix = 0; ix < derived_state_.x_pos_tmp_.getLength(); ++ix) {
diff --git a/src/model_header.template.h b/src/model_header.template.h
index af05c8ccc5..932fdeb1a0 100644
--- a/src/model_header.template.h
+++ b/src/model_header.template.h
@@ -121,6 +121,7 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER {
                   TPL_NZ,                                  // nz
                   TPL_NZTRUE,                              // nztrue
                   TPL_NEVENT,                              // nevent
+                  TPL_NEVENT_SOLVER,                       // nevent_solver
                   TPL_NSPL,                                // nspl
                   TPL_NOBJECTIVE,                          // nobjective
                   TPL_NW,                                  // nw
@@ -146,7 +147,8 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER {
               true,                                        // pythonGenerated
               TPL_NDXDOTDP_EXPLICIT,                       // ndxdotdp_explicit
               TPL_NDXDOTDX_EXPLICIT,                       // ndxdotdx_explicit
-              TPL_W_RECURSION_DEPTH                        // w_recursion_depth
+              TPL_W_RECURSION_DEPTH,                       // w_recursion_depth
+              {TPL_STATE_INDEPENDENT_EVENTS}               // state-independent events
           ) {
                  root_initial_values_ = std::vector<bool>(
                      rootInitialValues.begin(), rootInitialValues.end()
diff --git a/src/solver.cpp b/src/solver.cpp
index c114623050..22e1723640 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -159,7 +159,7 @@ void Solver::setup(
     /* activates stability limit detection */
     setStabLimDet(stldet_);
 
-    rootInit(model->ne);
+    rootInit(model->ne_solver);
 
     if (nx() == 0)
         return;
diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp
index 7157302c9e..b53be68c2e 100644
--- a/src/solver_cvodes.cpp
+++ b/src/solver_cvodes.cpp
@@ -13,6 +13,8 @@
 #include <colamd.h>
 #include <klu.h>
 
+#include <sstream>
+
 #define ZERO RCONST(0.0)
 #define ONE RCONST(1.0)
 #define FOUR RCONST(4.0)
@@ -490,13 +492,16 @@ void CVodeSolver::reInitPostProcess(
     if (status == CV_ROOT_RETURN)
         throw CvodeException(
             status,
-            "CVode returned a root after "
-            "reinitialization. The initial step-size after the event or "
-            "heaviside function is too small. To fix this, increase absolute "
+            "CVode returned a root after reinitialization. "
+            "The initial step-size after the event or "
+            "Heaviside function is too small. To fix this, increase absolute "
             "and relative tolerances!"
         );
-    if (status != CV_SUCCESS)
-        throw CvodeException(status, "reInitPostProcess");
+    if (status != CV_SUCCESS) {
+        std::stringstream msg;
+        msg<<"tout: "<<tout<<", t: "<<*t<<".";
+        throw CvodeException(status, "reInitPostProcess", msg.str().c_str());
+    }
 
     cv_mem->cv_nst = nst_tmp + 1;
     if (cv_mem->cv_adjMallocDone == SUNTRUE) {
@@ -515,7 +520,7 @@ void CVodeSolver::reInitPostProcess(
         dt_mem[cv_mem->cv_nst % ca_mem->ca_nsteps]->t = *t;
         ca_mem->ca_IMstore(cv_mem, dt_mem[cv_mem->cv_nst % ca_mem->ca_nsteps]);
 
-        /* Set t1 field of the current ckeck point structure
+        /* Set t1 field of the current check point structure
          for the case in which there will be no future
          check points */
         ca_mem->ck_mem->ck_t1 = *t;
@@ -1066,9 +1071,17 @@ static int froot(realtype t, N_Vector x, realtype* root, void* user_data) {
     auto model = dynamic_cast<Model_ODE*>(typed_udata->first);
     Expects(model);
 
-    model->froot(t, x, gsl::make_span<realtype>(root, model->ne));
+    if (model->ne != model->ne_solver) {
+        // temporary buffer to store all root function values, not only the ones
+        // tracked by the solver
+        static std::vector<realtype> root_buffer(model->ne, 0.0);
+        model->froot(t, x, root_buffer);
+        std::copy_n(root_buffer.begin(), model->ne_solver, root);
+    } else {
+        model->froot(t, x, gsl::make_span<realtype>(root, model->ne_solver));
+    }
     return model->checkFinite(
-        gsl::make_span<realtype>(root, model->ne), ModelQuantity::root
+        gsl::make_span<realtype>(root, model->ne_solver), ModelQuantity::root
     );
 }
 
diff --git a/tests/cpp/unittests/testExpData.cpp b/tests/cpp/unittests/testExpData.cpp
index 416a41227b..d6e1a6fff2 100644
--- a/tests/cpp/unittests/testExpData.cpp
+++ b/tests/cpp/unittests/testExpData.cpp
@@ -4,8 +4,6 @@
 #include <amici/model_ode.h>
 #include <amici/symbolic_functions.h>
 
-#include <cmath>
-#include <cstring>
 #include <exception>
 #include <vector>
 
@@ -49,6 +47,7 @@ class ExpDataTest : public ::testing::Test {
             nz,        // nz
             nz,        // nztrue
             nmaxevent, // ne
+            0,         // ne_solver
             0,         // nspl
             0,         // nJ
             0,         // nw
diff --git a/tests/cpp/unittests/testMisc.cpp b/tests/cpp/unittests/testMisc.cpp
index 80d2c3bc36..14af3a7a82 100644
--- a/tests/cpp/unittests/testMisc.cpp
+++ b/tests/cpp/unittests/testMisc.cpp
@@ -65,6 +65,7 @@ class ModelTest : public ::testing::Test {
             nz,        // nz
             nz,        // nztrue
             nmaxevent, // ne
+            0,         // ne_solver
             0,         // nspl
             0,         // nJ
             0,         // nw
@@ -303,6 +304,7 @@ class SolverTest : public ::testing::Test {
             nz,        // nz
             nz,        // nztrue
             ne,        // ne
+            0,         // ne_solver
             0,         // nspl
             0,         // nJ
             0,         // nw
diff --git a/tests/cpp/unittests/testSerialization.cpp b/tests/cpp/unittests/testSerialization.cpp
index a516de0880..f59f04d9c7 100644
--- a/tests/cpp/unittests/testSerialization.cpp
+++ b/tests/cpp/unittests/testSerialization.cpp
@@ -142,6 +142,7 @@ TEST(ModelSerializationTest, ToFile)
             nz,        // nz
             nz,        // nztrue
             ne,        // ne
+            0,         // ne_solver
             0,         // nspl
             0,         // nJ
             9,         // nw
@@ -207,6 +208,7 @@ TEST(ReturnDataSerializationTest, ToString)
             nz,        // nz
             nz,        // nztrue
             ne,        // ne
+            0,         // ne_solver
             0,         // nspl
             0,         // nJ
             9,         // nw