diff --git a/.github/workflows/test_petab_test_suite.yml b/.github/workflows/test_petab_test_suite.yml
index f089b544b9..4594317582 100644
--- a/.github/workflows/test_petab_test_suite.yml
+++ b/.github/workflows/test_petab_test_suite.yml
@@ -71,7 +71,8 @@ jobs:
       - name: Run PEtab-related unit tests
         run: |
           source ./build/venv/bin/activate \
-            && pytest --cov-report=xml --cov=./ python/tests/test_*petab*.py
+            && pytest --cov-report=xml:coverage.xml \
+              --cov=./ python/tests/test_*petab*.py
 
       # run test models
       - name: Run PEtab test suite
@@ -79,12 +80,15 @@ jobs:
         run: |
           source ./build/venv/bin/activate \
           && AMICI_PARALLEL_COMPILE=2 pytest -v \
-            --cov-report=xml --cov-append --cov=amici tests/petab_test_suite/
+            --cov-report=xml:coverage.xml \
+            --cov-append \
+            --cov=amici \
+            tests/petab_test_suite/
 
       - name: Codecov
-        uses: codecov/codecov-action@v2
+        uses: codecov/codecov-action@v3.1.0
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
-          file: ./coverage.xml
+          file: coverage.xml
           flags: petab
           fail_ci_if_error: true
diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml
index d7fbb7d5ae..69b09e9dde 100644
--- a/.github/workflows/test_python_cplusplus.yml
+++ b/.github/workflows/test_python_cplusplus.yml
@@ -104,12 +104,13 @@ jobs:
         scripts/runNotebook.sh documentation/GettingStarted.ipynb
 
     - name: Codecov Python
-      uses: codecov/codecov-action@v2
+      uses: codecov/codecov-action@v3.1.0
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
-        file: ./build/coverage_py.xml
+        file: build/coverage_py.xml
         flags: python
         fail_ci_if_error: true
+        verbose: true
 
     - name: lcov
       run: |
@@ -117,15 +118,15 @@ jobs:
         -d ${AMICI_DIR}/build/CMakeFiles/amici.dir/src \
         -b ${AMICI_DIR} -c -o coverage_cpp.info \
         && lcov --compat-libtool --no-external \
-        -d ${AMICI_DIR}/python/sdist/build/temp.linux-x86_64-$(python3 --version | sed -E 's/.*([0-9]+\.[0-9]+)\..*/\1/')/amici/src \
+        -d ${AMICI_DIR}/python/sdist/build/temp.linux-x86_64-$(python -c "import sys; print(sys.implementation.cache_tag)")/amici/src \
         -b ${AMICI_DIR}/python/sdist -c -o coverage_py.info \
         && lcov -a coverage_cpp.info -a coverage_py.info -o coverage.info
 
     - name: Codecov CPP
-      uses: codecov/codecov-action@v2
+      uses: codecov/codecov-action@v3.1.0
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
-        file: ./coverage.info
+        file: coverage.info
         flags: cpp
         fail_ci_if_error: true
 
diff --git a/.github/workflows/test_sbml_semantic_test_suite.yml b/.github/workflows/test_sbml_semantic_test_suite.yml
index 458f010fd0..7901ef8640 100644
--- a/.github/workflows/test_sbml_semantic_test_suite.yml
+++ b/.github/workflows/test_sbml_semantic_test_suite.yml
@@ -53,9 +53,9 @@ jobs:
         path: tests/amici-semantic-results
 
     - name: Codecov SBMLSuite
-      uses: codecov/codecov-action@v2
+      uses: codecov/codecov-action@v3.1.0
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
-        file: ./coverage_SBMLSuite.xml
+        file: coverage_SBMLSuite.xml
         flags: sbmlsuite
         fail_ci_if_error: true
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4bde64a0b1..07e7c3c6f8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,30 @@
 
 ## v0.X Series
 
+### v0.11.29 (2022-05-06)
+
+## What's Changed
+
+Features:
+* Performance: Limit newton step convergence check by @FFroehlich in
+  https://github.com/AMICI-dev/AMICI/pull/1780
+* More informative NaN/Inf warnings by @dweindl in
+  https://github.com/AMICI-dev/AMICI/pull/1640
+* SBML import can now handle initial events by @FFroehlich in
+  https://github.com/AMICI-dev/AMICI/pull/1789
+
+Fixes:
+* Avoid error if no measurements in PEtab problem; fixed type handling in
+  PEtab parameter mapping by @dilpath in
+  https://github.com/AMICI-dev/AMICI/pull/1783
+* Fixed substitution of expressions in root and stau by @dilpath in
+  https://github.com/AMICI-dev/AMICI/pull/1784
+* Workaround for PEtab problems with state-dependent noise models by @dweindl
+  in https://github.com/AMICI-dev/AMICI/pull/1791
+
+**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.11.28...v0.11.29
+
+
 ### v0.11.28 (2022-04-08)
 
 New features:
diff --git a/documentation/python_interface.rst b/documentation/python_interface.rst
index b4137c27e3..02cc22af35 100644
--- a/documentation/python_interface.rst
+++ b/documentation/python_interface.rst
@@ -26,7 +26,7 @@ AMICI can import :term:`SBML` models via the
 Status of SBML support in Python-AMICI
 ++++++++++++++++++++++++++++++++++++++
 
-Python-AMICI currently **passes 1014 out of the 1821 (~56%) test cases** from
+Python-AMICI currently **passes 1030 out of the 1821 (~57%) test cases** from
 the semantic
 `SBML Test Suite <https://github.com/sbmlteam/sbml-test-suite/>`_
 (`current status <https://github.com/AMICI-dev/AMICI/actions>`_).
diff --git a/include/amici/forwardproblem.h b/include/amici/forwardproblem.h
index 712377d959..38c3ba1785 100644
--- a/include/amici/forwardproblem.h
+++ b/include/amici/forwardproblem.h
@@ -290,9 +290,11 @@ class ForwardProblem {
      *
      * @param tlastroot pointer to the timepoint of the last event
      * @param seflag Secondary event flag
+     * @param initial_event initial event flag
      */
 
-    void handleEvent(realtype *tlastroot,bool seflag);
+    void handleEvent(realtype *tlastroot, bool seflag,
+                     bool initial_event);
 
     /**
      * @brief Extract output information for events
diff --git a/include/amici/misc.h b/include/amici/misc.h
index ee069b5224..329ca9a48d 100644
--- a/include/amici/misc.h
+++ b/include/amici/misc.h
@@ -187,6 +187,17 @@ class ContextManager{
     ContextManager(ContextManager &&other) = delete;
 };
 
+
+/**
+ * @brief Convert a flat index to a pair of row/column indices,
+ * assuming row-major order.
+ * @param flat_idx flat index
+ * @param num_cols number of columns of referred to matrix
+ * @return row index, column index
+ */
+auto unravel_index(size_t flat_idx, size_t num_cols)
+    -> std::pair<size_t, size_t>;
+
 } // namespace amici
 
 #endif // AMICI_MISC_H
diff --git a/include/amici/model.h b/include/amici/model.h
index fac088af2c..1897f9b295 100644
--- a/include/amici/model.h
+++ b/include/amici/model.h
@@ -34,6 +34,42 @@ void serialize(Archive &ar, amici::Model &m, unsigned int version);
 
 namespace amici {
 
+/**
+ * @brief Describes the various model quantities
+ */
+enum class ModelQuantity {
+    J,
+    JB,
+    Jv,
+    JvB,
+    JDiag,
+    sx,
+    sy,
+    ssigmay,
+    xdot,
+    sxdot,
+    xBdot,
+    x0_rdata,
+    x0,
+    x_rdata,
+    x,
+    dwdw,
+    dwdx,
+    dwdp,
+    y,
+    dydp,
+    dydx,
+    w,
+    root,
+    qBdot,
+    qBdot_ss,
+    xBdot_ss,
+    JSparseB_ss,
+};
+
+extern const std::map<ModelQuantity, std::string> model_quantity_to_str;
+
+
 /**
  * @brief The Model class represents an AMICI ODE/DAE model.
  *
@@ -61,7 +97,8 @@ class Model : public AbstractModel, public ModelDimensions {
           SimulationParameters simulation_parameters,
           amici::SecondOrderMode o2mode,
           std::vector<amici::realtype> idlist,
-          std::vector<int> z2event, bool pythonGenerated = false,
+          std::vector<int> z2event,
+          bool pythonGenerated = false,
           int ndxdotdp_explicit = 0, int ndxdotdx_explicit = 0,
           int w_recursion_depth = 0);
 
@@ -166,9 +203,11 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param sdx Reference to time derivative of state sensitivities (DAE only)
      * @param computeSensitivities Flag indicating whether sensitivities are to
      * be computed
+     * @param roots_found boolean indicators indicating whether roots were found at t0 by this fun
      */
     void initialize(AmiVector &x, AmiVector &dx, AmiVectorArray &sx,
-                    AmiVectorArray &sdx, bool computeSensitivities);
+                    AmiVectorArray &sdx, bool computeSensitivities,
+                    std::vector<int> &roots_found);
 
     /**
      * @brief Initialize model properties.
@@ -200,8 +239,10 @@ class Model : public AbstractModel, public ModelDimensions {
      *
      * @param x Reference to state variables
      * @param dx Reference to time derivative of states (DAE only)
+     * @param roots_found boolean indicators indicating whether roots were found at t0 by this fun
      */
