diff --git a/include/amici/event.h b/include/amici/event.h
index 7973396e52..af1d0a98ca 100644
--- a/include/amici/event.h
+++ b/include/amici/event.h
@@ -56,7 +56,7 @@ class Event {
      * @brief Get the initial value of the root function
      * @return The value of the root function at t_0.
      */
-    bool get_initial_value() { return initial_value_; }
+    bool get_initial_value() const { return initial_value_; }
 
     /**
      * @brief Get the priority of the event assignments
diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h
index 4d0604ebe4..b871328cd3 100644
--- a/include/amici/model_dae.h
+++ b/include/amici/model_dae.h
@@ -148,7 +148,7 @@ class Model_DAE : public Model {
 
     void
     fJv(realtype t, AmiVector const& x, AmiVector const& dx,
-        AmiVector const& xdot, AmiVector const& v, AmiVector& nJv,
+        AmiVector const& xdot, AmiVector const& v, AmiVector& Jv,
         realtype cj) override;
 
     /**
diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h
index 201d9ce5b0..f40014ed5e 100644
--- a/include/amici/model_ode.h
+++ b/include/amici/model_ode.h
@@ -141,7 +141,7 @@ class Model_ODE : public Model {
 
     void
     fJv(realtype t, AmiVector const& x, AmiVector const& dx,
-        AmiVector const& xdot, AmiVector const& v, AmiVector& nJv,
+        AmiVector const& xdot, AmiVector const& v, AmiVector& Jv,
         realtype cj) override;
 
     /**
diff --git a/include/amici/solver.h b/include/amici/solver.h
index 6714e61251..1640ccd0af 100644
--- a/include/amici/solver.h
+++ b/include/amici/solver.h
@@ -24,12 +24,10 @@ class Solver;
 } // namespace amici
 
 // for serialization friend in Solver
-namespace boost {
-namespace serialization {
+namespace boost::serialization {
 template <class Archive>
 void serialize(Archive& ar, amici::Solver& s, unsigned int version);
-}
-} // namespace boost
+} // namespace boost::serialization
 
 namespace amici {
 
diff --git a/include/amici/vector.h b/include/amici/vector.h
index 7775cf651b..e6092dd9a7 100644
--- a/include/amici/vector.h
+++ b/include/amici/vector.h
@@ -16,12 +16,10 @@ class AmiVector;
 }
 
 // for serialization friend
-namespace boost {
-namespace serialization {
+namespace boost::serialization {
 template <class Archive>
 void serialize(Archive& ar, amici::AmiVector& v, unsigned int version);
-}
-} // namespace boost
+} // namespace boost::serialization
 
 namespace amici {
 
@@ -51,7 +49,7 @@ class AmiVector {
      * @param length number of elements in vector
      * @param sunctx SUNDIALS context
      */
-    explicit AmiVector(long int const length, SUNContext sunctx)
+    explicit AmiVector(long int const length, SUNContext const sunctx)
         : vec_(static_cast<decltype(vec_)::size_type>(length), 0.0)
         , nvec_(N_VMake_Serial(length, vec_.data(), sunctx)) {}
 
@@ -63,7 +61,7 @@ class AmiVector {
      * @param rvec vector from which the data will be moved
      * @param sunctx SUNDIALS context
      */
-    explicit AmiVector(std::vector<realtype> rvec, SUNContext sunctx)
+    explicit AmiVector(std::vector<realtype> rvec, SUNContext const sunctx)
         : vec_(std::move(rvec))
         , nvec_(N_VMake_Serial(
               gsl::narrow<long int>(vec_.size()), vec_.data(), sunctx
@@ -74,7 +72,7 @@ class AmiVector {
      * @param rvec vector from which the data will be copied
      * @param sunctx SUNDIALS context
      */
-    explicit AmiVector(gsl::span<realtype const> rvec, SUNContext sunctx)
+    explicit AmiVector(gsl::span<realtype const> const rvec, SUNContext const sunctx)
         : AmiVector(std::vector(rvec.begin(), rvec.end()), sunctx) {}
 
     /**
diff --git a/src/backwardproblem.cpp b/src/backwardproblem.cpp
index 9454ff6d3f..0984e8d3a1 100644
--- a/src/backwardproblem.cpp
+++ b/src/backwardproblem.cpp
@@ -204,7 +204,7 @@ BwdSimWorkspace::BwdSimWorkspace(
     , xQB_(model_->nJ * model_->nplist(), solver->getSunContext()) {}
 
 void EventHandlingBwdSimulator::run(
-    realtype t_start, realtype t_end, realtype it,
+    realtype const t_start, realtype const t_end, realtype it,
     std::vector<realtype> const& timepoints, std::vector<realtype> const* dJydx,
     std::vector<realtype> const* dJzdx
 ) {
@@ -311,7 +311,7 @@ SteadyStateBackwardProblem::SteadyStateBackwardProblem(
     , solver_(&solver)
     , ws_(ws) {}
 
-void SteadyStateBackwardProblem::run(realtype t0) {
+void SteadyStateBackwardProblem::run(realtype const t0) {
     newton_solver_.reinitialize();
 
     // initialize quadratures
@@ -319,7 +319,7 @@ void SteadyStateBackwardProblem::run(realtype t0) {
     ws_->xQB_.zero();
 
     // Compute quadratures, track computation time
-    CpuTimer cpu_timer;
+    CpuTimer const cpu_timer;
 
     compute_steady_state_quadrature(t0);
     cpu_timeB_ = cpu_timer.elapsed_milliseconds();
@@ -333,7 +333,7 @@ AmiVector const& SteadyStateBackwardProblem::getAdjointQuadrature() const {
     return ws_->xQB_;
 }
 
-void SteadyStateBackwardProblem::compute_steady_state_quadrature(realtype t0) {
+void SteadyStateBackwardProblem::compute_steady_state_quadrature(realtype const t0) {
     // This routine computes the quadratures:
     //     xQB = Integral[ xB(x(t), t, p) * dxdot/dp(x(t), t, p) | dt ]
     // As we're in steady state, we have x(t) = x_ss (x_steadystate), hence
@@ -399,7 +399,7 @@ void SteadyStateBackwardProblem::compute_quadrature_by_lin_solve() {
     }
 }
 
-void SteadyStateBackwardProblem::compute_quadrature_by_simulation(realtype t0) {
+void SteadyStateBackwardProblem::compute_quadrature_by_simulation(realtype const t0) {
     // If the Jacobian is singular, the integral over xB must be computed
     // by usual integration over time, but simplifications can be applied:
     // x is not time-dependent, no forward trajectory is needed.
@@ -409,7 +409,7 @@ void SteadyStateBackwardProblem::compute_quadrature_by_simulation(realtype t0) {
     // xQ was written in getQuadratureByLinSolve() -> set to zero
     xQ_.zero();
 
-    auto sim_solver = std::unique_ptr<Solver>(solver_->clone());
+    auto const sim_solver = std::unique_ptr<Solver>(solver_->clone());
     sim_solver->logger = solver_->logger;
     sim_solver->setSensitivityMethod(SensitivityMethod::none);
     sim_solver->setSensitivityOrder(SensitivityOrder::none);
@@ -454,7 +454,7 @@ void SteadyStateBackwardProblem::run_simulation(Solver const& solver) {
     AmiVector xQBdot(model_->nplist(), solver.getSunContext());
 
     int const convergence_check_frequency = newton_step_conv_ ? 25 : 1;
-    auto max_steps = (solver.getMaxStepsBackwardProblem() > 0)
+    auto const max_steps = (solver.getMaxStepsBackwardProblem() > 0)
                          ? solver.getMaxStepsBackwardProblem()
                          : solver.getMaxSteps() * 100;
 
diff --git a/src/cblas.cpp b/src/cblas.cpp
index 6a47e59d21..b0b63d8e2f 100644
--- a/src/cblas.cpp
+++ b/src/cblas.cpp
@@ -55,7 +55,7 @@ void amici_dgemv(
 }
 
 void amici_daxpy(
-    int n, double alpha, double const* x, int const incx, double* y, int incy
+    int const n, double const alpha, double const* x, int const incx, double* y, int const incy
 ) {
     BLAS_FUNC(daxpy)(n, alpha, x, incx, y, incy);
 }
diff --git a/src/edata.cpp b/src/edata.cpp
index 2c279d7d85..b59a6e38c8 100644
--- a/src/edata.cpp
+++ b/src/edata.cpp
@@ -10,7 +10,7 @@
 
 namespace amici {
 
-ExpData::ExpData(int nytrue, int nztrue, int nmaxevent)
+ExpData::ExpData(int const nytrue, int const nztrue, int const nmaxevent)
     : nytrue_(nytrue)
     , nztrue_(nztrue)
     , nmaxevent_(nmaxevent) {
@@ -18,7 +18,7 @@ ExpData::ExpData(int nytrue, int nztrue, int nmaxevent)
 }
 
 ExpData::ExpData(
-    int nytrue, int nztrue, int nmaxevent, std::vector<realtype> ts
+    int const nytrue, int const nztrue, int const nmaxevent, std::vector<realtype> ts
 )
     : SimulationParameters(ts)
     , nytrue_(nytrue)
@@ -28,7 +28,7 @@ ExpData::ExpData(
 }
 
 ExpData::ExpData(
-    int nytrue, int nztrue, int nmaxevent, std::vector<realtype> ts,
+    int const nytrue, int const nztrue, int const nmaxevent, std::vector<realtype> ts,
     std::vector<realtype> fixedParameters
 )
     : SimulationParameters(ts)
@@ -40,7 +40,7 @@ ExpData::ExpData(
 }
 
 ExpData::ExpData(
-    int nytrue, int nztrue, int nmaxevent, std::vector<realtype> ts,
+    int const nytrue, int const nztrue, int const nmaxevent, std::vector<realtype> ts,
     std::vector<realtype> const& observedData,
     std::vector<realtype> const& observedDataStdDev,
     std::vector<realtype> const& observedEvents,
@@ -69,7 +69,7 @@ ExpData::ExpData(Model const& model)
 }
 
 ExpData::ExpData(
-    ReturnData const& rdata, realtype sigma_y, realtype sigma_z, int seed
+    ReturnData const& rdata, realtype const sigma_y, realtype const sigma_z, int const seed
 )
     : ExpData(
           rdata, std::vector<realtype>(rdata.nytrue * rdata.nt, sigma_y),
@@ -78,18 +78,18 @@ ExpData::ExpData(
 
 ExpData::ExpData(
     ReturnData const& rdata, std::vector<realtype> sigma_y,
-    std::vector<realtype> sigma_z, int seed
+    std::vector<realtype> sigma_z, int const seed
 )
     : ExpData(rdata.nytrue, rdata.nztrue, rdata.nmaxevent, rdata.ts) {
-    if (sigma_y.size() != (unsigned)nytrue_
-        && sigma_y.size() != (unsigned)nytrue_ * nt())
+    if (sigma_y.size() != static_cast<unsigned>(nytrue_)
+        && sigma_y.size() != static_cast<unsigned>(nytrue_) * nt())
         throw AmiException(
             "Dimension of sigma_y must be %d or %d, was %d", nytrue_,
             nytrue_ * nt(), sigma_y.size()
         );
 
-    if (sigma_z.size() != (unsigned)nztrue_
-        && sigma_z.size() != (unsigned)nztrue_ * nmaxevent_)
+    if (sigma_z.size() != static_cast<unsigned>(nztrue_)
+        && sigma_z.size() != static_cast<unsigned>(nztrue_) * nmaxevent_)
         throw AmiException(
             "Dimension of sigma_z must be %d or %d, was %d", nztrue_,
             nztrue_ * nmaxevent_, sigma_z.size()
@@ -104,7 +104,7 @@ ExpData::ExpData(
 
     for (int iy = 0; iy < nytrue_; ++iy) {
         for (int it = 0; it < nt(); ++it) {
-            sigma = sigma_y.size() == (unsigned)nytrue_
+            sigma = sigma_y.size() == static_cast<unsigned>(nytrue_)
                         ? sigma_y.at(iy)
                         : sigma_y.at(iy + nytrue_ * it);
             std::normal_distribution<> e{0, sigma};
@@ -116,7 +116,7 @@ ExpData::ExpData(
 
     for (int iz = 0; iz < nztrue_; ++iz) {
         for (int ie = 0; ie < nmaxevent_; ++ie) {
-            sigma = sigma_z.size() == (unsigned)nztrue_
+            sigma = sigma_z.size() == static_cast<unsigned>(nztrue_)
                         ? sigma_z.at(iz)
                         : sigma_z.at(iz + nztrue_ * ie);
             std::normal_distribution<> e{0, sigma};
@@ -143,12 +143,12 @@ std::vector<realtype> const& ExpData::getTimepoints() const { return ts_; }
 
 int ExpData::nt() const { return gsl::narrow<int>(ts_.size()); }
 
-realtype ExpData::getTimepoint(int it) const { return ts_.at(it); }
+realtype ExpData::getTimepoint(int const it) const { return ts_.at(it); }
 
 void ExpData::setObservedData(std::vector<realtype> const& observedData) {
     checkDataDimension(observedData, "observedData");
 
-    if (observedData.size() == (unsigned)nt() * nytrue_)
+    if (observedData.size() == static_cast<unsigned>(nt()) * nytrue_)
         observed_data_ = observedData;
     else if (observedData.empty())
         observed_data_.clear();
@@ -157,7 +157,7 @@ void ExpData::setObservedData(std::vector<realtype> const& observedData) {
 void ExpData::setObservedData(
     std::vector<realtype> const& observedData, int iy
 ) {
-    if (observedData.size() != (unsigned)nt())
+    if (observedData.size() != static_cast<unsigned>(nt()))
         throw AmiException(
             "Input observedData did not match dimensions nt (%i), was %i", nt(),
             observedData.size()
@@ -167,7 +167,7 @@ void ExpData::setObservedData(
         observed_data_.at(iy + it * nytrue_) = observedData.at(it);
 }
 
-bool ExpData::isSetObservedData(int it, int iy) const {
+bool ExpData::isSetObservedData(int const it, int const iy) const {
     return !observed_data_.empty()
            && !isNaN(observed_data_.at(it * nytrue_ + iy));
 }
@@ -176,7 +176,7 @@ std::vector<realtype> const& ExpData::getObservedData() const {
     return observed_data_;
 }
 
-realtype const* ExpData::getObservedDataPtr(int it) const {
+realtype const* ExpData::getObservedDataPtr(int const it) const {
     if (!observed_data_.empty())
         return &observed_data_.at(it * nytrue_);
 
@@ -189,7 +189,7 @@ void ExpData::setObservedDataStdDev(
     checkDataDimension(observedDataStdDev, "observedDataStdDev");
     checkSigmaPositivity(observedDataStdDev, "observedDataStdDev");
 
-    if (observedDataStdDev.size() == (unsigned)nt() * nytrue_)
+    if (observedDataStdDev.size() == static_cast<unsigned>(nt()) * nytrue_)
         observed_data_std_dev_ = observedDataStdDev;
     else if (observedDataStdDev.empty())
         observed_data_std_dev_.clear();
@@ -201,9 +201,9 @@ void ExpData::setObservedDataStdDev(realtype const stdDev) {
 }
 
 void ExpData::setObservedDataStdDev(
-    std::vector<realtype> const& observedDataStdDev, int iy
+    std::vector<realtype> const& observedDataStdDev, int const iy
 ) {
-    if (observedDataStdDev.size() != (unsigned)nt())
+    if (observedDataStdDev.size() != static_cast<unsigned>(nt()))
         throw AmiException(
             "Input observedDataStdDev did not match dimensions nt (%i), was %i",
             nt(), observedDataStdDev.size()
@@ -215,13 +215,13 @@ void ExpData::setObservedDataStdDev(
             = observedDataStdDev.at(it);
 }
 
-void ExpData::setObservedDataStdDev(realtype const stdDev, int iy) {
+void ExpData::setObservedDataStdDev(realtype const stdDev, int const iy) {
     checkSigmaPositivity(stdDev, "stdDev");
     for (int it = 0; it < nt(); ++it)
         observed_data_std_dev_.at(iy + it * nytrue_) = stdDev;
 }
 
-bool ExpData::isSetObservedDataStdDev(int it, int iy) const {
+bool ExpData::isSetObservedDataStdDev(int const it, int const iy) const {
     return !observed_data_std_dev_.empty()
            && !isNaN(observed_data_std_dev_.at(it * nytrue_ + iy));
 }
@@ -230,7 +230,7 @@ std::vector<realtype> const& ExpData::getObservedDataStdDev() const {
     return observed_data_std_dev_;
 }
 
-realtype const* ExpData::getObservedDataStdDevPtr(int it) const {
+realtype const* ExpData::getObservedDataStdDevPtr(int const it) const {
     if (!observed_data_std_dev_.empty())
         return &observed_data_std_dev_.at(it * nytrue_);
 
@@ -240,16 +240,16 @@ realtype const* ExpData::getObservedDataStdDevPtr(int it) const {
 void ExpData::setObservedEvents(std::vector<realtype> const& observedEvents) {
     checkEventsDimension(observedEvents, "observedEvents");
 
-    if (observedEvents.size() == (unsigned)nmaxevent_ * nztrue_)
+    if (observedEvents.size() == static_cast<unsigned>(nmaxevent_) * nztrue_)
         observed_events_ = observedEvents;
     else if (observedEvents.empty())
         observed_events_.clear();
 }
 
 void ExpData::setObservedEvents(
-    std::vector<realtype> const& observedEvents, int iz
+    std::vector<realtype> const& observedEvents, int const iz
 ) {
-    if (observedEvents.size() != (unsigned)nmaxevent_) {
+    if (observedEvents.size() != static_cast<unsigned>(nmaxevent_)) {
         throw AmiException(
             "Input observedEvents did not match dimensions nmaxevent (%i), was "
             "%i",
@@ -261,7 +261,7 @@ void ExpData::setObservedEvents(
         observed_events_.at(iz + ie * nztrue_) = observedEvents.at(ie);
 }
 
-bool ExpData::isSetObservedEvents(int ie, int iz) const {
+bool ExpData::isSetObservedEvents(int const ie, int const  iz) const {
     return !observed_events_.empty()
            && !isNaN(observed_events_.at(ie * nztrue_ + iz));
 }
@@ -270,7 +270,7 @@ std::vector<realtype> const& ExpData::getObservedEvents() const {
     return observed_events_;
 }
 
-realtype const* ExpData::getObservedEventsPtr(int ie) const {
+realtype const* ExpData::getObservedEventsPtr(int const  ie) const {
     if (!observed_events_.empty())
         return &observed_events_.at(ie * nztrue_);
 
@@ -295,7 +295,7 @@ void ExpData::setObservedEventsStdDev(realtype const stdDev) {
 }
 
 void ExpData::setObservedEventsStdDev(
-    std::vector<realtype> const& observedEventsStdDev, int iz
+    std::vector<realtype> const& observedEventsStdDev, int const iz
 ) {
     if (observedEventsStdDev.size() != (unsigned)nmaxevent_)
         throw AmiException(
@@ -310,14 +310,14 @@ void ExpData::setObservedEventsStdDev(
             = observedEventsStdDev.at(ie);
 }
 
-void ExpData::setObservedEventsStdDev(realtype const stdDev, int iz) {
+void ExpData::setObservedEventsStdDev(realtype const stdDev, int const iz) {
     checkSigmaPositivity(stdDev, "stdDev");
 
     for (int ie = 0; ie < nmaxevent_; ++ie)
         observed_events_std_dev_.at(iz + ie * nztrue_) = stdDev;
 }
 
-bool ExpData::isSetObservedEventsStdDev(int ie, int iz) const {
+bool ExpData::isSetObservedEventsStdDev(int const ie, int const iz) const {
     if (!observed_events_std_dev_.empty()) // avoid out of bound memory access
         return !isNaN(observed_events_std_dev_.at(ie * nztrue_ + iz));
 
@@ -328,7 +328,7 @@ std::vector<realtype> const& ExpData::getObservedEventsStdDev() const {
     return observed_events_std_dev_;
 }
 
-realtype const* ExpData::getObservedEventsStdDevPtr(int ie) const {
+realtype const* ExpData::getObservedEventsStdDevPtr(int const ie) const {
     if (!observed_events_std_dev_.empty())
         return &observed_events_std_dev_.at(ie * nztrue_);
 
@@ -360,7 +360,7 @@ void ExpData::applyEventDimension() {
 void ExpData::checkDataDimension(
     std::vector<realtype> const& input, char const* fieldname
 ) const {
-    if (input.size() != (unsigned)nt() * nytrue_ && !input.empty())
+    if (input.size() != static_cast<unsigned>(nt()) * nytrue_ && !input.empty())
         throw AmiException(
             "Input %s did not match dimensions nt (%i) x nytrue (%i), was %i",
             fieldname, nt(), nytrue_, input.size()
@@ -370,7 +370,7 @@ void ExpData::checkDataDimension(
 void ExpData::checkEventsDimension(
     std::vector<realtype> const& input, char const* fieldname
 ) const {
-    if (input.size() != (unsigned)nmaxevent_ * nztrue_ && !input.empty())
+    if (input.size() != static_cast<unsigned>(nmaxevent_) * nztrue_ && !input.empty())
         throw AmiException(
             "Input %s did not match dimensions nt (%i) x nytrue (%i), was %i",
             fieldname, nmaxevent_, nztrue_, input.size()
diff --git a/src/exception.cpp b/src/exception.cpp
index 984be52821..85c2f0ebd1 100644
--- a/src/exception.cpp
+++ b/src/exception.cpp
@@ -49,17 +49,17 @@ IDAException::IDAException(
           extra ? extra : ""
       ) {}
 
-IntegrationFailure::IntegrationFailure(int code, realtype t)
+IntegrationFailure::IntegrationFailure(int const code, realtype const t)
     : AmiException("AMICI failed to integrate the forward problem")
     , error_code(code)
     , time(t) {}
 
-IntegrationFailureB::IntegrationFailureB(int code, realtype t)
+IntegrationFailureB::IntegrationFailureB(int const code, realtype const t)
     : AmiException("AMICI failed to integrate the backward problem")
     , error_code(code)
     , time(t) {}
 
-NewtonFailure::NewtonFailure(int code, char const* function)
+NewtonFailure::NewtonFailure(int const code, char const* function)
     : AmiException(
           "NewtonSolver routine %s failed with error code %i", function, code
       ) {
diff --git a/src/hdf5.cpp b/src/hdf5.cpp
index 4ef03de2ec..53c46cc8f2 100644
--- a/src/hdf5.cpp
+++ b/src/hdf5.cpp
@@ -84,7 +84,7 @@ void createGroup(
     H5::H5File const& file, std::string const& groupPath, bool recursively
 ) {
 #if H5_VERSION_GE(1, 10, 6)
-    H5::LinkCreatPropList lcpl;
+    H5::LinkCreatPropList const lcpl;
     lcpl.setCreateIntermediateGroup(recursively);
     file.createGroup(groupPath.c_str(), lcpl);
 #else
@@ -113,7 +113,7 @@ std::unique_ptr<ExpData> readSimulationExpData(
     std::string const& hdf5Filename, std::string const& hdf5Root,
     Model const& model
 ) {
-    H5::H5File file(hdf5Filename.c_str(), H5F_ACC_RDONLY);
+    H5::H5File const file(hdf5Filename.c_str(), H5F_ACC_RDONLY);
 
     hsize_t m, n;
 
@@ -125,7 +125,7 @@ std::unique_ptr<ExpData> readSimulationExpData(
 
     if (model.ny * model.nt() > 0) {
         if (locationExists(file, hdf5Root + "/Y")) {
-            auto my = getDoubleDataset2D(file, hdf5Root + "/Y", m, n);
+            auto const my = getDoubleDataset2D(file, hdf5Root + "/Y", m, n);
             checkMeasurementDimensionsCompatible(m, n, model);
             edata->setObservedData(my);
         } else {
@@ -135,7 +135,7 @@ std::unique_ptr<ExpData> readSimulationExpData(
         }
 
         if (locationExists(file, hdf5Root + "/Sigma_Y")) {
-            auto sigmay = getDoubleDataset2D(file, hdf5Root + "/Sigma_Y", m, n);
+            auto const sigmay = getDoubleDataset2D(file, hdf5Root + "/Sigma_Y", m, n);
             checkMeasurementDimensionsCompatible(m, n, model);
             edata->setObservedDataStdDev(sigmay);
         } else {
@@ -148,7 +148,7 @@ std::unique_ptr<ExpData> readSimulationExpData(
 
     if (model.nz * model.nMaxEvent() > 0) {
         if (locationExists(file, hdf5Root + "/Z")) {
-            auto mz = getDoubleDataset2D(file, hdf5Root + "/Z", m, n);
+            auto const mz = getDoubleDataset2D(file, hdf5Root + "/Z", m, n);
             checkEventDimensionsCompatible(m, n, model);
             edata->setObservedEvents(mz);
         } else {
@@ -719,7 +719,7 @@ struct LogItemCStr {
 };
 
 void writeLogItemsToHDF5(
-    H5::H5File const& file, std::vector<amici::LogItem> const& logItems,
+    H5::H5File const& file, std::vector<LogItem> const& logItems,
     std::string const& hdf5Location
 ) {
     if (logItems.empty())
@@ -727,7 +727,7 @@ void writeLogItemsToHDF5(
 
     try {
         hsize_t dims[1] = {logItems.size()};
-        H5::DataSpace dataspace(1, dims);
+        const H5::DataSpace dataspace(1, dims);
 
         // works on Ubuntu, but segfaults on macos:
         /*
@@ -871,11 +871,11 @@ int getIntScalarAttribute(
 
 void createAndWriteInt1DDataset(
     H5::H5File const& file, std::string const& datasetName,
-    gsl::span<int const> buffer
+    gsl::span<int const> const buffer
 ) {
     hsize_t size = buffer.size();
     H5::DataSpace dataspace(1, &size);
-    auto dataset = file.createDataSet(
+    auto const dataset = file.createDataSet(
         datasetName.c_str(), H5::PredType::NATIVE_INT, dataspace
     );
     dataset.write(buffer.data(), H5::PredType::NATIVE_INT);
@@ -885,9 +885,9 @@ void createAndWriteDouble1DDataset(
     const H5::H5File& file, std::string const& datasetName,
     gsl::span<double const> buffer
 ) {
-    hsize_t size = buffer.size();
+    hsize_t const size = buffer.size();
     H5::DataSpace dataspace(1, &size);
-    auto dataset = file.createDataSet(
+    auto const dataset = file.createDataSet(
         datasetName.c_str(), H5::PredType::NATIVE_DOUBLE, dataspace
     );
     dataset.write(buffer.data(), H5::PredType::NATIVE_DOUBLE);
@@ -895,7 +895,7 @@ void createAndWriteDouble1DDataset(
 
 void createAndWriteDouble2DDataset(
     const H5::H5File& file, std::string const& datasetName,
-    gsl::span<double const> buffer, hsize_t const m, hsize_t const n
+    gsl::span<double const> const buffer, hsize_t const m, hsize_t const n
 ) {
     Expects(buffer.size() == m * n);
     hsize_t const adims[]{m, n};
@@ -908,7 +908,7 @@ void createAndWriteDouble2DDataset(
 
 void createAndWriteInt2DDataset(
     H5::H5File const& file, std::string const& datasetName,
-    gsl::span<int const> buffer, hsize_t const m, hsize_t const n
+    gsl::span<int const> const buffer, hsize_t const m, hsize_t const n
 ) {
     Expects(buffer.size() == m * n);
     hsize_t const adims[]{m, n};
@@ -921,7 +921,7 @@ void createAndWriteInt2DDataset(
 
 void createAndWriteDouble3DDataset(
     H5::H5File const& file, std::string const& datasetName,
-    gsl::span<double const> buffer, hsize_t const m, hsize_t const n,
+    gsl::span<double const> const buffer, hsize_t const m, hsize_t const n,
     hsize_t const o
 ) {
     Expects(buffer.size() == m * n * o);
@@ -958,7 +958,7 @@ void writeSolverSettingsToHDF5(
     Solver const& solver, std::string const& hdf5Filename,
     std::string const& hdf5Location
 ) {
-    auto file = createOrOpenForWriting(hdf5Filename);
+    auto const file = createOrOpenForWriting(hdf5Filename);
 
     writeSolverSettingsToHDF5(solver, file, hdf5Location);
 }
diff --git a/src/interface_matlab.cpp b/src/interface_matlab.cpp
index 5a6c587a7b..8ac331c484 100644
--- a/src/interface_matlab.cpp
+++ b/src/interface_matlab.cpp
@@ -615,7 +615,7 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, mxArray const* prhs[]) {
             "AMICI:mex:setup",
             "Incorrect number of input arguments (must be at least 7)!"
         );
-    };
+    }
 
     auto model = amici::generic_model::getModel();
     auto solver = model->getSolver();
diff --git a/src/misc.cpp b/src/misc.cpp
index 19df59fdad..3a1524a23f 100644
--- a/src/misc.cpp
+++ b/src/misc.cpp
@@ -20,7 +20,7 @@
 
 namespace amici {
 
-void writeSlice(AmiVector const& s, gsl::span<realtype> b) {
+void writeSlice(AmiVector const& s, gsl::span<realtype> const b) {
     writeSlice(s.getVector(), b);
 }
 
@@ -40,8 +40,8 @@ double getUnscaledParameter(
 }
 
 void unscaleParameters(
-    gsl::span<realtype const> bufferScaled,
-    gsl::span<ParameterScaling const> pscale, gsl::span<realtype> bufferUnscaled
+    gsl::span<realtype const> const bufferScaled,
+    gsl::span<ParameterScaling const> const pscale, gsl::span<realtype> const bufferUnscaled
 ) {
     Expects(bufferScaled.size() == pscale.size());
     Expects(bufferScaled.size() == bufferUnscaled.size());
@@ -52,7 +52,7 @@ void unscaleParameters(
     }
 }
 
-double getScaledParameter(double unscaledParameter, ParameterScaling scaling) {
+double getScaledParameter(double const unscaledParameter, ParameterScaling const scaling) {
     switch (scaling) {
     case ParameterScaling::log10:
         return log10(unscaledParameter);
@@ -66,8 +66,8 @@ double getScaledParameter(double unscaledParameter, ParameterScaling scaling) {
 }
 
 void scaleParameters(
-    gsl::span<realtype const> bufferUnscaled,
-    gsl::span<ParameterScaling const> pscale, gsl::span<realtype> bufferScaled
+    gsl::span<realtype const> const bufferUnscaled,
+    gsl::span<ParameterScaling const> const pscale, gsl::span<realtype> const bufferScaled
 ) {
     Expects(bufferScaled.size() == pscale.size());
     Expects(bufferScaled.size() == bufferUnscaled.size());
@@ -127,7 +127,7 @@ std::string backtraceString(int const maxFrames, int const first_frame) {
     return trace_buf.str();
 }
 
-std::string regexErrorToString(std::regex_constants::error_type err_type) {
+std::string regexErrorToString(std::regex_constants::error_type const err_type) {
     switch (err_type) {
     case std::regex_constants::error_collate:
         return "error_collate";
@@ -177,7 +177,7 @@ std::string printfToString(char const* fmt, va_list ap) {
     return str;
 }
 
-std::pair<size_t, size_t> unravel_index(size_t flat_idx, size_t num_cols) {
+std::pair<size_t, size_t> unravel_index(size_t const flat_idx, size_t const num_cols) {
     return {flat_idx / num_cols, flat_idx % num_cols};
 }
 
diff --git a/src/model.template.cpp b/src/model.template.cpp
index b00f19b52f..eb158d1745 100644
--- a/src/model.template.cpp
+++ b/src/model.template.cpp
@@ -1,9 +1,7 @@
 #include <amici/defines.h>
 #include <array>
 
-namespace amici {
-
-namespace model_TPL_MODELNAME {
+namespace amici::model_TPL_MODELNAME {
 
 // clang-format off
 
@@ -57,6 +55,4 @@ std::array<int, TPL_NX_SOLVER> stateIdxsSolver = {
 
 // clang-format on
 
-} // namespace model_TPL_MODELNAME
-
-} // namespace amici
+} // namespace amici::model_TPL_MODELNAME
diff --git a/src/model_dae.cpp b/src/model_dae.cpp
index b49c52fd2c..bc32545aa5 100644
--- a/src/model_dae.cpp
+++ b/src/model_dae.cpp
@@ -11,7 +11,7 @@ void Model_DAE::fJ(
 }
 
 void Model_DAE::fJ(
-    realtype t, realtype cj, const_N_Vector x, const_N_Vector dx,
+    realtype const t, realtype const cj, const_N_Vector x, const_N_Vector dx,
     const_N_Vector /*xdot*/, SUNMatrix J
 ) {
     fJSparse(t, cj, x, dx, derived_state_.J_);
@@ -28,9 +28,9 @@ void Model_DAE::fJSparse(
 }
 
 void Model_DAE::fJSparse(
-    realtype t, realtype cj, const_N_Vector x, const_N_Vector dx, SUNMatrix J
+    realtype const t, realtype const cj, const_N_Vector x, const_N_Vector dx, SUNMatrix J
 ) {
-    auto x_pos = computeX_pos(x);
+    auto const x_pos = computeX_pos(x);
     fdwdx(t, N_VGetArrayPointerConst(x_pos), false);
     if (pythonGenerated) {
         auto JSparse = SUNMatrixWrapper(J);
@@ -84,8 +84,8 @@ void Model_DAE::fJv(
 }
 
 void Model_DAE::fJv(
-    realtype t, const_N_Vector x, const_N_Vector dx, const_N_Vector v,
-    N_Vector Jv, realtype cj
+    realtype const t, const_N_Vector x, const_N_Vector dx, const_N_Vector v,
+    N_Vector Jv, realtype const cj
 ) {
     N_VConst(0.0, Jv);
     fJSparse(t, cj, x, dx, derived_state_.J_);
@@ -95,16 +95,16 @@ void Model_DAE::fJv(
 
 void Model_DAE::froot(
     realtype const t, AmiVector const& x, AmiVector const& dx,
-    gsl::span<realtype> root
+    gsl::span<realtype> const root
 ) {
     froot(t, x.getNVector(), dx.getNVector(), root);
 }
 
 void Model_DAE::froot(
-    realtype t, const_N_Vector x, const_N_Vector dx, gsl::span<realtype> root
+    realtype const t, const_N_Vector x, const_N_Vector dx, gsl::span<realtype> root
 ) {
     std::ranges::fill(root, 0.0);
-    auto x_pos = computeX_pos(x);
+    auto const x_pos = computeX_pos(x);
     froot(
         root.data(), t, N_VGetArrayPointerConst(x_pos),
         state_.unscaledParameters.data(), state_.fixedParameters.data(),
@@ -119,9 +119,9 @@ void Model_DAE::fxdot(
 }
 
 void Model_DAE::fxdot(
-    realtype t, const_N_Vector x, const_N_Vector dx, N_Vector xdot
+    realtype const t, const_N_Vector x, const_N_Vector dx, N_Vector xdot
 ) {
-    auto x_pos = computeX_pos(x);
+    auto const x_pos = computeX_pos(x);
     fw(t, N_VGetArrayPointerConst(x), false);
     N_VConst(0.0, xdot);
     fxdot(
@@ -148,7 +148,7 @@ void Model_DAE::fdxdotdw(
 ) {
     derived_state_.dxdotdw_.zero();
     if (nw > 0 && derived_state_.dxdotdw_.capacity()) {
-        auto x_pos = computeX_pos(x);
+        auto const x_pos = computeX_pos(x);
 
         fdxdotdw_colptrs(derived_state_.dxdotdw_);
         fdxdotdw_rowvals(derived_state_.dxdotdw_);
@@ -164,7 +164,7 @@ void Model_DAE::fdxdotdw(
 void Model_DAE::fdxdotdp(
     realtype const t, const_N_Vector const x, const_N_Vector const dx
 ) {
-    auto x_pos = computeX_pos(x);
+    auto const x_pos = computeX_pos(x);
     fdwdp(t, N_VGetArrayPointerConst(x_pos));
 
     if (pythonGenerated) {
@@ -210,7 +210,7 @@ void Model_DAE::fdxdotdp(
     }
 }
 
-void Model_DAE::fM(realtype t, const_N_Vector x) {
+void Model_DAE::fM(realtype const t, const_N_Vector x) {
     if (pythonGenerated) {
         /*
          * non-algebraic states in python generated code always have factor
@@ -230,7 +230,7 @@ void Model_DAE::fM(realtype t, const_N_Vector x) {
         assert(ndiff == derived_state_.MSparse_.capacity());
     } else {
         derived_state_.M_.zero();
-        auto x_pos = computeX_pos(x);
+        auto const x_pos = computeX_pos(x);
         fM(derived_state_.M_.data(), t, N_VGetArrayPointerConst(x_pos),
            state_.unscaledParameters.data(), state_.fixedParameters.data());
     }
@@ -376,7 +376,7 @@ void Model_DAE::fJB(
 }
 
 void Model_DAE::fJB(
-    realtype t, realtype cj, const_N_Vector x, const_N_Vector dx,
+    realtype const t, realtype const cj, const_N_Vector x, const_N_Vector dx,
     const_N_Vector /*xB*/, const_N_Vector /*dxB*/, SUNMatrix JB
 ) {
     fJSparse(t, cj, x, dx, derived_state_.J_);
@@ -386,7 +386,7 @@ void Model_DAE::fJB(
 }
 
 void Model_DAE::fJSparseB(
-    realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
+    realtype const t, realtype const cj, AmiVector const& x, AmiVector const& dx,
     AmiVector const& xB, AmiVector const& dxB, AmiVector const& /*xBdot*/,
     SUNMatrix JB
 ) {
@@ -397,7 +397,7 @@ void Model_DAE::fJSparseB(
 }
 
 void Model_DAE::fJSparseB(
-    realtype t, realtype cj, const_N_Vector x, const_N_Vector dx,
+    realtype const t, realtype const cj, const_N_Vector x, const_N_Vector dx,
     const_N_Vector /*xB*/, const_N_Vector /*dxB*/, SUNMatrix JB
 ) {
     fJSparse(t, cj, x, dx, derived_state_.J_);
@@ -407,7 +407,7 @@ void Model_DAE::fJSparseB(
 }
 
 void Model_DAE::fJvB(
-    realtype t, const_N_Vector x, const_N_Vector dx, const_N_Vector xB,
+    realtype const t, const_N_Vector x, const_N_Vector dx, const_N_Vector xB,
     const_N_Vector dxB, const_N_Vector vB, N_Vector JvB, realtype cj
 ) {
     N_VConst(0.0, JvB);
@@ -417,7 +417,7 @@ void Model_DAE::fJvB(
 }
 
 void Model_DAE::fxBdot(
-    realtype t, const_N_Vector x, const_N_Vector dx, const_N_Vector xB,
+    realtype const t, const_N_Vector x, const_N_Vector dx, const_N_Vector xB,
     const_N_Vector dxB, N_Vector xBdot
 ) {
     N_VConst(0.0, xBdot);
@@ -428,7 +428,7 @@ void Model_DAE::fxBdot(
 }
 
 void Model_DAE::fqBdot(
-    realtype t, const_N_Vector x, const_N_Vector dx, const_N_Vector xB,
+    realtype const t, const_N_Vector x, const_N_Vector dx, const_N_Vector xB,
     const_N_Vector /*dxB*/, N_Vector qBdot
 ) {
     N_VConst(0.0, qBdot);
diff --git a/src/model_ode.cpp b/src/model_ode.cpp
index b4908dbe63..525ce8983d 100644
--- a/src/model_ode.cpp
+++ b/src/model_ode.cpp
@@ -12,7 +12,7 @@ void Model_ODE::fJ(
 }
 
 void Model_ODE::fJ(
-    realtype t, const_N_Vector x, const_N_Vector /*xdot*/, SUNMatrix J
+    realtype const t, const_N_Vector x, const_N_Vector /*xdot*/, SUNMatrix J
 ) {
     fJSparse(t, x, derived_state_.J_);
     derived_state_.J_.refresh();
@@ -27,7 +27,7 @@ void Model_ODE::fJSparse(
     fJSparse(t, x.getNVector(), J);
 }
 
-void Model_ODE::fJSparse(realtype t, const_N_Vector x, SUNMatrix J) {
+void Model_ODE::fJSparse(realtype const t, const_N_Vector x, SUNMatrix J) {
     auto const x_pos = computeX_pos(x);
     fdwdx(t, N_VGetArrayPointerConst(x_pos), false);
     if (pythonGenerated) {
@@ -74,7 +74,7 @@ void Model_ODE::
 }
 
 void Model_ODE::fJv(
-    const_N_Vector v, N_Vector Jv, realtype t, const_N_Vector x
+    const_N_Vector v, N_Vector Jv, realtype const t, const_N_Vector x
 ) {
     N_VConst(0.0, Jv);
     fJSparse(t, x, derived_state_.J_);
@@ -89,8 +89,8 @@ void Model_ODE::froot(
     froot(t, x.getNVector(), root);
 }
 
-void Model_ODE::froot(realtype t, const_N_Vector x, gsl::span<realtype> root) {
-    auto x_pos = computeX_pos(x);
+void Model_ODE::froot(realtype const t, const_N_Vector x, gsl::span<realtype> root) {
+    auto const x_pos = computeX_pos(x);
     std::ranges::fill(root, 0.0);
     froot(
         root.data(), t, N_VGetArrayPointerConst(x_pos),
@@ -106,8 +106,8 @@ void Model_ODE::fxdot(
     fxdot(t, x.getNVector(), xdot.getNVector());
 }
 
-void Model_ODE::fxdot(realtype t, const_N_Vector x, N_Vector xdot) {
-    auto x_pos = computeX_pos(x);
+void Model_ODE::fxdot(realtype const t, const_N_Vector x, N_Vector xdot) {
+    auto const x_pos = computeX_pos(x);
     fw(t, N_VGetArrayPointerConst(x_pos));
     N_VConst(0.0, xdot);
     fxdot(
@@ -130,7 +130,7 @@ void Model_ODE::fJDiag(
 void Model_ODE::fdxdotdw(realtype const t, const_N_Vector x) {
     derived_state_.dxdotdw_.zero();
     if (nw > 0 && derived_state_.dxdotdw_.capacity()) {
-        auto x_pos = computeX_pos(x);
+        auto const x_pos = computeX_pos(x);
 
         fdxdotdw_colptrs(derived_state_.dxdotdw_);
         fdxdotdw_rowvals(derived_state_.dxdotdw_);
@@ -143,7 +143,7 @@ void Model_ODE::fdxdotdw(realtype const t, const_N_Vector x) {
 }
 
 void Model_ODE::fdxdotdp(realtype const t, const_N_Vector x) {
-    auto x_pos = computeX_pos(x);
+    auto const x_pos = computeX_pos(x);
     fdwdp(t, N_VGetArrayPointerConst(x_pos), false);
 
     if (pythonGenerated) {
@@ -193,7 +193,7 @@ void Model_ODE::
 }
 
 std::unique_ptr<Solver> Model_ODE::getSolver() {
-    return std::unique_ptr<Solver>(new amici::CVodeSolver());
+    return std::unique_ptr<Solver>(new CVodeSolver());
 }
 
 void Model_ODE::fJSparse(
diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp
index 0b94f2a494..2fa0a9adee 100644
--- a/src/newton_solver.cpp
+++ b/src/newton_solver.cpp
@@ -12,7 +12,7 @@
 namespace amici {
 
 NewtonSolver::NewtonSolver(
-    Model const& model, LinearSolver const linsol_type, SUNContext sunctx
+    Model const& model, LinearSolver const linsol_type, SUNContext const sunctx
 )
     : xdot_(model.nx_solver, sunctx)
     , x_(model.nx_solver, sunctx)
diff --git a/src/simulation_parameters.cpp b/src/simulation_parameters.cpp
index 0262b5f6b4..c56b2fa3db 100644
--- a/src/simulation_parameters.cpp
+++ b/src/simulation_parameters.cpp
@@ -26,7 +26,7 @@ bool operator==(SimulationParameters const& a, SimulationParameters const& b) {
 
 void SimulationParameters::
     reinitializeAllFixedParameterDependentInitialStatesForPresimulation(
-        int nx_rdata
+        int const nx_rdata
     ) {
     reinitialization_state_idxs_presim.resize(nx_rdata);
     std::iota(
@@ -37,7 +37,7 @@ void SimulationParameters::
 
 void SimulationParameters::
     reinitializeAllFixedParameterDependentInitialStatesForSimulation(
-        int nx_rdata
+        int const nx_rdata
     ) {
     reinitialization_state_idxs_sim.resize(nx_rdata);
     std::iota(
@@ -47,7 +47,7 @@ void SimulationParameters::
 }
 
 void SimulationParameters::reinitializeAllFixedParameterDependentInitialStates(
-    int nx_rdata
+    int const nx_rdata
 ) {
     reinitializeAllFixedParameterDependentInitialStatesForPresimulation(
         nx_rdata
diff --git a/src/solver.cpp b/src/solver.cpp
index aebf3a1cba..147f4f5753 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -15,8 +15,8 @@ namespace amici {
 
 /* Error handler passed to SUNDIALS. */
 void wrapErrHandlerFn(
-    [[maybe_unused]] int line, char const* func, char const* file,
-    char const* msg, SUNErrCode err_code, void* err_user_data,
+    [[maybe_unused]] int const line, char const* func, char const* file,
+    char const* msg, SUNErrCode const err_code, void* err_user_data,
     [[maybe_unused]] SUNContext sunctx
 ) {
     constexpr int BUF_SIZE = 250;
@@ -33,7 +33,7 @@ void wrapErrHandlerFn(
     // we need a matlab-compatible message ID
     // i.e. colon separated and only  [A-Za-z0-9_]
     // see https://mathworks.com/help/matlab/ref/mexception.html
-    std::filesystem::path path(file);
+    std::filesystem::path const path(file);
     auto file_stem = path.stem().string();
 
     // Error code to string. Remove 'AMICI_' prefix.
@@ -139,7 +139,7 @@ void Solver::apply_max_num_steps_B() const {
 
 int Solver::run(realtype const tout) const {
     setStopTime(tout);
-    CpuTimer cpu_timer;
+    CpuTimer const cpu_timer;
     int status = AMICI_SUCCESS;
 
     apply_max_num_steps();
@@ -672,8 +672,8 @@ void Solver::applyQuadTolerancesASA(int const which) const {
     if (sensi_ < SensitivityOrder::first)
         return;
 
-    realtype quad_rtol = isNaN(quad_rtol_) ? rtol_ : quad_rtol_;
-    realtype quad_atol = isNaN(quad_atol_) ? atol_ : quad_atol_;
+    realtype const quad_rtol = isNaN(quad_rtol_) ? rtol_ : quad_rtol_;
+    realtype const quad_atol = isNaN(quad_atol_) ? atol_ : quad_atol_;
 
     /* Enable Quadrature Error Control */
     setQuadErrConB(which, !std::isinf(quad_atol) && !std::isinf(quad_rtol));
@@ -740,7 +740,7 @@ void Solver::setSensitivityMethodPreequilibration(
 }
 
 void Solver::checkSensitivityMethod(
-    SensitivityMethod const sensi_meth, bool preequilibration
+    SensitivityMethod const sensi_meth, bool const preequilibration
 ) const {
     if (rdata_mode_ == RDataReporting::residuals
         && sensi_meth == SensitivityMethod::adjoint)
@@ -752,7 +752,7 @@ void Solver::checkSensitivityMethod(
         resetMutableMemory(nx(), nplist(), nquad());
 }
 
-void Solver::setMaxNonlinIters(int max_nonlin_iters) {
+void Solver::setMaxNonlinIters(int const max_nonlin_iters) {
     if (max_nonlin_iters < 0)
         throw AmiException("max_nonlin_iters must be a non-negative number");
 
@@ -761,7 +761,7 @@ void Solver::setMaxNonlinIters(int max_nonlin_iters) {
 
 int Solver::getMaxNonlinIters() const { return max_nonlin_iters_; }
 
-void Solver::setMaxConvFails(int max_conv_fails) {
+void Solver::setMaxConvFails(int const max_conv_fails) {
     if (max_conv_fails < 0)
         throw AmiException("max_conv_fails must be a non-negative number");
 
@@ -813,7 +813,7 @@ double Solver::getNewtonDampingFactorLowerBound() const {
     return newton_damping_factor_lower_bound_;
 }
 
-void Solver::setNewtonDampingFactorLowerBound(double dampingFactorLowerBound) {
+void Solver::setNewtonDampingFactorLowerBound(double const dampingFactorLowerBound) {
     newton_damping_factor_lower_bound_ = dampingFactorLowerBound;
 }
 
@@ -848,7 +848,7 @@ double Solver::getAbsoluteTolerance() const {
     return static_cast<double>(atol_);
 }
 
-void Solver::setAbsoluteTolerance(double atol) {
+void Solver::setAbsoluteTolerance(double const atol) {
     if (atol < 0)
         throw AmiException("atol must be a non-negative number");
 
@@ -960,11 +960,11 @@ double Solver::getSteadyStateToleranceFactor() const {
     return static_cast<double>(ss_tol_factor_);
 }
 
-void Solver::setSteadyStateToleranceFactor(double const ss_tol_factor) {
-    if (ss_tol_factor < 0)
+void Solver::setSteadyStateToleranceFactor(double const factor) {
+    if (factor < 0)
         throw AmiException("ss_tol_factor must be a non-negative number");
 
-    ss_tol_factor_ = static_cast<realtype>(ss_tol_factor);
+    ss_tol_factor_ = static_cast<realtype>(factor);
 }
 
 double Solver::getRelativeToleranceSteadyState() const {
@@ -1036,13 +1036,13 @@ long int Solver::getMaxSteps() const { return maxsteps_; }
 
 double Solver::getMaxTime() const { return maxtime_.count(); }
 
-void Solver::setMaxTime(double maxtime) {
+void Solver::setMaxTime(double const maxtime) {
     maxtime_ = std::chrono::duration<double>(maxtime);
 }
 
 void Solver::startTimer() const { simulation_timer_.reset(); }
 
-bool Solver::timeExceeded(int interval) const {
+bool Solver::timeExceeded(int const interval) const {
     thread_local int eval_counter = 0;
 
     // 0 means infinite time
@@ -1139,7 +1139,7 @@ void Solver::setStabilityLimitFlag(bool const stldet) {
 
 LinearSolver Solver::getLinearSolver() const { return linsol_; }
 
-void Solver::setLinearSolver(LinearSolver linsol) {
+void Solver::setLinearSolver(LinearSolver const linsol) {
     if (solver_memory_)
         resetMutableMemory(nx(), nplist(), nquad());
     linsol_ = linsol;
@@ -1234,7 +1234,7 @@ int Solver::nx() const { return x_.getLength(); }
 
 int Solver::nquad() const { return xQB_.getLength(); }
 
-bool Solver::getInitDone() const { return initialized_; };
+bool Solver::getInitDone() const { return initialized_; }
 
 bool Solver::getSensInitDone() const { return sens_initialized_; }
 
@@ -1252,7 +1252,7 @@ bool Solver::getQuadInitDoneB(int const which) const {
 
 bool Solver::getQuadInitDone() const { return quad_initialized_; }
 
-void Solver::setInitDone() const { initialized_ = true; };
+void Solver::setInitDone() const { initialized_ = true; }
 
 void Solver::setSensInitDone() const { sens_initialized_ = true; }
 
diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp
index 63d8af121c..1c14604bb9 100644
--- a/src/solver_idas.cpp
+++ b/src/solver_idas.cpp
@@ -223,7 +223,7 @@ void IDASolver::qbinit(int const which, AmiVector const& xQB0) const {
         throw IDAException(status, "IDAQuadInitB");
 }
 
-void IDASolver::rootInit(int ne) const {
+void IDASolver::rootInit(int const ne) const {
     int status = IDARootInit(solver_memory_.get(), ne, froot);
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDARootInit");
@@ -370,7 +370,7 @@ void IDASolver::setUserData() const {
         throw IDAException(status, "IDASetUserData");
 }
 
-void IDASolver::setUserDataB(int which) const {
+void IDASolver::setUserDataB(int const which) const {
     int status = IDASetUserDataB(solver_memory_.get(), which, &user_data);
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDASetUserDataB");
@@ -610,14 +610,14 @@ void IDASolver::getB(int const which) const {
         throw IDAException(status, "IDAGetB");
 }
 
-void IDASolver::getDkyB(realtype const t, int k, int const which) const {
+void IDASolver::getDkyB(realtype const t, int const k, int const which) const {
     int status = IDAGetDky(
         IDAGetAdjIDABmem(solver_memory_.get(), which), t, k, dky_.getNVector()
     );
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDAGetB");
 }
-void IDASolver::getQuadB(int which) const {
+void IDASolver::getQuadB(int const which) const {
     realtype tDummy = 0;
     int status
         = IDAGetQuadB(solver_memory_.get(), which, &tDummy, xQB_.getNVector());
@@ -631,7 +631,7 @@ void IDASolver::getQuad(realtype& t) const {
         throw IDAException(status, "IDAGetQuad");
 }
 
-void IDASolver::getQuadDkyB(realtype const t, int k, int const which) const {
+void IDASolver::getQuadDkyB(realtype const t, int const k, int const which) const {
     int status = IDAGetQuadDky(
         IDAGetAdjIDABmem(solver_memory_.get(), which), t, k, xQB_.getNVector()
     );
@@ -710,9 +710,9 @@ void IDASolver::quadSStolerancesB(
 }
 
 void IDASolver::quadSStolerances(
-    realtype const reltolQB, realtype const abstolQB
+    realtype const reltolQ, realtype const abstolQ
 ) const {
-    int status = IDAQuadSStolerances(solver_memory_.get(), reltolQB, abstolQB);
+    int status = IDAQuadSStolerances(solver_memory_.get(), reltolQ, abstolQ);
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDAQuadSStolerances");
 }
@@ -854,9 +854,8 @@ Model const* IDASolver::getModel() const {
         throw AmiException(
             "Solver has not been allocated, information is not available"
         );
-    auto ida_mem = static_cast<IDAMem>(solver_memory_.get());
-    auto user_data = static_cast<user_data_type*>(ida_mem->ida_user_data);
-    if (user_data)
+    auto const ida_mem = static_cast<IDAMem>(solver_memory_.get());
+    if (auto user_data = static_cast<user_data_type*>(ida_mem->ida_user_data))
         return user_data->first;
     return nullptr;
 }
diff --git a/src/splinefunctions.cpp b/src/splinefunctions.cpp
index ead4329d94..ed7b43d53b 100644
--- a/src/splinefunctions.cpp
+++ b/src/splinefunctions.cpp
@@ -490,8 +490,8 @@ void HermiteSpline::compute_slope_sensitivities_by_fd(
     int nplist, int spline_offset, int ip, gsl::span<realtype> dvaluesdp,
     gsl::span<realtype> dslopesdp
 ) {
-    int last = n_nodes() - 1;
-    int node_offset = spline_offset + ip;
+    int const last = n_nodes() - 1;
+    int const node_offset = spline_offset + ip;
 
     // Left boundary (first node)
     switch (first_node_bc_) {
diff --git a/src/sundials_linsol_wrapper.cpp b/src/sundials_linsol_wrapper.cpp
index 5befe85499..96dd480515 100644
--- a/src/sundials_linsol_wrapper.cpp
+++ b/src/sundials_linsol_wrapper.cpp
@@ -15,7 +15,7 @@
 
 namespace amici {
 
-SUNLinSolWrapper::SUNLinSolWrapper(SUNLinearSolver linsol)
+SUNLinSolWrapper::SUNLinSolWrapper(SUNLinearSolver const linsol)
     : solver_(linsol) {}
 
 SUNLinSolWrapper::SUNLinSolWrapper(
diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp
index 808dd9ec2a..74689bd5e1 100644
--- a/src/sundials_matrix_wrapper.cpp
+++ b/src/sundials_matrix_wrapper.cpp
@@ -148,22 +148,22 @@ SUNMatrixWrapper& SUNMatrixWrapper::operator=(SUNMatrixWrapper&& other) {
     return *this;
 }
 
-void SUNMatrixWrapper::reallocate(sunindextype NNZ) {
+void SUNMatrixWrapper::reallocate(sunindextype nnz) {
     if (sparsetype() != CSC_MAT && sparsetype() != CSR_MAT)
         throw std::invalid_argument(
             "Invalid sparsetype. Must be CSC_MAT or "
             "CSR_MAT."
         );
 
-    if (int ret = SUNSparseMatrix_Reallocate(matrix_, NNZ) != SUN_SUCCESS)
+    if (int ret = SUNSparseMatrix_Reallocate(matrix_, nnz) != SUN_SUCCESS)
         throw std::runtime_error(
             "SUNSparseMatrix_Reallocate failed with error code "
             + std::to_string(ret) + "."
         );
 
     update_ptrs();
-    capacity_ = NNZ;
-    assert((NNZ && columns() && rows()) || !matrix_);
+    capacity_ = nnz;
+    assert((nnz && columns() && rows()) || !matrix_);
 }
 
 void SUNMatrixWrapper::realloc() {
@@ -220,7 +220,7 @@ int SUNMatrixWrapper::sparsetype() const {
 
 void SUNMatrixWrapper::scale(realtype a) {
     if (matrix_) {
-        auto cap = capacity();
+        auto const cap = capacity();
         for (sunindextype idx = 0; idx < cap; ++idx)
             data_[idx] *= a;
     }
@@ -441,9 +441,6 @@ void SUNMatrixWrapper::sparse_add(
 
     sunindextype nnz = 0; // this keeps track of the nonzero index in C
 
-    sunindextype ccol;
-    sunindextype cidx;
-
     // first call, make sure that matrix is initialized with no capacity
     if (!capacity())
         reallocate(A.num_nonzeros() + B.num_nonzeros());
@@ -451,7 +448,7 @@ void SUNMatrixWrapper::sparse_add(
     auto w = std::vector<sunindextype>(rows());
     auto x = std::vector<realtype>(rows());
 
-    for (ccol = 0; ccol < columns(); ccol++) {
+    for (sunindextype ccol = 0; ccol < columns(); ccol++) {
         set_indexptr(ccol, nnz); /* column j of C starts here */
         nnz = A.scatter(
             ccol, alpha, w.data(), gsl::make_span(x), ccol + 1, this, nnz
@@ -460,7 +457,7 @@ void SUNMatrixWrapper::sparse_add(
             ccol, beta, w.data(), gsl::make_span(x), ccol + 1, this, nnz
         );
         // no reallocation should happen here
-        for (cidx = get_indexptr(ccol); cidx < nnz; cidx++) {
+        for (sunindextype cidx = get_indexptr(ccol); cidx < nnz; cidx++) {
             auto x_idx = get_indexval(cidx);
             assert(x_idx >= 0 && gsl::narrow<std::size_t>(x_idx) < x.size());
             set_data(cidx, x[x_idx]); // copy data to C
diff --git a/src/vector.cpp b/src/vector.cpp
index d93769fced..217a020f9c 100644
--- a/src/vector.cpp
+++ b/src/vector.cpp
@@ -29,21 +29,21 @@ void AmiVector::minus() {
     std::ranges::transform(vec_, vec_.begin(), std::negate<realtype>());
 }
 
-void AmiVector::set(realtype val) { std::ranges::fill(vec_, val); }
+void AmiVector::set(realtype const val) { std::ranges::fill(vec_, val); }
 
-realtype& AmiVector::operator[](int pos) {
+realtype& AmiVector::operator[](int const pos) {
     return vec_.at(gsl::narrow<decltype(vec_)::size_type>(pos));
 }
 
-realtype const& AmiVector::operator[](int pos) const {
+realtype const& AmiVector::operator[](int const pos) const {
     return vec_.at(gsl::narrow<decltype(vec_)::size_type>(pos));
 }
 
-realtype& AmiVector::at(int pos) {
+realtype& AmiVector::at(int const pos) {
     return vec_.at(gsl::narrow<decltype(vec_)::size_type>(pos));
 }
 
-realtype const& AmiVector::at(int pos) const {
+realtype const& AmiVector::at(int const pos) const {
     return vec_.at(gsl::narrow<decltype(vec_)::size_type>(pos));
 }
 
@@ -57,7 +57,7 @@ void AmiVector::copy(AmiVector const& other) {
     std::ranges::copy(other.vec_, vec_.begin());
 }
 
-void AmiVector::synchroniseNVector(SUNContext sunctx) {
+void AmiVector::synchroniseNVector(SUNContext const sunctx) {
     if (nvec_)
         N_VDestroy_Serial(nvec_);
     if (sunctx) {
@@ -73,7 +73,7 @@ AmiVector::~AmiVector() {
 }
 
 AmiVectorArray::AmiVectorArray(
-    long int length_inner, long int length_outer, SUNContext sunctx
+    long int const length_inner, long int const length_outer, SUNContext const sunctx
 )
     : vec_array_(length_outer, AmiVector(length_inner, sunctx)) {
     nvec_array_.resize(length_outer);
@@ -101,29 +101,29 @@ AmiVectorArray::AmiVectorArray(AmiVectorArray const& vaold)
 
 realtype* AmiVectorArray::data(int pos) { return vec_array_.at(pos).data(); }
 
-realtype const* AmiVectorArray::data(int pos) const {
+realtype const* AmiVectorArray::data(int const pos) const {
     return vec_array_.at(pos).data();
 }
 
-realtype& AmiVectorArray::at(int ipos, int jpos) {
+realtype& AmiVectorArray::at(int const ipos, int const jpos) {
     return vec_array_.at(jpos).at(ipos);
 }
 
-realtype const& AmiVectorArray::at(int ipos, int jpos) const {
+realtype const& AmiVectorArray::at(int const ipos, int const jpos) const {
     return vec_array_.at(jpos).at(ipos);
 }
 
 N_Vector* AmiVectorArray::getNVectorArray() { return nvec_array_.data(); }
 
-N_Vector AmiVectorArray::getNVector(int pos) { return nvec_array_.at(pos); }
+N_Vector AmiVectorArray::getNVector(int const pos) { return nvec_array_.at(pos); }
 
-const_N_Vector AmiVectorArray::getNVector(int pos) const {
+const_N_Vector AmiVectorArray::getNVector(int const pos) const {
     return nvec_array_.at(pos);
 }
 
-AmiVector& AmiVectorArray::operator[](int pos) { return vec_array_.at(pos); }
+AmiVector& AmiVectorArray::operator[](int const pos) { return vec_array_.at(pos); }
 
-AmiVector const& AmiVectorArray::operator[](int pos) const {
+AmiVector const& AmiVectorArray::operator[](int const pos) const {
     return vec_array_.at(pos);
 }
 
diff --git a/src/wrapfunctions.template.h b/src/wrapfunctions.template.h
index 33c8a9d819..e7ba14ce79 100644
--- a/src/wrapfunctions.template.h
+++ b/src/wrapfunctions.template.h
@@ -5,9 +5,7 @@
 
 #include "amici/model.h"
 
-namespace amici {
-namespace generic_model {
-
+namespace amici::generic_model {
 
 /**
  * @brief Wrapper function to instantiate the linked Amici model without knowing
@@ -16,10 +14,7 @@ namespace generic_model {
  */
 std::unique_ptr<amici::Model> getModel();
 
-
-} // namespace generic_model
-
-} // namespace amici
+} // namespace amici::generic_model
 
 
 #endif /* _amici_wrapfunctions_h */
diff --git a/swig/stdvec2numpy.h b/swig/stdvec2numpy.h
index 2a67a8f478..b051e9ac9a 100644
--- a/swig/stdvec2numpy.h
+++ b/swig/stdvec2numpy.h
@@ -10,7 +10,7 @@ PyObject* stdVec2ndarray(std::vector<double>& vec, int dim1) {
     if (vec.size() != (unsigned) dim1) throw std::runtime_error("Size mismatch in stdVec2ndarray");
     npy_intp dims[1] = { dim1 };
     PyObject * array = PyArray_SimpleNewFromData(1, dims,  NPY_DOUBLE, vec.data());
-    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");;
+    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");
     return array;
 }
 
@@ -25,7 +25,7 @@ PyObject* stdVec2ndarray(std::vector<double>& vec, int dim1, int dim2) {
     if (vec.size() != (unsigned) dim1 * dim2) throw std::runtime_error("Size mismatch in stdVec2ndarray");
     npy_intp dims[2] = { dim1, dim2 };
     PyObject * array = PyArray_SimpleNewFromData(2, dims,  NPY_DOUBLE, vec.data());
-    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");;
+    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");
     return array;
 }
 
@@ -41,7 +41,7 @@ PyObject* stdVec2ndarray(std::vector<double>& vec, int dim1, int dim2, int dim3)
     if (vec.size() != (unsigned) dim1 * dim2 * dim3) throw std::runtime_error("Size mismatch in stdVec2ndarray");
     npy_intp dims[3] = { dim1, dim2, dim3 };
     PyObject * array = PyArray_SimpleNewFromData(3, dims,  NPY_DOUBLE, vec.data());
-    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");;
+    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");
     return array;
 }
 
@@ -58,7 +58,7 @@ PyObject* stdVec2ndarray(std::vector<double>& vec, int dim1, int dim2, int dim3,
     if (vec.size() != (unsigned) dim1 * dim2 * dim3 * dim4) throw std::runtime_error("Size mismatch in stdVec2ndarray");
     npy_intp dims[4] = { dim1, dim2, dim3, dim4 };
     PyObject * array = PyArray_SimpleNewFromData(4, dims,  NPY_DOUBLE, vec.data());
-    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");;
+    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");
     return array;
 }
 
@@ -73,7 +73,7 @@ PyObject* stdVec2ndarray(std::vector<int>& vec, int dim1) {
     if (vec.size() != (unsigned) dim1) throw std::runtime_error("Size mismatch in stdVec2ndarray");
     npy_intp dims[1] = { dim1 };
     PyObject * array = PyArray_SimpleNewFromData(1, dims,  NPY_INT, vec.data());
-    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");;
+    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");
     return array;
 }
 
@@ -88,7 +88,7 @@ PyObject* stdVec2ndarray(std::vector<int>& vec, int dim1, int dim2) {
     if (vec.size() != (unsigned) dim1 * dim2) throw std::runtime_error("Size mismatch in stdVec2ndarray");
     npy_intp dims[2] = { dim1, dim2 };
     PyObject * array = PyArray_SimpleNewFromData(2, dims,  NPY_INT, vec.data());
-    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");;
+    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");
     return array;
 }
 
@@ -104,7 +104,7 @@ PyObject* stdVec2ndarray(std::vector<int>& vec, int dim1, int dim2, int dim3) {
     if (vec.size() != (unsigned) dim1 * dim2 * dim3) throw std::runtime_error("Size mismatch in stdVec2ndarray");
     npy_intp dims[3] = { dim1, dim2, dim3 };
     PyObject * array = PyArray_SimpleNewFromData(3, dims,  NPY_INT, vec.data());
-    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");;
+    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");
     return array;
 }
 
@@ -121,7 +121,7 @@ PyObject* stdVec2ndarray(std::vector<int>& vec, int dim1, int dim2, int dim3, in
     if (vec.size() != (unsigned) dim1 * dim2 * dim3 * dim4) throw std::runtime_error("Size mismatch in stdVec2ndarray");
     npy_intp dims[4] = { dim1, dim2, dim3, dim4 };
     PyObject * array = PyArray_SimpleNewFromData(4, dims,  NPY_INT, vec.data());
-    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");;
+    if (!array) throw std::runtime_error("Unknown failure in stdVec2ndarray");
     return array;
 }
 
diff --git a/tests/cpp/testfunctions.h b/tests/cpp/testfunctions.h
index dcd26b0747..ecb3ea933a 100644
--- a/tests/cpp/testfunctions.h
+++ b/tests/cpp/testfunctions.h
@@ -204,7 +204,7 @@ void verifyReturnData(const std::string &hdffile, const std::string &resultPath,
                       const ReturnData *rdata, const Model *model, double atol,
                       double rtol);
 
-void verifyReturnDataSensitivities(const H5::H5File &file_id,
+void verifyReturnDataSensitivities(const H5::H5File &file,
                                    const std::string &resultPath,
                                    const ReturnData *rdata, const Model *model,
                                    double atol, double rtol);