-    void initHeaviside(const AmiVector &x, const AmiVector &dx);
+    void initEvents(const AmiVector &x, const AmiVector &dx,
+                    std::vector<int> &roots_found);
 
     /**
      * @brief Get number of parameters wrt to which sensitivities are computed.
@@ -1201,18 +1242,43 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void updateHeavisideB(const int *rootsfound);
 
+
+    /**
+     * @brief Check if the given array has only finite elements.
+     *
+     * For (1D) spans.
+     *
+     * @param array
+     * @param model_quantity The model quantity `array` corresponds to
+     * @return
+     */
+    int checkFinite(gsl::span<const realtype> array,
+                           ModelQuantity model_quantity) const;
+    /**
+     * @brief Check if the given array has only finite elements.
+     *
+     * For flattened 2D arrays.
+     *
+     * @param array Flattened matrix
+     * @param model_quantity The model quantity `array` corresponds to
+     * @param num_cols Number of columns of the non-flattened matrix
+     * @return
+     */
+    int checkFinite(gsl::span<const realtype> array,
+                    ModelQuantity model_quantity,
+                    size_t num_cols) const;
+
     /**
      * @brief Check if the given array has only finite elements.
      *
-     * If not, try to give hints by which other fields this could be caused.
+     * For SUNMatrix.
      *
-     * @param array Array to check
-     * @param fun Name of the function that generated the values (for more
-     * informative messages).
-     * @return `amici::AMICI_RECOVERABLE_ERROR` if a NaN/Inf value was found,
-     * `amici::AMICI_SUCCESS` otherwise
+     * @param m Matrix to check
+     * @param model_quantity The model quantity `m` corresponds to
+     * @param t current timepoint
+     * @return
      */
-    int checkFinite(gsl::span<const realtype> array, const char *fun) const;
+    int checkFinite(SUNMatrix m, ModelQuantity model_quantity, realtype t) const;
 
     /**
      * @brief Set whether the result of every call to `Model::f*` should be
@@ -1803,6 +1869,11 @@ class Model : public AbstractModel, public ModelDimensions {
     /** vector of bools indicating whether state variables are to be assumed to
      * be positive */
     std::vector<bool> state_is_non_negative_;
+    
+    /** Vector of booleans indicating the initial boolean value for every event trigger function. Events at t0
+     * can only trigger if the initial value is set to `false`. Must be specified during model compilation by
+     * setting the `initialValue` attribute of an event trigger. */
+    std::vector<bool> root_initial_values_;
 
     /** boolean indicating whether any entry in stateIsNonNegative is `true` */
     bool any_state_non_negative_ {false};
diff --git a/include/amici/sundials_matrix_wrapper.h b/include/amici/sundials_matrix_wrapper.h
index b94667e39e..52d6c3377e 100644
--- a/include/amici/sundials_matrix_wrapper.h
+++ b/include/amici/sundials_matrix_wrapper.h
@@ -532,6 +532,16 @@ class SUNMatrixWrapper {
     bool ownmat = true;
 };
 
+
+/**
+ * @brief Convert a flat index to a pair of row/column indices.
+ * @param i flat index
+ * @param m referred to matrix
+ * @return row index, column index
+ */
+auto unravel_index(sunindextype i, SUNMatrix m)
+    -> std::pair<sunindextype, sunindextype>;
+
 } // namespace amici
 
 namespace gsl {
diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py
index a1f40409b6..b40e025524 100644
--- a/python/amici/ode_export.py
+++ b/python/amici/ode_export.py
@@ -822,7 +822,7 @@ def transform_dxdt_to_concentration(species_id, dxdt):
             else:
                 args += ['value']
             if symbol_name == SymbolId.EVENT:
-                args += ['state_update', 'event_observable']
+                args += ['state_update', 'event_observable', 'initial_value']
             if symbol_name == SymbolId.OBSERVABLE:
                 args += ['transformation']
 
@@ -1297,7 +1297,7 @@ def parse_events(self) -> None:
             expr.set_val(self._process_heavisides(expr.get_val(), roots))
 
         # remove all possible Heavisides from roots, which may arise from
-        # the substitution of `'w'` in `_get_unique_root`
+        # the substitution of `'w'` in `_collect_heaviside_roots`
         for root in roots:
             root.set_val(self._process_heavisides(root.get_val(), roots))
 
@@ -1589,34 +1589,43 @@ def _compute_equation(self, name: str) -> None:
         elif name == 'stau':
             self._eqs[name] = [
                 -self.eq('sroot')[ie, :] / self.eq('drootdt_total')[ie]
+                if not self.eq('drootdt_total')[ie].is_zero else
+                sp.zeros(*self.eq('sroot')[ie, :].shape)
                 for ie in range(self.num_events())
             ]
 
         elif name == 'deltasx':
             event_eqs = []
             for ie, event in enumerate(self._events):
-                if event._state_update is not None:
-                    # ====== chain rule for the state variables ===============
-                    # get xdot with expressions back-substituted
-                    tmp_eq = smart_multiply(
+
+                tmp_eq = sp.zeros(self.num_states_solver(), self.num_par())
+
+                # only add stau part if trigger is time-dependent
+                if not self.eq('drootdt_total')[ie].is_zero:
+                    tmp_eq += smart_multiply(
                         (self.sym('xdot_old') - self.sym('xdot')),
                         self.eq('stau')[ie])
-                    # construct an enhanced state sensitivity, which accounts
-                    # for the time point sensitivity as well
+
+                # only add deltax part if there is state update
+                if event._state_update is not None:
+                    # partial derivative for the parameters
+                    tmp_eq += self.eq('ddeltaxdp')[ie]
+
+                    # initial part of chain rule state variables
                     tmp_dxdp = self.sym('sx') * sp.ones(1, self.num_par())
-                    tmp_dxdp += smart_multiply(self.sym('xdot'),
-                                               self.eq('stau')[ie])
+
+                    # only add stau part if trigger is time-dependent
+                    if not self.eq('drootdt_total')[ie].is_zero:
+                        # chain rule for the time point
+                        tmp_eq += smart_multiply(self.eq('ddeltaxdt')[ie],
+                                                 self.eq('stau')[ie])
+
+                        # additional part of chain rule state variables
+                        tmp_dxdp += smart_multiply(self.sym('xdot'),
+                                                   self.eq('stau')[ie])
+                    # finish chain rule for the state variables
                     tmp_eq += smart_multiply(self.eq('ddeltaxdx')[ie],
                                              tmp_dxdp)
-                    # ====== chain rule for the time point ====================
-                    tmp_eq += smart_multiply(self.eq('ddeltaxdt')[ie],
-                                             self.eq('stau')[ie])
-                    # ====== partial derivative for the parameters ============
-                    tmp_eq += self.eq('ddeltaxdp')[ie]
-                else:
-                    tmp_eq = smart_multiply(
-                        (self.eq('xdot_old') - self.eq('xdot')),
-                        self.eq('stau')[ie])
 
                 event_eqs.append(tmp_eq)
 
@@ -2066,15 +2075,6 @@ def _get_unique_root(
             unique identifier for root, or ``None`` if the root is not
             time-dependent
         """
-        # substitute 'w' expressions into root expressions now, to avoid
-        # rewriting '{model_name}_root.cpp' and '{model_name}_stau.cpp' headers
-        # to include 'w.h'
-        w_sorted = toposort_symbols(dict(zip(
-            [expr.get_id() for expr in self._expressions],
-            [expr.get_val() for expr in self._expressions],
-        )))
-        root_found = root_found.subs(w_sorted)
-
         if not self._expr_is_time_dependent(root_found):
             return None
 
@@ -2115,6 +2115,18 @@ def _collect_heaviside_roots(
             elif arg.has(sp.Heaviside):
                 root_funs.extend(self._collect_heaviside_roots(arg.args))
 
+        # substitute 'w' expressions into root expressions now, to avoid
+        # rewriting '{model_name}_root.cpp' and '{model_name}_stau.cpp' headers
+        # to include 'w.h'
+        w_sorted = toposort_symbols(dict(zip(
+            [expr.get_id()  for expr in self._expressions],
+            [expr.get_val() for expr in self._expressions],
+        )))
+        root_funs = [
+            r.subs(w_sorted)
+            for r in root_funs
+        ]
+
         return root_funs
 
     def _process_heavisides(
@@ -2934,6 +2946,11 @@ def _write_model_header_cpp(self) -> None:
             'W_RECURSION_DEPTH': self.model._w_recursion_depth,
             'QUADRATIC_LLH': 'true'
                 if self.model._has_quadratic_nllh else 'false',
+            'ROOT_INITIAL_VALUES':
+                ', '.join([
+                    'true' if event.get_initial_value() else 'false'
+                    for event in self.model._events
+                ])
         }
 
         for func_name, func_info in self.functions.items():
diff --git a/python/amici/ode_model.py b/python/amici/ode_model.py
index a5885216f3..ed494cf3ea 100644
--- a/python/amici/ode_model.py
+++ b/python/amici/ode_model.py
@@ -2,17 +2,7 @@
 
 
 import sympy as sp
-import numpy as np
-import re
-import shutil
-import subprocess
-import sys
-import os
-import copy
 import numbers
-import logging
-import itertools
-import contextlib
 
 try:
     import pysb
@@ -20,25 +10,11 @@
     pysb = None
 
 from typing import (
-    Callable, Optional, Union, List, Dict, Tuple, SupportsFloat, Sequence,
-    Set, Any
+    Optional, Union, Dict, SupportsFloat, Set
 )
-from dataclasses import dataclass
-from string import Template
-from sympy.matrices.immutable import ImmutableDenseMatrix
-from sympy.matrices.dense import MutableDenseMatrix
-from sympy.logic.boolalg import BooleanAtom
-from itertools import chain
-from .cxxcodeprinter import AmiciCxxCodePrinter, get_switch_statement
-
-from . import (
-    amiciSwigPath, amiciSrcPath, amiciModulePath, __version__, __commit__,
-    sbml_import
-)
-from .logging import get_logger, log_execution_time, set_log_level
-from .constants import SymbolId
-from .import_utils import smart_subs_dict, toposort_symbols, \
-    ObservableTransformation, generate_measurement_symbol, RESERVED_SYMBOLS
+
+from .import_utils import ObservableTransformation, \
+    generate_measurement_symbol, RESERVED_SYMBOLS
 from .import_utils import cast_to_sym
 
 __all__ = [
@@ -46,6 +22,7 @@
     'ModelQuantity', 'Observable', 'Parameter', 'SigmaY', 'State'
 ]
 
+
 class ModelQuantity:
     """
     Base class for model components
@@ -166,14 +143,6 @@ def __init__(self,
         self._ncoeff: sp.Expr = coefficients[state_id]
         super(ConservationLaw, self).__init__(identifier, name, value)
 
-    def get_state(self) -> sp.Symbol:
-        """
-        Get the identifier of the state that this conservation law replaces
-
-        :return: identifier of the state
-        """
-        return self._state_id
-
     def get_ncoeff(self, state_id) -> Union[sp.Expr, int, float]:
         """
         Computes the normalized coefficient a_i/a_j where i is the index of
@@ -211,10 +180,6 @@ class State(ModelQuantity):
         algebraic formula that defines the temporal derivative of this state
 
     """
-
-    _dt: Union[sp.Expr, None] = None
-    _conservation_law: Union[sp.Expr, None] = None
-
     def __init__(self,
                  identifier: sp.Symbol,
                  name: str,
@@ -276,7 +241,7 @@ def get_dt(self) -> sp.Expr:
         """
         return self._dt
 
-    def get_free_symbols(self) -> Set[sp.Symbol]:
+    def get_free_symbols(self) -> Set[sp.Basic]:
         """
         Gets the set of free symbols in time derivative and initial conditions
 
@@ -307,8 +272,9 @@ def get_x_rdata(self):
 
     def get_dx_rdata_dx_solver(self, state_id):
         """
-        Returns the expression that allows computation of ``dx_rdata_dx_solver`` for this
-        state, accounting for conservation laws.
+        Returns the expression that allows computation of
+        ``dx_rdata_dx_solver`` for this state, accounting for conservation
+        laws.
 
         :return: dx_rdata_dx_solver expression
         """
@@ -514,7 +480,8 @@ def __init__(self,
                  name: str,
                  value: sp.Expr,
                  state_update: Union[sp.Expr, None],
-                 event_observable: Union[sp.Expr, None]):
+                 event_observable: Union[sp.Expr, None],
+                 initial_value: Optional[bool] = True):
         """
         Create a new Event instance.
 
@@ -534,15 +501,30 @@ def __init__(self,
         :param event_observable:
             formula a potential observable linked to the event
             (None for Heaviside functions, empty events without observable)
+
+        :param initial_value:
+            initial boolean value of the trigger function at t0. If set to
+            `False`, events may trigger at ``t==t0``, otherwise not.
         """
         super(Event, self).__init__(identifier, name, value)
         # add the Event specific components
         self._state_update = state_update
         self._observable = event_observable
+        self._initial_value = initial_value
+
+    def get_initial_value(self) -> bool:
+        """
+        Return the initial value for the root function.
+
+        :return:
+            initial value formula
+        """
+        return self._initial_value
 
     def __eq__(self, other):
         """
         Check equality of events at the level of trigger/root functions, as we
         need to collect unique root functions for ``roots.cpp``
         """
-        return self.get_val() == other.get_val()
+        return self.get_val() == other.get_val() and \
+            (self.get_initial_value() == other.get_initial_value())
diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py
index f1183b451f..bc520f9d20 100644
--- a/python/amici/petab_import.py
+++ b/python/amici/petab_import.py
@@ -40,6 +40,38 @@
 PREEQ_INDICATOR_ID = 'preequilibration_indicator'
 
 
+def _add_global_parameter(sbml_model: libsbml.Model,
+                          parameter_id: str,
+                          parameter_name: str = None,
+                          constant: bool = False,
+                          units: str = 'dimensionless',
+                          value: float = 0.0) -> libsbml.Parameter:
+    """Add new global parameter to SBML model
+
+    Arguments:
+        sbml_model: SBML model
+        parameter_id: ID of the new parameter
+        parameter_name: Name of the new parameter
+        constant: Is parameter constant?
+        units: SBML unit ID
+        value: parameter value
+
+    Returns:
+        The created parameter
+    """
+
+    if parameter_name is None:
+        parameter_name = parameter_id
+
+    p = sbml_model.createParameter()
+    p.setId(parameter_id)
+    p.setName(parameter_name)
+    p.setConstant(constant)
+    p.setValue(value)
+    p.setUnits(units)
+    return p
+
+
 def get_fixed_parameters(
         sbml_model: 'libsbml.Model',
         condition_df: Optional[pd.DataFrame] = None,
@@ -501,11 +533,12 @@ def import_model_sbml(
                            key=lambda symbol: symbol.name)
         for free_sym in free_syms:
             sym = str(free_sym)
-            if sbml_model.getElementBySId(sym) is None and sym != 'time':
+            if sbml_model.getElementBySId(sym) is None and sym != 'time' \
+                    and sym not in observables:
                 output_parameters[sym] = None
     logger.debug(f"Adding output parameters to model: {output_parameters}")
     for par in output_parameters.keys():
-        petab.add_global_parameter(sbml_model, par)
+        _add_global_parameter(sbml_model, par)
     # <EndWorkAround>
 
     # TODO: to parameterize initial states or compartment sizes, we currently
@@ -624,10 +657,15 @@ def get_observation_model(
         observables[oid] = {'name': name, 'formula': formula_obs}
         sigmas[oid] = formula_noise
 
-    # Replace observableIds occurring in error model definition
+    # PEtab does currently not allow observables in noiseFormula and AMICI
+    #  cannot handle states in sigma expressions. Therefore, where possible,
+    #  replace species occurring in error model definition by observableIds.
+    replacements = {
+        sp.sympify(observable['formula']): sp.Symbol(observable_id)
+        for observable_id, observable in observables.items()
+    }
     for observable_id, formula in sigmas.items():
-        repl = sp.sympify(formula).subs(
-            observable_id, observables[observable_id]['formula'])
+        repl = sp.sympify(formula).subs(replacements)
         sigmas[observable_id] = str(repl)
 
     noise_distrs = petab_noise_distributions_to_amici(observable_df)
diff --git a/python/amici/petab_objective.py b/python/amici/petab_objective.py
index 1e6dfb0e7f..0cb8b4d9d2 100644
--- a/python/amici/petab_objective.py
+++ b/python/amici/petab_objective.py
@@ -154,8 +154,11 @@ def simulate_petab(
     # Log results
     sim_cond = petab_problem.get_simulation_conditions_from_measurement_df()
     for i, rdata in enumerate(rdatas):
-        logger.debug(f"Condition: {sim_cond.iloc[i, :].values}, status: "
-                     f"{rdata['status']}, llh: {rdata['llh']}")
+        sim_cond_id = "N/A" if sim_cond.empty else sim_cond.iloc[i, :].values
+        logger.debug(
+            f"Condition: {sim_cond_id}, status: {rdata['status']}, "
+            f"llh: {rdata['llh']}"
+        )
 
     return {
         LLH: llh,
@@ -230,7 +233,7 @@ def create_parameterized_edatas(
 
 def create_parameter_mapping(
         petab_problem: petab.Problem,
-        simulation_conditions: Union[pd.DataFrame, Dict],
+        simulation_conditions: Union[pd.DataFrame, List[Dict]],
         scaled_parameters: bool,
         amici_model: AmiciModel,
 ) -> ParameterMapping:
@@ -254,6 +257,8 @@ def create_parameter_mapping(
     if simulation_conditions is None:
         simulation_conditions = \
             petab_problem.get_simulation_conditions_from_measurement_df()
+    if isinstance(simulation_conditions, list):
+        simulation_conditions = pd.DataFrame(data=simulation_conditions)
 
     # Because AMICI globalizes all local parameters during model import,
     # we need to do that here as well to prevent parameter mapping errors
diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index 616f82e102..8b38e06f87 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -484,13 +484,6 @@ def check_event_support(self) -> None:
             if trigger_sbml.getMath() is None:
                 logger.warning(f'Event {event_id} trigger has no trigger '
                                'expression, so a dummy trigger will be set.')
-            if not trigger_sbml.getInitialValue():
-                # True: event not executed if triggered at time == 0
-                # (corresponding to AMICI default). Raise if set to False.
-                raise SBMLException(
-                    f'Event {event_id} has a trigger that has an initial '
-                    'value of False. This is currently not supported in AMICI.'
-                )
 
             if not trigger_sbml.getPersistent():
                 raise SBMLException(
@@ -1001,10 +994,16 @@ def _convert_event_assignment_parameter_targets_to_species(self):
                     parameter_def = \
                         self.symbols[symbol_id].pop(parameter_target)
             if parameter_def is None:
-                raise AssertionError(
-                    'Unexpected error. The parameter target of an event '
-                    'assignment could not be found.'
+                # this happens for parameters that have initial assignments
+                # or are assignment rule targets
+                par = self.sbml.getElementBySId(str(parameter_target))
+                ia_init = self._get_element_initial_assignment(
+                    par.getId()
                 )
+                parameter_def = {
+                    'name': par.getName() if par.isSetName() else par.getId(),
+                    'value': par.getValue() if ia_init is None else ia_init
+                }
             # Fixed parameters are added as species such that they can be
             # targets of events.
             self.symbols[SymbolId.SPECIES][parameter_target] = {
@@ -1140,6 +1139,9 @@ def get_empty_bolus_value() -> sp.Float:
                 'value': trigger,
                 'state_update': sp.MutableDenseMatrix(bolus),
                 'event_observable': None,
+                'initial_value':
+                    trigger_sbml.getInitialValue() if trigger_sbml is not None
+                    else True,
             }
 
     @log_execution_time('processing SBML observables', logger)
diff --git a/src/amici.cpp b/src/amici.cpp
index 5be90dad00..cc20bd5927 100644
--- a/src/amici.cpp
+++ b/src/amici.cpp
@@ -179,7 +179,7 @@ AmiciApplication::runAmiciSimulation(Solver& solver,
                 throw;
             warningF("AMICI:simulation",
                      "AMICI forward simulation failed at t = %f: "
-                     "Maximum time exceeed.\n",
+                     "Maximum time exceeded.\n",
                      ex.time);
         } else {
             rdata->status = ex.error_code;
@@ -199,7 +199,7 @@ AmiciApplication::runAmiciSimulation(Solver& solver,
             warningF(
                 "AMICI:simulation",
                 "AMICI backward simulation failed when trying to solve until "
-                "t = %f: Maximum time exceeed.\n",
+                "t = %f: Maximum time exceeded.\n",
                 ex.time);
 
         } else {
@@ -305,6 +305,8 @@ AmiciApplication::errorF(const char* identifier, const char* format, ...) const
 int
 AmiciApplication::checkFinite(gsl::span<const realtype> array, const char* fun)
 {
+    Expects(array.size()
+            <= static_cast<unsigned>(std::numeric_limits<int>::max()));
 
     for (int idx = 0; idx < (int)array.size(); idx++) {
         if (isNaN(array[idx])) {
diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp
index 8147edce05..773ffc5c4b 100644
--- a/src/forwardproblem.cpp
+++ b/src/forwardproblem.cpp
@@ -52,15 +52,21 @@ void ForwardProblem::workForwardProblem() {
     if (!preequilibrated_)
         model->initialize(x_, dx_, sx_, sdx_,
                           solver->getSensitivityOrder() >=
-                          SensitivityOrder::first);
-    else if (model->ne)
-        model->initHeaviside(x_, dx_);
+                          SensitivityOrder::first,
+                          roots_found_);
+    else if (model->ne) {
+        model->initEvents(x_, dx_, roots_found_);
+    }
 
     /* compute initial time and setup solver for (pre-)simulation */
     auto t0 = model->t0();
     if (presimulate)
         t0 -= edata->t_presim;
     solver->setup(t0, model, x_, dx_, sx_, sdx_);
+    
+    if (model->ne && std::any_of(roots_found_.begin(), roots_found_.end(),
+                                 [](int rf){return rf==1;}))
+        handleEvent(&t0, false, true);
 
     /* perform presimulation if necessary */
     if (presimulate) {
@@ -69,9 +75,14 @@ void ForwardProblem::workForwardProblem() {
                                " is currently not implemented.");
         handlePresimulation();
         t_ = model->t0();
-        if (model->ne)
-            model->initHeaviside(x_, dx_);
+        if (model->ne) {
+            model->initEvents(x_, dx_, roots_found_);
+            if (std::any_of(roots_found_.begin(), roots_found_.end(),
+                            [](int rf){return rf==1;}))
+                handleEvent(&t0, false, true);
+        }
     }
+    
     /* when computing adjoint sensitivity analysis with presimulation,
      we need to store sx after the reinitialization after preequilibration
      but before reinitialization after presimulation. As presimulation with ASA
@@ -92,7 +103,6 @@ void ForwardProblem::workForwardProblem() {
     if (solver->computingFSA() || (solver->computingASA() && !presimulate ))
         sx_ = solver->getStateSensitivity(model->t0());
 
-
     /* store initial state and sensitivity*/
     initial_state_ = getSimulationState();
 
@@ -114,7 +124,7 @@ void ForwardProblem::workForwardProblem() {
                     /* clustering of roots => turn off rootfinding */
                     solver->turnOffRootFinding();
                 } else if (status == AMICI_ROOT_RETURN) {
-                    handleEvent(&tlastroot_, false);
+                    handleEvent(&tlastroot_, false, false);
                 }
             }
         }
@@ -138,7 +148,8 @@ void ForwardProblem::handlePresimulation()
 }
 
 
-void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
+void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag,
+                                 const bool initial_event) {
     /* store Heaviside information at event occurrence */
     model->froot(t_, x_, dx_, rootvals_);
 
@@ -146,14 +157,14 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
     discs_.push_back(t_);
 
     /* extract and store which events occurred */
-    if (!seflag) {
+    if (!seflag && !initial_event) {
         solver->getRootInfo(roots_found_.data());
     }
     root_idx_.push_back(roots_found_);
 
     rval_tmp_ = rootvals_;
 
-    if (!seflag) {
+    if (!seflag && !initial_event) {
         /* only check this in the first event fired, otherwise this will always
          * be true */
         if (t_ == *tlastroot) {
@@ -165,24 +176,24 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
         *tlastroot = t_;
     }
 
-    if(model->nz>0)
+    if(model->nz > 0)
         storeEvent();
 
     /* if we need to do forward sensitivities later on we need to store the old
      * x and the old xdot */
     if (solver->getSensitivityOrder() >= SensitivityOrder::first) {
         /* store x and xdot to compute jump in sensitivities */
-        x_old_ = x_;
+        x_old_.copy(x_);
     }
     if (solver->computingFSA()) {
         model->fxdot(t_, x_, dx_, xdot_);
-        xdot_old_ = xdot_;
-        dx_old_ = dx_;
+        xdot_old_.copy(xdot_);
+        dx_old_.copy(dx_);
         /* compute event-time derivative only for primary events, we get
          * into trouble with multiple simultaneously firing events here (but
          * is this really well defined then?), in that case just use the
          * last ie and hope for the best. */
-        if (!seflag) {
+        if (!seflag && !initial_event) {
             for (int ie = 0; ie < model->ne; ie++) {
                 if (roots_found_.at(ie) == 1) {
                     /* only consider transitions false -> true */
@@ -190,6 +201,8 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
                 }
             }
         }
+        if (initial_event) // t0 has no parameter dependency
+            std::fill(stau_.begin(), stau_.end(), 0.0);
     } else if (solver->computingASA()) {
         /* store x to compute jump in discontinuity */
         x_disc_.push_back(x_);
@@ -197,7 +210,8 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
         xdot_old_disc_.push_back(xdot_old_);
     }
 
-    model->updateHeaviside(roots_found_);
+    if (!initial_event)
+        model->updateHeaviside(roots_found_);
 
     applyEventBolus();
 
@@ -239,7 +253,7 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) {
                                  "Secondary event was triggered. Depending on "
                                  "the bolus of the secondary event, forward "
                                  "sensitivities can be incorrect.");
-        handleEvent(tlastroot, true);
+        handleEvent(tlastroot, true, false);
     }
 
     /* only reinitialise in the first event fired */
diff --git a/src/misc.cpp b/src/misc.cpp
index 60e9d3ea18..9b4c650cba 100644
--- a/src/misc.cpp
+++ b/src/misc.cpp
@@ -178,4 +178,9 @@ std::string printfToString(const char *fmt, va_list ap) {
     return str;
 }
 
+std::pair<size_t, size_t> unravel_index(size_t flat_idx, size_t num_cols)
+{
+    return {flat_idx / num_cols, flat_idx % num_cols};
+}
+
 } // namespace amici
diff --git a/src/model.ODE_template.cpp b/src/model.ODE_template.cpp
index c5fb93198d..049124881d 100644
--- a/src/model.ODE_template.cpp
+++ b/src/model.ODE_template.cpp
@@ -53,6 +53,10 @@ std::array<int, TPL_NX_SOLVER> stateIdxsSolver = {
     TPL_STATE_IDXS_SOLVER_INITIALIZER_LIST
 };
 
+std::array<bool, TPL_NEVENT> rootInitialValues = {
+    TPL_ROOT_INITIAL_VALUES
+};
+
 } // namespace model_TPL_MODELNAME
 
 } // namespace amici
diff --git a/src/model.cpp b/src/model.cpp
index 1edc211dde..e6097fa88f 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -15,6 +15,39 @@
 
 namespace amici {
 
+/**
+ * @brief Maps ModelQuantity items to their string value
+ */
+const std::map<ModelQuantity, std::string> model_quantity_to_str {
+    {ModelQuantity::J, "J"},
+    {ModelQuantity::JB, "JB"},
+    {ModelQuantity::Jv, "Jv"},
+    {ModelQuantity::JvB, "JvB"},
+    {ModelQuantity::JDiag, "JDiag"},
+    {ModelQuantity::sx, "sx"},
+    {ModelQuantity::sy, "sy"},
+    {ModelQuantity::ssigmay, "ssigmay"},
+    {ModelQuantity::xdot, "xdot"},
+    {ModelQuantity::sxdot, "sxdot"},
+    {ModelQuantity::xBdot, "xBdot"},
+    {ModelQuantity::x0, "x0"},
+    {ModelQuantity::x0_rdata, "x0_rdata"},
+    {ModelQuantity::x, "x"},
+    {ModelQuantity::x_rdata, "x_rdata"},
+    {ModelQuantity::dwdw, "dwdw"},
+    {ModelQuantity::dwdx, "dwdx"},
+    {ModelQuantity::dwdp, "dwdp"},
+    {ModelQuantity::y, "y"},
+    {ModelQuantity::dydp, "dydp"},
+    {ModelQuantity::dydx, "dydx"},
+    {ModelQuantity::w, "w"},
+    {ModelQuantity::root, "root"},
+    {ModelQuantity::qBdot, "qBdot"},
+    {ModelQuantity::qBdot_ss, "qBdot_ss"},
+    {ModelQuantity::xBdot_ss, "xBdot_ss"},
+    {ModelQuantity::JSparseB_ss, "JSparseB_ss"},
+};
+
 /**
  * @brief local helper function to get parameters
  * @param ids vector of name/ids of (fixed)Parameters
@@ -97,7 +130,8 @@ static int setValueByIdRegex(std::vector<std::string> const &ids,
 
 Model::Model(ModelDimensions const & model_dimensions,
              SimulationParameters simulation_parameters,
-             SecondOrderMode o2mode, std::vector<realtype> idlist, std::vector<int> z2event,
+             SecondOrderMode o2mode, std::vector<realtype> idlist,
+             std::vector<int> z2event,
              const bool pythonGenerated, const int ndxdotdp_explicit,
              const int ndxdotdx_explicit, const int w_recursion_depth)
     : ModelDimensions(model_dimensions), pythonGenerated(pythonGenerated),
@@ -120,6 +154,8 @@ Model::Model(ModelDimensions const & model_dimensions,
                       simulation_parameters_.pscale, state_.unscaledParameters);
     state_.fixedParameters = simulation_parameters_.fixedParameters;
     state_.plist = simulation_parameters_.plist;
+          
+    root_initial_values_.resize(ne, true);
 
     /* If Matlab wrapped: dxdotdp is a full AmiVector,
        if Python wrapped: dxdotdp_explicit and dxdotdp_implicit are CSC matrices
@@ -208,7 +244,8 @@ bool operator==(const ModelDimensions &a, const ModelDimensions &b) {
 
 
 void Model::initialize(AmiVector &x, AmiVector &dx, AmiVectorArray &sx,
-                       AmiVectorArray & /*sdx*/, bool computeSensitivities) {
+                       AmiVectorArray & /*sdx*/, bool computeSensitivities,
+                       std::vector<int> &roots_found) {
     initializeStates(x);
     if (computeSensitivities)
         initializeStateSensitivities(sx, x);
@@ -218,7 +255,7 @@ void Model::initialize(AmiVector &x, AmiVector &dx, AmiVectorArray &sx,
         fsdx0();
 
     if (ne)
-        initHeaviside(x, dx);
+        initEvents(x, dx, roots_found);
 }
 
 void Model::initializeB(AmiVector &xB, AmiVector &dxB, AmiVector &xQB,
@@ -265,14 +302,18 @@ void Model::initializeStateSensitivities(AmiVectorArray &sx,
     }
 }
 
-void Model::initHeaviside(AmiVector const &x, AmiVector const &dx) {
+void Model::initEvents(AmiVector const &x, AmiVector const &dx,
+                       std::vector<int> &roots_found) {
     std::vector<realtype> rootvals(ne, 0.0);
     froot(simulation_parameters_.tstart_, x, dx, rootvals);
+    std::fill(roots_found.begin(), roots_found.end(), 0);
     for (int ie = 0; ie < ne; ie++) {
         if (rootvals.at(ie) < 0) {
             state_.h.at(ie) = 0.0;
         } else {
             state_.h.at(ie) = 1.0;
+            if (!root_initial_values_.at(ie)) // only false->true triggers event
+                roots_found.at(ie) = 1;
         }
     }
 }
@@ -843,7 +884,7 @@ void Model::getObservableSensitivity(gsl::span<realtype> sy, const realtype t,
     writeSlice(derived_state_.dydp_, sy);
 
     if (always_check_finite_)
-        checkFinite(sy, "sy");
+        checkFinite(sy, ModelQuantity::sy, nplist());
 }
 
 void Model::getObservableSigma(gsl::span<realtype> sigmay, const int it,
@@ -875,7 +916,7 @@ void Model::getObservableSigmaSensitivity(gsl::span<realtype> ssigmay,
     }
 
     if (always_check_finite_)
-        checkFinite(ssigmay, "ssigmay");
+        checkFinite(ssigmay, ModelQuantity::ssigmay, nplist());
 }
 
 void Model::addObservableObjective(realtype &Jy, const int it,
@@ -1218,18 +1259,242 @@ void Model::updateHeavisideB(const int *rootsfound) {
     }
 }
 
+int Model::checkFinite(gsl::span<const realtype> array,
+                       ModelQuantity model_quantity) const
+{
+    auto it = std::find_if(array.begin(), array.end(),
+                           [](realtype x){return !std::isfinite(x);});
+    if(it == array.end()) {
+        return AMICI_SUCCESS;
+    }
+
+    // there is some issue - produce a meaningful message
+    auto flat_index = it - array.begin();
+
+    std::string msg_id;
+    std::string non_finite_type;
+    if (std::isnan(array[flat_index])) {
+        msg_id = "AMICI:NaN";
+        non_finite_type = "NaN";
+    } else if (std::isinf(array[flat_index])) {
+        msg_id = "AMICI:Inf";
+        non_finite_type = "Inf";
+    }
+    std::string element_id = std::to_string(flat_index);
+
+    switch (model_quantity) {
+    case ModelQuantity::xdot:
+    case ModelQuantity::xBdot:
+    case ModelQuantity::x0:
+    case ModelQuantity::x:
+    case ModelQuantity::x_rdata:
+    case ModelQuantity::x0_rdata:
+    case ModelQuantity::Jv:
+    case ModelQuantity::JvB:
+    case ModelQuantity::JDiag:
+        if(hasStateIds()) {
+            element_id = getStateIdsSolver()[flat_index];
+        }
+        break;
+    case ModelQuantity::y:
+        if(hasObservableIds()) {
+            element_id = getObservableIds()[flat_index];
+        }
+        break;
+    case ModelQuantity::w:
+        if(hasExpressionIds()) {
+            element_id = getExpressionIds()[flat_index];
+        }
+        break;
+    default:
+        break;
+    }
 
-int Model::checkFinite(gsl::span<const realtype> array, const char *fun) const {
-    auto result = app->checkFinite(array, fun);
-
-    if (result != AMICI_SUCCESS) {
-        app->checkFinite(state_.fixedParameters, "k");
-        app->checkFinite(state_.unscaledParameters, "p");
+    std::string model_quantity_str;
+    try {
+        model_quantity_str = model_quantity_to_str.at(model_quantity);
+    } catch (std::out_of_range const&) {
+        // Missing model quantity string - terminate if this is a debug build,
+        // but show the quantity number if non-debug.
+        gsl_ExpectsDebug(false);
+        model_quantity_str = std::to_string(static_cast<int>(model_quantity));
+    }
+    app->warningF(msg_id.c_str(),
+                  "AMICI encountered a %s value for %s[%i] (%s)",
+                  non_finite_type.c_str(),
+                  model_quantity_str.c_str(),
+                  static_cast<int>(flat_index),
+                  element_id.c_str()
+                  );
+
+    // check upstream
+    app->checkFinite(state_.fixedParameters, "k");
+    app->checkFinite(state_.unscaledParameters, "p");
+    app->checkFinite(simulation_parameters_.ts_, "t");
+    if(!always_check_finite_) {
+        // don't check twice if always_check_finite_ is true
         app->checkFinite(derived_state_.w_, "w");
-        app->checkFinite(simulation_parameters_.ts_, "t");
     }
 
-    return result;
+    return AMICI_RECOVERABLE_ERROR;
+}
+
+int Model::checkFinite(gsl::span<const realtype> array,
+                       ModelQuantity model_quantity, size_t num_cols) const
+{
+    auto it = std::find_if(array.begin(), array.end(),
+                           [](realtype x){return !std::isfinite(x);});
+    if(it == array.end()) {
+        return AMICI_SUCCESS;
+    }
+
+    // there is some issue - produce a meaningful message
+    auto flat_index = it - array.begin();
+    sunindextype row, col;
+    std::tie(row, col) = unravel_index(flat_index, num_cols);
+    std::string msg_id;
+    std::string non_finite_type;
+    if (std::isnan(array[flat_index])) {
+        msg_id = "AMICI:NaN";
+        non_finite_type = "NaN";
+    } else if (std::isinf(array[flat_index])) {
+        msg_id = "AMICI:Inf";
+        non_finite_type = "Inf";
+    }
+    std::string row_id = std::to_string(row);
+    std::string col_id = std::to_string(col);
+
+    switch (model_quantity) {
+    case ModelQuantity::sy:
+    case ModelQuantity::ssigmay:
+    case ModelQuantity::dydp:
+        row_id += " " + getObservableIds()[row];
+        col_id += " " + getParameterIds()[plist(col)];
+        break;
+    case ModelQuantity::dydx:
+        row_id += " " + getObservableIds()[row];
+        col_id += " " + getStateIdsSolver()[col];
+        break;
+    default:
+        break;
+    }
+
+    std::string model_quantity_str;
+    try {
+        model_quantity_str = model_quantity_to_str.at(model_quantity);
+    } catch (std::out_of_range const&) {
+        // Missing model quantity string - terminate if this is a debug build,
+        // but show the quantity number if non-debug.
+        gsl_ExpectsDebug(false);
+        model_quantity_str = std::to_string(static_cast<int>(model_quantity));
+    }
+
+    app->warningF(msg_id.c_str(),
+                  "AMICI encountered a %s value for %s[%i] (%s, %s)",
+                  non_finite_type.c_str(),
+                  model_quantity_str.c_str(),
+                  static_cast<int>(flat_index),
+                  row_id.c_str(),
+                  col_id.c_str()
+                  );
+
+    // check upstream
+    app->checkFinite(state_.fixedParameters, "k");
+    app->checkFinite(state_.unscaledParameters, "p");
+    app->checkFinite(simulation_parameters_.ts_, "t");
+    app->checkFinite(derived_state_.w_, "w");
+
+    return AMICI_RECOVERABLE_ERROR;
+}
+
+int Model::checkFinite(SUNMatrix m, ModelQuantity model_quantity, realtype t) const
+{
+    // check flat array, to see if there are any issues
+    // (faster, in particular for sparse arrays)
+    auto m_flat = gsl::make_span(m);
+    auto it = std::find_if(m_flat.begin(), m_flat.end(),
+                    [](realtype x){return !std::isfinite(x);});
+    if(it == m_flat.end()) {
+        return AMICI_SUCCESS;
+    }
+
+    // there is some issue - produce a meaningful message
+    auto flat_index = it - m_flat.begin();
+    sunindextype row, col;
+    std::tie(row, col) = unravel_index(flat_index, m);
+    std::string msg_id;
+    std::string non_finite_type;
+    if (std::isnan(m_flat[flat_index])) {
+        msg_id = "AMICI:NaN";
+        non_finite_type = "NaN";
+    } else if (std::isinf(m_flat[flat_index])) {
+        msg_id = "AMICI:Inf";
+        non_finite_type = "Inf";
+    } else {
+        throw std::runtime_error(
+            "Value is not finite, but neither infinite nor NaN.");
+    }
+    std::string row_id = std::to_string(row);
+    std::string col_id = std::to_string(col);
+
+    switch (model_quantity) {
+    case ModelQuantity::J:
+    case ModelQuantity::JB:
+        if(hasStateIds()) {
+            auto state_ids = getStateIdsSolver();
+            row_id += " " + state_ids[row];
+            col_id += " " + state_ids[col];
+        }
+        break;
+    case ModelQuantity::dwdx:
+        if(hasExpressionIds())
+            row_id += " " + getExpressionIds()[row];
+        if(hasStateIds())
+            col_id += " " + getStateIdsSolver()[col];
+        break;
+    case ModelQuantity::dwdw:
+        if(hasExpressionIds()) {
+            auto expr_ids = getExpressionIds();
+            row_id += " " + expr_ids[row];
+            col_id += " " + expr_ids[col];
+        }
+        break;
+    case ModelQuantity::dwdp:
+        if(hasExpressionIds())
+            row_id += " " + getExpressionIds()[row];
+        if(hasParameterIds())
+            col_id += " " + getParameterIds()[plist(col)];
+        break;
+    default:
+        break;
+    }
+
+    std::string model_quantity_str;
+    try {
+        model_quantity_str = model_quantity_to_str.at(model_quantity);
+    } catch (std::out_of_range const&) {
+        // Missing model quantity string - terminate if this is a debug build,
+        // but show the quantity number if non-debug.
+        gsl_ExpectsDebug(false);
+        model_quantity_str = std::to_string(static_cast<int>(model_quantity));
+    }
+    app->warningF(msg_id.c_str(),
+                  "AMICI encountered a %s value for %s[%i] (%s, %s) at t=%g",
+                  non_finite_type.c_str(),
+                  model_quantity_str.c_str(),
+                  static_cast<int>(flat_index),
+                  row_id.c_str(),
+                  col_id.c_str(),
+                  t
+                  );
+
+    // check upstream
+    app->checkFinite(state_.fixedParameters, "k");
+    app->checkFinite(state_.unscaledParameters, "p");
+    app->checkFinite(simulation_parameters_.ts_, "t");
+    app->checkFinite(derived_state_.w_, "w");
+
+    return AMICI_RECOVERABLE_ERROR;
 }
 
 void Model::setAlwaysCheckFinite(bool alwaysCheck) {
@@ -1250,8 +1515,8 @@ void Model::fx0(AmiVector &x) {
               state_.fixedParameters.data());
 
     if (always_check_finite_) {
-        checkFinite(derived_state_.x_rdata_, "x0 x_rdata");
-        checkFinite(x.getVector(), "x0 x");
+        checkFinite(derived_state_.x_rdata_, ModelQuantity::x0_rdata);
+        checkFinite(x.getVector(), ModelQuantity::x0);
     }
 }
 
@@ -1333,7 +1598,7 @@ void Model::fx_rdata(AmiVector &x_rdata, const AmiVector &x) {
     fx_rdata(x_rdata.data(), computeX_pos(x), state_.total_cl.data(),
              state_.unscaledParameters.data(), state_.fixedParameters.data());
     if (always_check_finite_)
-        checkFinite(x_rdata.getVector(), "x_rdata");
+        checkFinite(x_rdata.getVector(), ModelQuantity::x_rdata);
 }
 
 void Model::fsx_rdata(AmiVectorArray &sx_rdata, const AmiVectorArray &sx,
@@ -1411,7 +1676,7 @@ void Model::fy(const realtype t, const AmiVector &x) {
        state_.h.data(), derived_state_.w_.data());
 
     if (always_check_finite_) {
-        app->checkFinite(gsl::make_span(derived_state_.y_.data(), ny), "y");
+        checkFinite(gsl::make_span(derived_state_.y_.data(), ny), ModelQuantity::y);
     }
 }
 
@@ -1441,7 +1706,7 @@ void Model::fdydp(const realtype t, const AmiVector &x) {
         }
 
     if (always_check_finite_) {
-        app->checkFinite(derived_state_.dydp_, "dydp");
+        checkFinite(derived_state_.dydp_, ModelQuantity::dydp, nplist());
     }
 }
 
@@ -1461,7 +1726,7 @@ void Model::fdydx(const realtype t, const AmiVector &x) {
           derived_state_.dwdx_.data());
 
     if (always_check_finite_) {
-        app->checkFinite(derived_state_.dydx_, "dydx");
+        checkFinite(derived_state_.dydx_, ModelQuantity::dydx, ny);
     }
 }
 
@@ -2076,7 +2341,7 @@ void Model::fw(const realtype t, const realtype *x) {
        state_.fixedParameters.data(), state_.h.data(), state_.total_cl.data());
 
     if (always_check_finite_) {
-        app->checkFinite(derived_state_.w_, "w");
+        checkFinite(derived_state_.w_, ModelQuantity::w);
     }
 }
 
@@ -2116,7 +2381,7 @@ void Model::fdwdp(const realtype t, const realtype *x) {
     }
 
     if (always_check_finite_) {
-        app->checkFinite(gsl::make_span(derived_state_.dwdp_.get()), "dwdp");
+        checkFinite(derived_state_.dwdp_.get(), ModelQuantity::dwdp, t);
     }
 }
 
@@ -2157,7 +2422,7 @@ void Model::fdwdx(const realtype t, const realtype *x) {
     }
 
     if (always_check_finite_) {
-        app->checkFinite(gsl::make_span(derived_state_.dwdx_.get()), "dwdx");
+        checkFinite(derived_state_.dwdx_.get(), ModelQuantity::dwdx, t);
     }
 }
 
@@ -2172,7 +2437,7 @@ void Model::fdwdw(const realtype t, const realtype *x) {
           derived_state_.w_.data(), state_.total_cl.data());
 
     if (always_check_finite_) {
-        app->checkFinite(gsl::make_span(dwdw_.get()), "dwdw");
+        checkFinite(dwdw_.get(), ModelQuantity::dwdw, t);
     }
 }
 
diff --git a/src/model_dae.cpp b/src/model_dae.cpp
index bed9b4c8c9..2e3d6d06ee 100644
--- a/src/model_dae.cpp
+++ b/src/model_dae.cpp
@@ -87,7 +87,7 @@ void Model_DAE::fJDiag(const realtype t, AmiVector &JDiag,
     fJSparse(t, 0.0, x.getNVector(), dx.getNVector(), derived_state_.J_.get());
     derived_state_.J_.refresh();
     derived_state_.J_.to_diag(JDiag.getNVector());
-    if (!checkFinite(JDiag.getVector(), "Jacobian"))
+    if (checkFinite(JDiag.getVector(), ModelQuantity::JDiag) != AMICI_SUCCESS)
         throw AmiException("Evaluation of fJDiag failed!");
 }
 
diff --git a/src/model_header.ODE_template.h b/src/model_header.ODE_template.h
index 27e664bd8f..8a885a61ba 100644
--- a/src/model_header.ODE_template.h
+++ b/src/model_header.ODE_template.h
@@ -27,6 +27,7 @@ extern std::array<const char*, TPL_NX_RDATA> stateIds;
 extern std::array<const char*, TPL_NY> observableIds;
 extern std::array<const char*, TPL_NW> expressionIds;
 extern std::array<int, TPL_NX_SOLVER> stateIdxsSolver;
+extern std::array<bool, TPL_NEVENT> rootInitialValues;
 
 TPL_JY_DEF
 TPL_DJYDSIGMA_DEF
@@ -129,7 +130,11 @@ class Model_TPL_MODELNAME : public amici::Model_ODE {
               TPL_NDXDOTDP_EXPLICIT,                       // ndxdotdp_explicit
               TPL_NDXDOTDX_EXPLICIT,                       // ndxdotdx_explicit
               TPL_W_RECURSION_DEPTH                        // w_recursion_depth
-          ) {}
+          ) {
+                 root_initial_values_ = std::vector<bool>(
+                     rootInitialValues.begin(), rootInitialValues.end()
+                 );
+          }
 
     /**
      * @brief Clone this model instance.
diff --git a/src/model_ode.cpp b/src/model_ode.cpp
index 02905ef2c0..65552b8579 100644
--- a/src/model_ode.cpp
+++ b/src/model_ode.cpp
@@ -105,7 +105,7 @@ void Model_ODE::fJDiag(const realtype t, AmiVector &JDiag,
                        const realtype /*cj*/, const AmiVector &x,
                        const AmiVector & /*dx*/) {
     fJDiag(t, JDiag.getNVector(), x.getNVector());
-    if (checkFinite(JDiag.getVector(), "Jacobian") != AMICI_SUCCESS)
+    if (checkFinite(JDiag.getVector(), ModelQuantity::JDiag) != AMICI_SUCCESS)
         throw AmiException("Evaluation of fJDiag failed!");
 }
 
diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp
index c37b1207cd..74998fdf66 100644
--- a/src/solver_cvodes.cpp
+++ b/src/solver_cvodes.cpp
@@ -807,7 +807,6 @@ const Model *CVodeSolver::getModel() const {
 
 /**
  * @brief Jacobian of xdot with respect to states x
- * @param N number of state variables
  * @param t timepoint
  * @param x Vector with the states
  * @param xdot Vector with the right hand side
@@ -827,13 +826,12 @@ static int fJ(realtype t, N_Vector x, N_Vector xdot, SUNMatrix J,
     Expects(model);
 
     model->fJ(t, x, xdot, J);
-    return model->checkFinite(gsl::make_span(J), "Jacobian");
+    return model->checkFinite(J, ModelQuantity::J, t);
 }
 
 
 /**
  * @brief Jacobian of xBdot with respect to adjoint state xB
- * @param NeqBdot number of adjoint state variables
  * @param t timepoint
  * @param x Vector with the states
  * @param xB Vector with the adjoint states
@@ -854,7 +852,7 @@ static int fJB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot,
     Expects(model);
 
     model->fJB(t, x, xB, xBdot, JB);
-    return model->checkFinite(gsl::make_span(JB), "Jacobian");
+    return model->checkFinite(gsl::make_span(JB), ModelQuantity::JB);
 }
 
 
@@ -879,7 +877,7 @@ static int fJSparse(realtype t, N_Vector x, N_Vector /*xdot*/,
     Expects(model);
 
     model->fJSparse(t, x, J);
-    return model->checkFinite(gsl::make_span(J), "Jacobian");
+    return model->checkFinite(J, ModelQuantity::J, t);
 }
 
 
@@ -905,7 +903,7 @@ static int fJSparseB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot,
     Expects(model);
 
     model->fJSparseB(t, x, xB, xBdot, JB);
-    return model->checkFinite(gsl::make_span(JB), "Jacobian");
+    return model->checkFinite(gsl::make_span(JB), ModelQuantity::JB);
 }
 
 
@@ -967,7 +965,7 @@ static int fJv(N_Vector v, N_Vector Jv, realtype t, N_Vector x,
     Expects(model);
 
     model->fJv(v, Jv, t, x);
-    return model->checkFinite(gsl::make_span(Jv), "Jacobian");
+    return model->checkFinite(gsl::make_span(Jv), ModelQuantity::Jv);
 }
 
 
@@ -993,7 +991,7 @@ static int fJvB(N_Vector vB, N_Vector JvB, realtype t, N_Vector x,
     Expects(model);
 
     model->fJvB(vB, JvB, t, x, xB);
-    return model->checkFinite(gsl::make_span(JvB), "Jacobian");
+    return model->checkFinite(gsl::make_span(JvB), ModelQuantity::JvB);
 }
 
 
@@ -1014,7 +1012,7 @@ static int froot(realtype t, N_Vector x, realtype *root,
 
     model->froot(t, x, gsl::make_span<realtype>(root, model->ne));
     return model->checkFinite(gsl::make_span<realtype>(root, model->ne),
-                              "root function");
+                              ModelQuantity::root);
 }
 
 
@@ -1038,7 +1036,8 @@ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data) {
         return AMICI_MAX_TIME_EXCEEDED;
     }
 
-    if (t > 1e200 && !model->checkFinite(gsl::make_span(x), "fxdot")) {
+    if (t > 1e200
+        && !model->checkFinite(gsl::make_span(x), ModelQuantity::xdot)) {
         /* when t is large (typically ~1e300), CVODES may pass all NaN x
            to fxdot from which we typically cannot recover. To save time
            on normal execution, we do not always want to check finiteness
@@ -1047,7 +1046,7 @@ static int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data) {
     }
 
     model->fxdot(t, x, xdot);
-    return model->checkFinite(gsl::make_span(xdot), "fxdot");
+    return model->checkFinite(gsl::make_span(xdot), ModelQuantity::xdot);
 }
 
 
@@ -1074,7 +1073,7 @@ static int fxBdot(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot,
     }
 
     model->fxBdot(t, x, xB, xBdot);
-    return model->checkFinite(gsl::make_span(xBdot), "fxBdot");
+    return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot);
 }
 
 
@@ -1095,7 +1094,7 @@ static int fqBdot(realtype t, N_Vector x, N_Vector xB, N_Vector qBdot,
     Expects(model);
 
     model->fqBdot(t, x, xB, qBdot);
-    return model->checkFinite(gsl::make_span(qBdot), "qBdot");
+    return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot);
 }
 
 
@@ -1116,7 +1115,7 @@ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector xBdot,
     Expects(model);
 
     model->fxBdot_ss(t, xB, xBdot);
-    return model->checkFinite(gsl::make_span(xBdot), "fxBdot_ss");
+    return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot_ss);
 }
 
 
@@ -1137,7 +1136,7 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector qBdot,
     Expects(model);
 
     model->fqBdot_ss(t, xB, qBdot);
-    return model->checkFinite(gsl::make_span(qBdot), "qBdot_ss");
+    return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot_ss);
 }
 
 /**
@@ -1161,7 +1160,8 @@ static int fJSparseB_ss(realtype /*t*/, N_Vector /*x*/, N_Vector xBdot,
     Expects(model);
 
     model->fJSparseB_ss(JB);
-    return model->checkFinite(gsl::make_span(xBdot), "JSparseB_ss");
+    return model->checkFinite(gsl::make_span(xBdot),
+                              ModelQuantity::JSparseB_ss);
 }
 
 
@@ -1189,7 +1189,7 @@ static int fsxdot(int /*Ns*/, realtype t, N_Vector x, N_Vector /*xdot*/,
     Expects(model);
 
     model->fsxdot(t, x, ip, sx, sxdot);
-    return model->checkFinite(gsl::make_span(sxdot), "sxdot");
+    return model->checkFinite(gsl::make_span(sxdot), ModelQuantity::sxdot);
 }
 
 bool operator==(const CVodeSolver &a, const CVodeSolver &b) {
diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp
index 0742f9c095..0039f6de0a 100644
--- a/src/solver_idas.cpp
+++ b/src/solver_idas.cpp
@@ -817,7 +817,7 @@ int fJ(realtype t, realtype cj, N_Vector x, N_Vector dx,
     auto model = dynamic_cast<Model_DAE *>(typed_udata->first);
     Expects(model);
     model->fJ(t, cj, x, dx, xdot, J);
-    return model->checkFinite(gsl::make_span(J), "Jacobian");
+    return model->checkFinite(J, ModelQuantity::J, t);
 }
 
 /**
@@ -846,7 +846,8 @@ int fJB(realtype t, realtype cj, N_Vector x, N_Vector dx,
     Expects(model);
 
     model->fJB(t, cj, x, dx, xB, dxB, JB);
-    return model->checkFinite(gsl::make_span(JB), "Jacobian");
+    return model->checkFinite(JB, ModelQuantity::JB, t);
+
 }
 
 /**
@@ -873,7 +874,7 @@ int fJSparse(realtype t, realtype cj, N_Vector x, N_Vector dx,
     Expects(model);
 
     model->fJSparse(t, cj, x, dx, J);
-    return model->checkFinite(gsl::make_span(J), "Jacobian");
+    return model->checkFinite(J, ModelQuantity::J, t);
 }
 
 /**
@@ -902,7 +903,7 @@ int fJSparseB(realtype t, realtype cj, N_Vector x, N_Vector dx,
     Expects(model);
 
     model->fJSparseB(t, cj, x, dx, xB, dxB, JB);
-    return model->checkFinite(gsl::make_span(JB), "Jacobian");
+    return model->checkFinite(JB, ModelQuantity::JB, t);
 }
 
 /**
@@ -973,7 +974,7 @@ int fJv(realtype t, N_Vector x, N_Vector dx, N_Vector /*xdot*/,
     Expects(model);
 
     model->fJv(t, x, dx, v, Jv, cj);
-    return model->checkFinite(gsl::make_span(Jv), "Jacobian");
+    return model->checkFinite(gsl::make_span(Jv), ModelQuantity::Jv);
 }
 
 /**
@@ -1003,7 +1004,7 @@ int fJvB(realtype t, N_Vector x, N_Vector dx, N_Vector xB,
     Expects(model);
 
     model->fJvB(t, x, dx, xB, dxB, vB, JvB, cj);
-    return model->checkFinite(gsl::make_span(JvB), "Jacobian");
+    return model->checkFinite(gsl::make_span(JvB), ModelQuantity::JvB);
 }
 
 /**
@@ -1024,7 +1025,7 @@ int froot(realtype t, N_Vector x, N_Vector dx, realtype *root,
 
     model->froot(t, x, dx, gsl::make_span<realtype>(root, model->ne));
     return model->checkFinite(gsl::make_span<realtype>(root, model->ne),
-                              "root function");
+                              ModelQuantity::root);
 }
 
 /**
@@ -1058,7 +1059,7 @@ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot,
     }
 
     model->fxdot(t, x, dx, xdot);
-    return model->checkFinite(gsl::make_span(xdot), "fxdot");
+    return model->checkFinite(gsl::make_span(xdot), ModelQuantity::xdot);
 }
 
 /**
@@ -1086,7 +1087,7 @@ int fxBdot(realtype t, N_Vector x, N_Vector dx, N_Vector xB,
     }
 
     model->fxBdot(t, x, dx, xB, dxB, xBdot);
-    return model->checkFinite(gsl::make_span(xBdot), "xBdot");
+    return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot);
 }
 
 /**
@@ -1109,7 +1110,7 @@ int fqBdot(realtype t, N_Vector x, N_Vector dx, N_Vector xB,
     Expects(model);
 
     model->fqBdot(t, x, dx, xB, dxB, qBdot);
-    return model->checkFinite(gsl::make_span(qBdot), "qBdot");
+    return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot);
 
 }
 
@@ -1132,7 +1133,7 @@ static int fxBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector xBdot,
     Expects(model);
 
     model->fxBdot_ss(t, xB, dxB, xBdot);
-    return model->checkFinite(gsl::make_span(xBdot), "xBdot_ss");
+    return model->checkFinite(gsl::make_span(xBdot), ModelQuantity::xBdot_ss);
 }
 
 
@@ -1154,7 +1155,7 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector qBdot,
     Expects(model);
 
     model->fqBdot_ss(t, xB, dxB, qBdot);
-    return model->checkFinite(gsl::make_span(qBdot), "qBdot_ss");
+    return model->checkFinite(gsl::make_span(qBdot), ModelQuantity::qBdot_ss);
 }
 
 /**
@@ -1181,7 +1182,8 @@ static int fqBdot_ss(realtype t, N_Vector xB, N_Vector dxB, N_Vector qBdot,
     Expects(model);
 
     model->fJSparseB_ss(JB);
-    return model->checkFinite(gsl::make_span(xBdot), "JSparseB_ss");
+    return model->checkFinite(gsl::make_span(xBdot),
+                              ModelQuantity::JSparseB_ss);
 }
 
 /**
@@ -1212,7 +1214,7 @@ int fsxdot(int /*Ns*/, realtype t, N_Vector x, N_Vector dx,
 
     for (int ip = 0; ip < model->nplist(); ip++) {
         model->fsxdot(t, x, dx, ip, sx[ip], sdx[ip], sxdot[ip]);
-        if (model->checkFinite(gsl::make_span(sxdot[ip]), "sxdot")
+        if (model->checkFinite(gsl::make_span(sxdot[ip]), ModelQuantity::sxdot)
                 != AMICI_SUCCESS)
             return AMICI_RECOVERABLE_ERROR;
     }
diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp
index df5104dd9e..a2acec9596 100644
--- a/src/steadystateproblem.cpp
+++ b/src/steadystateproblem.cpp
@@ -188,9 +188,11 @@ void SteadystateProblem::initializeForwardProblem(int it, const Solver &solver,
     /* process solver handling for pre- or postequilibration */
     if (it == -1) {
         /* solver was not run before, set up everything */
+        auto roots_found = std::vector<int>(model.ne, 0);
         model.initialize(state_.x, state_.dx, state_.sx, sdx_,
-                          solver.getSensitivityOrder() >=
-                              SensitivityOrder::first);
+                         solver.getSensitivityOrder() >=
+                         SensitivityOrder::first,
+                         roots_found);
         state_.t = model.t0();
         solver.setup(state_.t, &model, state_.x, state_.dx, state_.sx, sdx_);
     } else {
@@ -629,7 +631,8 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver,
 
     /* If run after Newton's method checks again if it converged */
     wrms_ = getWrms(model, sensitivityFlag);
-    int sim_steps = 0;
+    
+    int &sim_steps = backward ? numstepsB_ : numsteps_.at(1);
     
     int convergence_check_frequency = 1;
     
@@ -653,22 +656,16 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver,
             flagUpdatedState();
         }
         
-        try {
         /* Check for convergence */
-        wrms_ = getWrms(model, sensitivityFlag);
-        /* getWrms needs to be called before getWrmsFSA such that the linear
-         system is prepared for newton type convergence check */
-        if (wrms_ < conv_thresh && check_sensi_conv_ &&
-            sensitivityFlag == SensitivityMethod::forward &&
-            sim_steps % convergence_check_frequency == 0) {
-            updateSensiSimulation(solver);
-            wrms_ = getWrmsFSA(model);
-        }
-        
-        } catch (NewtonFailure const &) {
-            /* linear solves in getWrms failed */
-            numsteps_.at(1) = sim_steps;
-            throw;
+        if (sim_steps % convergence_check_frequency == 0) {
+            wrms_ = getWrms(model, sensitivityFlag);
+            /* getWrms needs to be called before getWrmsFSA such that the linear
+             system is prepared for newton type convergence check */
+            if (wrms_ < conv_thresh && check_sensi_conv_ &&
+                sensitivityFlag == SensitivityMethod::forward) {
+                updateSensiSimulation(solver);
+                wrms_ = getWrmsFSA(model);
+            }
         }
 
         if (wrms_ < conv_thresh)
@@ -676,12 +673,10 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver,
         /* increase counter, check for maxsteps */
         sim_steps++;
         if (sim_steps >= solver.getMaxSteps()) {
-            numsteps_.at(1) = sim_steps;
             throw NewtonFailure(AMICI_TOO_MUCH_WORK,
                                 "exceeded maximum number of steps");
         }
         if (state_.t >= 1e200) {
-            numsteps_.at(1) = sim_steps;
             throw NewtonFailure(AMICI_NO_STEADY_STATE,
                                 "simulated to late time"
                                 " point without convergence of RHS");
@@ -691,13 +686,6 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver,
     // if check_sensi_conv_ is deactivated, we still have to update sensis
     if (sensitivityFlag == SensitivityMethod::forward)
         updateSensiSimulation(solver);
-
-    /* store information about steps and sensitivities, if necessary */
-    if (backward) {
-        numstepsB_ = sim_steps;
-    } else {
-        numsteps_.at(1) = sim_steps;
-    }
 }
 
 std::unique_ptr<Solver>
diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp
index bd3fda4fb8..f3515c357a 100644
--- a/src/sundials_matrix_wrapper.cpp
+++ b/src/sundials_matrix_wrapper.cpp
@@ -747,5 +747,45 @@ void SUNMatrixWrapper::refresh() {
 
 SUNMatrix SUNMatrixWrapper::get() const { return matrix_; }
 
+std::pair<sunindextype, sunindextype> unravel_index(sunindextype i, SUNMatrix m)
+{
+    gsl_ExpectsDebug(i >= 0);
+    auto mat_id = SUNMatGetID(m);
+    if(mat_id == SUNMATRIX_DENSE) {
+        gsl_ExpectsDebug(i < SM_COLUMNS_D(m) * SM_ROWS_D(m));
+
+        auto num_rows = SM_ROWS_D(m);
+        // col-major
+        sunindextype row = i % num_rows;
+        sunindextype col = i / num_rows;
+
+        gsl_EnsuresDebug(row >= 0);
+        gsl_EnsuresDebug(row < SM_ROWS_D(m));
+        gsl_EnsuresDebug(col >= 0);
+        gsl_EnsuresDebug(col < SM_COLUMNS_D(m));
+
+        return {row, col};
+    }
+
+    if(mat_id == SUNMATRIX_SPARSE) {
+        gsl_ExpectsDebug(i < SM_NNZ_S(m));
+        sunindextype row = SM_INDEXVALS_S(m)[i];
+        sunindextype i_colptr = 0;
+        while(SM_INDEXPTRS_S(m)[i_colptr] < SM_NNZ_S(m)) {
+            if(SM_INDEXPTRS_S(m)[i_colptr + 1] > i) {
+                sunindextype col = i_colptr;
+                gsl_EnsuresDebug(row >= 0);
+                gsl_EnsuresDebug(row < SM_ROWS_S(m));
+                gsl_EnsuresDebug(col >= 0);
+                gsl_EnsuresDebug(col < SM_COLUMNS_S(m));
+                return {row, col};
+            }
+            ++i_colptr;
+        }
+    }
+
+    throw amici::AmiException("Unimplemented SUNMatrix type for unravel_index");
+}
+
 } // namespace amici
 
diff --git a/swig/amici.i b/swig/amici.i
index a754cd61b8..ccfbf65933 100644
--- a/swig/amici.i
+++ b/swig/amici.i
@@ -120,6 +120,7 @@ wrap_unique_ptr(ExpDataPtr, amici::ExpData)
 %ignore amici::ContextManager;
 %ignore amici::ModelState;
 %ignore amici::ModelStateDerived;
+%ignore amici::unravel_index;
 
 // Include before any other header which uses enums defined there
 %include "amici/defines.h"
diff --git a/swig/model.i b/swig/model.i
index 81466a5c3b..d2c9f2eabd 100644
--- a/swig/model.i
+++ b/swig/model.i
@@ -32,7 +32,7 @@ using namespace amici;
 %ignore getObservable;
 %ignore getObservableSensitivity;
 %ignore getExpression;
-%ignore initHeaviside;
+%ignore initEvents;
 %ignore initialize;
 %ignore initializeB;
 %ignore initializeStateSensitivities;
diff --git a/tests/cpp/unittests/testMisc.cpp b/tests/cpp/unittests/testMisc.cpp
index 278e93b2bc..c689090c8c 100644
--- a/tests/cpp/unittests/testMisc.cpp
+++ b/tests/cpp/unittests/testMisc.cpp
@@ -6,6 +6,7 @@
 #include <amici/solver_cvodes.h>
 #include <amici/solver_idas.h>
 #include <amici/symbolic_functions.h>
+#include <amici/misc.h>
 
 #include <cmath>
 #include <cstring>
@@ -646,4 +647,62 @@ TEST_F(SunMatrixWrapperTest, BlockTranspose)
                   SM_INDEXPTRS_S(B_sparse.get())[icol]);
 }
 
+TEST(UnravelIndex, UnravelIndex)
+{
+    EXPECT_EQ(unravel_index(0, 2), std::make_pair((size_t) 0, (size_t) 0));
+    EXPECT_EQ(unravel_index(1, 2), std::make_pair((size_t) 0, (size_t) 1));
+    EXPECT_EQ(unravel_index(2, 2), std::make_pair((size_t) 1, (size_t) 0));
+    EXPECT_EQ(unravel_index(3, 2), std::make_pair((size_t) 1, (size_t) 1));
+    EXPECT_EQ(unravel_index(4, 2), std::make_pair((size_t) 2, (size_t) 0));
+    EXPECT_EQ(unravel_index(5, 2), std::make_pair((size_t) 2, (size_t) 1));
+    EXPECT_EQ(unravel_index(6, 2), std::make_pair((size_t) 3, (size_t) 0));
+}
+
+TEST(UnravelIndex, UnravelIndexSunMatDense)
+{
+    SUNMatrixWrapper A = SUNMatrixWrapper(3, 2);
+
+    A.set_data(0, 0, 0);
+    A.set_data(1, 0, 1);
+    A.set_data(2, 0, 2);
+    A.set_data(0, 1, 3);
+    A.set_data(1, 1, 4);
+    A.set_data(2, 1, 5);
+
+    for(int i = 0; i < 6; ++i) {
+        auto idx = unravel_index(i, A.get());
+        EXPECT_EQ(A.get_data(idx.first, idx.second), i);
+    }
+}
+
+TEST(UnravelIndex, UnravelIndexSunMatSparse)
+{
+    SUNMatrixWrapper D = SUNMatrixWrapper(4, 2);
+
+    // [0, 3]
+    // [0, 0]
+    // [1, 0]
+    // [2, 0]
+    // data [1, 2, 3]
+    // colptrs [0, 2, 3]
+    // rowidxs [2, 3, 1]
+    D.set_data(0, 0, 0);
+    D.set_data(1, 0, 0);
+    D.set_data(2, 0, 1);
+    D.set_data(3, 0, 2);
+    D.set_data(0, 1, 3);
+    D.set_data(1, 1, 0);
+    D.set_data(2, 1, 0);
+    D.set_data(3, 1, 0);
+
+    auto S = SUNSparseFromDenseMatrix(D.get(), 1e-15, CSC_MAT);
+
+    EXPECT_EQ(unravel_index(0, S), std::make_pair((sunindextype) 2, (sunindextype) 0));
+    EXPECT_EQ(unravel_index(1, S), std::make_pair((sunindextype) 3, (sunindextype) 0));
+    EXPECT_EQ(unravel_index(2, S), std::make_pair((sunindextype) 0, (sunindextype) 1));
+
+    SUNMatDestroy(S);
+}
+
+
 } // namespace
diff --git a/tests/testSBMLSuite.py b/tests/testSBMLSuite.py
index 3187720a8f..767c4287fb 100755
--- a/tests/testSBMLSuite.py
+++ b/tests/testSBMLSuite.py
@@ -60,6 +60,8 @@ def test_sbml_testsuite_case(
     sensitivity_check_cases = {
         # parameter-dependent conservation laws
         '00783': 1.5e-2,
+        # initial events
+        '00995': 1e-3,
     }
 
     try:
diff --git a/version.txt b/version.txt
index 2b71711fca..36e40e13e0 100644
--- a/version.txt
+++ b/version.txt
@@ -1 +1 @@
-0.11.28
+0.11.29