diff --git a/.github/actions/setup-swig/action.yml b/.github/actions/setup-swig/action.yml
index 0eb8a04473..b4441cd05e 100644
--- a/.github/actions/setup-swig/action.yml
+++ b/.github/actions/setup-swig/action.yml
@@ -7,7 +7,7 @@ inputs:
   swig_version:
     description: 'Swig version to build'
     required: false
-    default: '4.1.1'
+    default: '4.2.0'
 
 runs:
   using: "composite"
diff --git a/.github/workflows/deploy_branch.yml b/.github/workflows/deploy_branch.yml
index 82d1452d37..ff1be3145d 100644
--- a/.github/workflows/deploy_branch.yml
+++ b/.github/workflows/deploy_branch.yml
@@ -13,11 +13,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 20
 
@@ -31,7 +31,7 @@ jobs:
         scripts/buildSdist.sh
 
     - name: "Upload artifact: sdist"
-      uses: actions/upload-artifact@v3
+      uses: actions/upload-artifact@v4
       with:
         name: sdist
         path: python/sdist/dist/amici-*.gz
diff --git a/.github/workflows/deploy_protected.yml b/.github/workflows/deploy_protected.yml
index cd9257818d..c520896a70 100644
--- a/.github/workflows/deploy_protected.yml
+++ b/.github/workflows/deploy_protected.yml
@@ -25,10 +25,10 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git archive -o container/amici.tar.gz --format=tar.gz HEAD
     - name: Set up QEMU
       uses: docker/setup-qemu-action@v2
diff --git a/.github/workflows/deploy_release.yml b/.github/workflows/deploy_release.yml
index adcfd88d41..de7a073cc2 100644
--- a/.github/workflows/deploy_release.yml
+++ b/.github/workflows/deploy_release.yml
@@ -14,13 +14,20 @@ jobs:
       matrix:
         python-version: [3.9]
 
+    environment:
+      name: pypi
+      url: https://pypi.org/p/amici
+
+    permissions:
+      id-token: write
+
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 20
 
@@ -40,8 +47,6 @@ jobs:
     - name: Publish a Python distribution to PyPI
       uses: pypa/gh-action-pypi-publish@release/v1
       with:
-        user: __token__
-        password: ${{ secrets.pypi_password }}
         packages-dir: python/sdist/dist
 
   bioSimulatorsUpdateCliAndDockerImage:
diff --git a/.github/workflows/test_benchmark_collection_models.yml b/.github/workflows/test_benchmark_collection_models.yml
index bf9509e5ce..ab29938a12 100644
--- a/.github/workflows/test_benchmark_collection_models.yml
+++ b/.github/workflows/test_benchmark_collection_models.yml
@@ -29,11 +29,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 20
 
@@ -49,7 +49,7 @@ jobs:
     - name: Install AMICI sdist
       run: |
         pip3 install --user petab[vis] && \
-        AMICI_PARALLEL_COMPILE=2 pip3 install -v --user \
+        AMICI_PARALLEL_COMPILE="" pip3 install -v --user \
             $(ls -t python/sdist/dist/amici-*.tar.gz | head -1)[petab,test,vis]
 
     # retrieve test models
@@ -57,7 +57,7 @@ jobs:
       run: |
         git clone --depth 1 https://github.com/benchmarking-initiative/Benchmark-Models-PEtab.git \
           && export BENCHMARK_COLLECTION="$(pwd)/Benchmark-Models-PEtab/Benchmark-Models/" \
-          && AMICI_PARALLEL_COMPILE=2 tests/benchmark-models/test_benchmark_collection.sh
+          && AMICI_PARALLEL_COMPILE="" tests/benchmark-models/test_benchmark_collection.sh
 
     # run gradient checks
     - name: Run Gradient Checks
@@ -66,9 +66,9 @@ jobs:
           && cd tests/benchmark-models && pytest ./test_petab_benchmark.py
 
     # upload results
-    - uses: actions/upload-artifact@v3
+    - uses: actions/upload-artifact@v4
       with:
-        name: computation times
+        name: computation-times-${{ matrix.python-version }}-${{ matrix.extract_subexpressions }}
         path: |
           tests/benchmark-models/computation_times.csv
           tests/benchmark-models/computation_times.png
diff --git a/.github/workflows/test_doc.yml b/.github/workflows/test_doc.yml
index 7a60aa99d5..a4fe10de22 100644
--- a/.github/workflows/test_doc.yml
+++ b/.github/workflows/test_doc.yml
@@ -24,11 +24,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - name: Set up doxygen
@@ -47,11 +47,11 @@ jobs:
 
     steps:
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
 
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - run: git fetch --prune --unshallow
 
       - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV
diff --git a/.github/workflows/test_install.yml b/.github/workflows/test_install.yml
index 166616b434..0c624a33e5 100644
--- a/.github/workflows/test_install.yml
+++ b/.github/workflows/test_install.yml
@@ -13,11 +13,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV
@@ -54,11 +54,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV
@@ -87,11 +87,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV
diff --git a/.github/workflows/test_matlab.yml b/.github/workflows/test_matlab.yml
index 914ce19f09..46185e9d45 100644
--- a/.github/workflows/test_matlab.yml
+++ b/.github/workflows/test_matlab.yml
@@ -15,7 +15,7 @@ jobs:
     runs-on: ubuntu-20.04
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV
diff --git a/.github/workflows/test_performance.yml b/.github/workflows/test_performance.yml
index de9dd686d9..c02245224c 100644
--- a/.github/workflows/test_performance.yml
+++ b/.github/workflows/test_performance.yml
@@ -27,11 +27,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 20
 
@@ -51,7 +51,7 @@ jobs:
 
     - name: Install AMICI sdist
       run: |
-        AMICI_PARALLEL_COMPILE=2 check_time.sh \
+        AMICI_PARALLEL_COMPILE="" check_time.sh \
             install_sdist pip3 install -v --user \
             $(ls -t python/sdist/dist/amici-*.tar.gz | head -1)
 
@@ -64,7 +64,7 @@ jobs:
         AMICI_IMPORT_NPROCS=2 check_time.sh petab_import python tests/performance/test.py import
 
     - name: "Upload artifact: CS_Signalling_ERBB_RAS_AKT_petab"
-      uses: actions/upload-artifact@v3
+      uses: actions/upload-artifact@v4
       with:
         name: model_performance_test
         path: model_performance_test
diff --git a/.github/workflows/test_petab_test_suite.yml b/.github/workflows/test_petab_test_suite.yml
index 9b48f21d95..1af341009f 100644
--- a/.github/workflows/test_petab_test_suite.yml
+++ b/.github/workflows/test_petab_test_suite.yml
@@ -26,11 +26,11 @@ jobs:
 
     steps:
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
 
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
         with:
           fetch-depth: 20
 
@@ -78,15 +78,14 @@ jobs:
       - name: Run PEtab test suite
         run: |
           source ./build/venv/bin/activate \
-          && AMICI_PARALLEL_COMPILE=2 pytest -v \
+          && AMICI_PARALLEL_COMPILE="" pytest -v \
             --cov-report=xml:coverage.xml \
             --cov-append \
             --cov=amici \
             tests/petab_test_suite/
 
       - name: Codecov
-        if: github.event_name == 'pull_request' || github.repository_owner == 'AMICI-dev'
-        uses: codecov/codecov-action@v3.1.0
+        uses: codecov/codecov-action@v4
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: coverage.xml
diff --git a/.github/workflows/test_pypi.yml b/.github/workflows/test_pypi.yml
index 68675c578a..fecda3edb0 100644
--- a/.github/workflows/test_pypi.yml
+++ b/.github/workflows/test_pypi.yml
@@ -36,7 +36,7 @@ jobs:
           fi
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
 
diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml
index 5aa47173b8..507283b156 100644
--- a/.github/workflows/test_python_cplusplus.yml
+++ b/.github/workflows/test_python_cplusplus.yml
@@ -19,11 +19,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV
@@ -66,8 +66,7 @@ jobs:
           ${AMICI_DIR}/python/tests/test_splines.py
 
     - name: Codecov Python
-      if: github.event_name == 'pull_request' || github.repository_owner == 'AMICI-dev'
-      uses: codecov/codecov-action@v3.1.0
+      uses: codecov/codecov-action@v4
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         file: build/coverage_py.xml
@@ -86,8 +85,7 @@ jobs:
         && lcov -a coverage_cpp.info -a coverage_py.info -o coverage.info
 
     - name: Codecov CPP
-      if: github.event_name == 'pull_request' || github.repository_owner == 'AMICI-dev'
-      uses: codecov/codecov-action@v3.1.0
+      uses: codecov/codecov-action@v4
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         file: coverage.info
@@ -114,11 +112,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV
@@ -141,8 +139,7 @@ jobs:
           ${AMICI_DIR}/python/tests/test_splines_short.py
 
     - name: Codecov Python
-      if: github.event_name == 'pull_request' || github.repository_owner == 'AMICI-dev'
-      uses: codecov/codecov-action@v3.1.0
+      uses: codecov/codecov-action@v4
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         file: build/coverage_py.xml
@@ -161,8 +158,7 @@ jobs:
         && lcov -a coverage_cpp.info -a coverage_py.info -o coverage.info
 
     - name: Codecov CPP
-      if: github.event_name == 'pull_request' || github.repository_owner == 'AMICI-dev'
-      uses: codecov/codecov-action@v3.1.0
+      uses: codecov/codecov-action@v4
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         file: coverage.info
@@ -190,11 +186,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV
@@ -224,11 +220,11 @@ jobs:
 
     steps:
     - name: Set up Python
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: 3.9
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - name: Install dependencies
@@ -266,11 +262,11 @@ jobs:
 
     steps:
     - name: Set up Python
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: 3.9
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - name: Install dependencies
diff --git a/.github/workflows/test_python_ver_matrix.yml b/.github/workflows/test_python_ver_matrix.yml
index 9290cd0c1a..5f9d96d295 100644
--- a/.github/workflows/test_python_ver_matrix.yml
+++ b/.github/workflows/test_python_ver_matrix.yml
@@ -33,11 +33,11 @@ jobs:
     - run: echo "BNGPATH=${AMICI_DIR}/ThirdParty/BioNetGen-2.7.0" >> $GITHUB_ENV
 
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 20
 
diff --git a/.github/workflows/test_sbml_semantic_test_suite.yml b/.github/workflows/test_sbml_semantic_test_suite.yml
index 2af34d3762..cf03a5d458 100644
--- a/.github/workflows/test_sbml_semantic_test_suite.yml
+++ b/.github/workflows/test_sbml_semantic_test_suite.yml
@@ -33,29 +33,28 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 1
 
     - name: Install apt dependencies
       uses: ./.github/actions/install-apt-dependencies
 
-    - run: AMICI_PARALLEL_COMPILE=2 ./scripts/installAmiciSource.sh
-    - run: AMICI_PARALLEL_COMPILE=2 ./scripts/run-SBMLTestsuite.sh ${{ matrix.cases }}
+    - run: AMICI_PARALLEL_COMPILE="" ./scripts/installAmiciSource.sh
+    - run: AMICI_PARALLEL_COMPILE="" ./scripts/run-SBMLTestsuite.sh ${{ matrix.cases }}
 
     - name: "Upload artifact: SBML semantic test suite results"
-      uses: actions/upload-artifact@v3
+      uses: actions/upload-artifact@v4
       with:
-        name: amici-semantic-results
+        name: amici-semantic-results-${{ matrix.cases }}
         path: tests/amici-semantic-results
 
     - name: Codecov SBMLSuite
-      if: github.event_name == 'pull_request' || github.repository_owner == 'AMICI-dev'
-      uses: codecov/codecov-action@v3.1.0
+      uses: codecov/codecov-action@v4
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         file: coverage_SBMLSuite.xml
diff --git a/.github/workflows/test_valgrind.yml b/.github/workflows/test_valgrind.yml
index 32eea2e0c0..034945d461 100644
--- a/.github/workflows/test_valgrind.yml
+++ b/.github/workflows/test_valgrind.yml
@@ -25,11 +25,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - name: Install apt dependencies
@@ -64,11 +64,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - name: Install apt dependencies
diff --git a/.github/workflows/test_windows.yml b/.github/workflows/test_windows.yml
index 9dc13efbea..fbd91ce73a 100644
--- a/.github/workflows/test_windows.yml
+++ b/.github/workflows/test_windows.yml
@@ -31,11 +31,11 @@ jobs:
 
     steps:
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - run: git fetch --prune --unshallow
 
     - shell: bash
@@ -61,7 +61,7 @@ jobs:
       shell: powershell
       run: scripts/installOpenBLAS
 
-    - uses: actions/upload-artifact@v3
+    - uses: actions/upload-artifact@v4
       with:
         name: OpenBLAS
         path: C:\BLAS\OpenBLAS
diff --git a/.readthedocs.yml b/.readthedocs.yml
index 38c2e8e41b..23cc6addeb 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -17,7 +17,6 @@ formats:
 python:
   install:
     - requirements: documentation/rtd_requirements.txt
-    - requirements: documentation/rtd_requirements2.txt
 build:
   os: "ubuntu-22.04"
   apt_packages:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1fa8153e35..45edeb8adf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,18 @@
 
 ## v0.X Series
 
+### v0.21.2 (2024-02-06)
+
+* Fixed `Solver` copyctor issues with swig4.2 that resulted in installation
+  errors
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2276
+* Fixed error when calling `amici.ExpData()`
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2280
+* Fixed invalid-type-error when loading an antimony model from file
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2281
+
+**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.21.1...v0.21.2
+
 ### v0.21.1 (2024-01-17)
 
 Fixed package configuration for PyPI upload. No further changes.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d6f08a5097..8b6c9140e9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -130,8 +130,9 @@ set(VENDORED_SUNDIALS_DIR ${CMAKE_CURRENT_SOURCE_DIR}/ThirdParty/sundials)
 set(VENDORED_SUNDIALS_BUILD_DIR ${VENDORED_SUNDIALS_DIR}/build)
 set(VENDORED_SUNDIALS_INSTALL_DIR ${VENDORED_SUNDIALS_BUILD_DIR})
 set(SUNDIALS_PRIVATE_INCLUDE_DIRS "${VENDORED_SUNDIALS_DIR}/src")
-find_package(SUNDIALS REQUIRED PATHS
-             "${VENDORED_SUNDIALS_INSTALL_DIR}/${CMAKE_INSTALL_LIBDIR}/cmake/sundials/")
+find_package(
+  SUNDIALS REQUIRED PATHS
+  "${VENDORED_SUNDIALS_INSTALL_DIR}/${CMAKE_INSTALL_LIBDIR}/cmake/sundials/")
 message(STATUS "Found SUNDIALS: ${SUNDIALS_DIR}")
 
 set(GSL_LITE_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/ThirdParty/gsl")
diff --git a/LICENSE.md b/LICENSE.md
index 4dadbc28a9..a086d45955 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -5,7 +5,7 @@
 AMICI is released under the 3-Clause BSD License (BSD-3-Clause) with the
 following terms:
 
-Copyright (c) 2015-2020, Fabian Fröhlich, Jan Hasenauer, Daniel Weindl and Paul Stapor
+Copyright (c) 2015-2024, AMICI Developers.
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
diff --git a/documentation/rtd_requirements.txt b/documentation/rtd_requirements.txt
index ae49f9dfbb..8d2c2100f9 100644
--- a/documentation/rtd_requirements.txt
+++ b/documentation/rtd_requirements.txt
@@ -13,7 +13,7 @@ sphinx-autodoc-typehints
 git+https://github.com/readthedocs/sphinx-hoverxref@main
 ipython>=8.13.2
 breathe>=4.35.0
-#exhale>=0.3.5
+exhale>=0.3.7
 -e git+https://github.com/mithro/sphinx-contrib-mithro#egg=sphinx-contrib-exhale-multiproject&subdirectory=sphinx-contrib-exhale-multiproject
 sphinxcontrib-matlabdomain>=0.20.0
 sphinxcontrib-napoleon>=0.7
diff --git a/documentation/rtd_requirements2.txt b/documentation/rtd_requirements2.txt
deleted file mode 100644
index 2c307f0ea4..0000000000
--- a/documentation/rtd_requirements2.txt
+++ /dev/null
@@ -1 +0,0 @@
--e git+https://github.com/svenevs/exhale.git@a1a8551321e246e3ab81f5456e04a8159804595b#egg=exhale
diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h
index bb824577b9..d88347078d 100644
--- a/include/amici/abstract_model.h
+++ b/include/amici/abstract_model.h
@@ -41,7 +41,7 @@ class AbstractModel {
      * @param root array to which values of the root function will be written
      */
     virtual void froot(
-        const realtype t, AmiVector const& x, AmiVector const& dx,
+        realtype const t, AmiVector const& x, AmiVector const& dx,
         gsl::span<realtype> root
     ) = 0;
 
@@ -54,7 +54,7 @@ class AbstractModel {
      * written
      */
     virtual void fxdot(
-        const realtype t, AmiVector const& x, AmiVector const& dx,
+        realtype const t, AmiVector const& x, AmiVector const& dx,
         AmiVector& xdot
     ) = 0;
 
@@ -70,7 +70,7 @@ class AbstractModel {
      * will be written
      */
     virtual void fsxdot(
-        const realtype t, AmiVector const& x, AmiVector const& dx, int ip,
+        realtype const t, AmiVector const& x, AmiVector const& dx, int ip,
         AmiVector const& sx, AmiVector const& sdx, AmiVector& sxdot
     ) = 0;
 
@@ -83,7 +83,7 @@ class AbstractModel {
      * written
      */
     virtual void fxBdot_ss(
-        const realtype t, AmiVector const& xB, AmiVector const& dxB,
+        realtype const t, AmiVector const& xB, AmiVector const& dxB,
         AmiVector& xBdot
     ) = 0;
 
@@ -105,7 +105,7 @@ class AbstractModel {
      * @param xBdot Vector with the adjoint state right hand side
      */
     virtual void writeSteadystateJB(
-        const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+        realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot
     ) = 0;
 
@@ -119,7 +119,7 @@ class AbstractModel {
      * @param J dense matrix to which values of the jacobian will be written
      */
     virtual void
-    fJ(const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+    fJ(realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
        AmiVector const& xdot, SUNMatrix J)
         = 0;
 
@@ -135,7 +135,7 @@ class AbstractModel {
      * @param JB dense matrix to which values of the jacobian will be written
      */
     virtual void
-    fJB(const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+    fJB(realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot,
         SUNMatrix JB)
         = 0;
@@ -150,7 +150,7 @@ class AbstractModel {
      * @param J sparse matrix to which values of the Jacobian will be written
      */
     virtual void fJSparse(
-        const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+        realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xdot, SUNMatrix J
     ) = 0;
 
@@ -166,7 +166,7 @@ class AbstractModel {
      * @param JB dense matrix to which values of the jacobian will be written
      */
     virtual void fJSparseB(
-        const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+        realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot,
         SUNMatrix JB
     ) = 0;
@@ -180,7 +180,7 @@ class AbstractModel {
      * @param dx time derivative of state (DAE only)
      */
     virtual void fJDiag(
-        const realtype t, AmiVector& Jdiag, realtype cj, AmiVector const& x,
+        realtype const t, AmiVector& Jdiag, realtype cj, AmiVector const& x,
         AmiVector const& dx
     ) = 0;
 
@@ -192,7 +192,7 @@ class AbstractModel {
      * @param dx time derivative of state (DAE only)
      */
     virtual void
-    fdxdotdp(const realtype t, AmiVector const& x, AmiVector const& dx)
+    fdxdotdp(realtype const t, AmiVector const& x, AmiVector const& dx)
         = 0;
 
     /**
@@ -206,7 +206,7 @@ class AbstractModel {
      * @param cj scaling factor (inverse of timestep, DAE only)
      */
     virtual void
-    fJv(const realtype t, AmiVector const& x, AmiVector const& dx,
+    fJv(realtype const t, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xdot, AmiVector const& v, AmiVector& nJv, realtype cj)
         = 0;
 
@@ -230,7 +230,7 @@ class AbstractModel {
      * @param k constant vector
      */
     virtual void
-    fx0(realtype* x0, const realtype t, realtype const* p, realtype const* k);
+    fx0(realtype* x0, realtype const t, realtype const* p, realtype const* k);
 
     /**
      * @brief Function indicating whether reinitialization of states depending
@@ -250,7 +250,7 @@ class AbstractModel {
      * based on provided constants / fixed parameters.
      */
     virtual void fx0_fixedParameters(
-        realtype* x0, const realtype t, realtype const* p, realtype const* k,
+        realtype* x0, realtype const t, realtype const* p, realtype const* k,
         gsl::span<int const> reinitialization_state_idxs
     );
 
@@ -266,7 +266,7 @@ class AbstractModel {
      * based on provided constants / fixed parameters.
      */
     virtual void fsx0_fixedParameters(
-        realtype* sx0, const realtype t, realtype const* x0, realtype const* p,
+        realtype* sx0, realtype const t, realtype const* x0, realtype const* p,
         realtype const* k, int ip,
         gsl::span<int const> reinitialization_state_idxs
     );
@@ -281,7 +281,7 @@ class AbstractModel {
      * @param ip sensitivity index
      */
     virtual void fsx0(
-        realtype* sx0, const realtype t, realtype const* x0, realtype const* p,
+        realtype* sx0, realtype const t, realtype const* x0, realtype const* p,
         realtype const* k, int ip
     );
 
@@ -308,7 +308,7 @@ class AbstractModel {
      * @param ie event index
      */
     virtual void fstau(
-        realtype* stau, const realtype t, realtype const* x, realtype const* p,
+        realtype* stau, realtype const t, realtype const* x, realtype const* p,
         realtype const* k, realtype const* h, realtype const* tcl,
         realtype const* sx, int ip, int ie
     );
@@ -324,7 +324,7 @@ class AbstractModel {
      * @param w repeating elements vector
      */
     virtual void
-    fy(realtype* y, const realtype t, realtype const* x, realtype const* p,
+    fy(realtype* y, realtype const t, realtype const* x, realtype const* p,
        realtype const* k, realtype const* h, realtype const* w);
 
     /**
@@ -340,7 +340,7 @@ class AbstractModel {
      * @param dwdp Recurring terms in xdot, parameter derivative
      */
     virtual void fdydp(
-        realtype* dydp, const realtype t, realtype const* x, realtype const* p,
+        realtype* dydp, realtype const t, realtype const* x, realtype const* p,
         realtype const* k, realtype const* h, int ip, realtype const* w,
         realtype const* dwdp
     );
@@ -362,7 +362,7 @@ class AbstractModel {
      * \f$
      */
     virtual void fdydp(
-        realtype* dydp, const realtype t, realtype const* x, realtype const* p,
+        realtype* dydp, realtype const t, realtype const* x, realtype const* p,
         realtype const* k, realtype const* h, int ip, realtype const* w,
         realtype const* tcl, realtype const* dtcldp, realtype const* spl,
         realtype const* sspl
@@ -380,7 +380,7 @@ class AbstractModel {
      * @param dwdx Recurring terms in xdot, state derivative
      */
     virtual void fdydx(
-        realtype* dydx, const realtype t, realtype const* x, realtype const* p,
+        realtype* dydx, realtype const t, realtype const* x, realtype const* p,
         realtype const* k, realtype const* h, realtype const* w,
         realtype const* dwdx
     );
@@ -396,7 +396,7 @@ class AbstractModel {
      * @param h Heaviside vector
      */
     virtual void
-    fz(realtype* z, int ie, const realtype t, realtype const* x,
+    fz(realtype* z, int ie, realtype const t, realtype const* x,
        realtype const* p, realtype const* k, realtype const* h);
 
     /**
@@ -412,7 +412,7 @@ class AbstractModel {
      * @param ip sensitivity index
      */
     virtual void
-    fsz(realtype* sz, int ie, const realtype t, realtype const* x,
+    fsz(realtype* sz, int ie, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h,
         realtype const* sx, int ip);
 
@@ -428,7 +428,7 @@ class AbstractModel {
      * @param h Heaviside vector
      */
     virtual void
-    frz(realtype* rz, int ie, const realtype t, realtype const* x,
+    frz(realtype* rz, int ie, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h);
 
     /**
@@ -444,7 +444,7 @@ class AbstractModel {
      * @param ip sensitivity index
      */
     virtual void fsrz(
-        realtype* srz, int ie, const realtype t, realtype const* x,
+        realtype* srz, int ie, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h,
         realtype const* sx, int ip
     );
@@ -462,7 +462,7 @@ class AbstractModel {
      * @param ip parameter index w.r.t. which the derivative is requested
      */
     virtual void fdzdp(
-        realtype* dzdp, int ie, const realtype t, realtype const* x,
+        realtype* dzdp, int ie, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h, int ip
     );
 
@@ -478,7 +478,7 @@ class AbstractModel {
      * @param h Heaviside vector
      */
     virtual void fdzdx(
-        realtype* dzdx, int ie, const realtype t, realtype const* x,
+        realtype* dzdx, int ie, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h
     );
 
@@ -495,7 +495,7 @@ class AbstractModel {
      * @param ip parameter index w.r.t. which the derivative is requested
      */
     virtual void fdrzdp(
-        realtype* drzdp, int ie, const realtype t, realtype const* x,
+        realtype* drzdp, int ie, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h, int ip
     );
 
@@ -510,7 +510,7 @@ class AbstractModel {
      * @param h Heaviside vector
      */
     virtual void fdrzdx(
-        realtype* drzdx, int ie, const realtype t, realtype const* x,
+        realtype* drzdx, int ie, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h
     );
 
@@ -527,7 +527,7 @@ class AbstractModel {
      * @param xdot_old previous model right hand side
      */
     virtual void fdeltax(
-        realtype* deltax, const realtype t, realtype const* x,
+        realtype* deltax, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h, int ie,
         realtype const* xdot, realtype const* xdot_old
     );
@@ -550,7 +550,7 @@ class AbstractModel {
      * @param tcl total abundances for conservation laws
      */
     virtual void fdeltasx(
-        realtype* deltasx, const realtype t, realtype const* x,
+        realtype* deltasx, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h,
         realtype const* w, int ip, int ie, realtype const* xdot,
         realtype const* xdot_old, realtype const* sx, realtype const* stau,
@@ -571,7 +571,7 @@ class AbstractModel {
      * @param xB current adjoint state
      */
     virtual void fdeltaxB(
-        realtype* deltaxB, const realtype t, realtype const* x,
+        realtype* deltaxB, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h, int ie,
         realtype const* xdot, realtype const* xdot_old, realtype const* xB
     );
@@ -591,7 +591,7 @@ class AbstractModel {
      * @param xB adjoint state
      */
     virtual void fdeltaqB(
-        realtype* deltaqB, const realtype t, realtype const* x,
+        realtype* deltaqB, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h, int ip, int ie,
         realtype const* xdot, realtype const* xdot_old, realtype const* xB
     );
@@ -605,7 +605,7 @@ class AbstractModel {
      * @param y model output at timepoint t
      */
     virtual void fsigmay(
-        realtype* sigmay, const realtype t, realtype const* p,
+        realtype* sigmay, realtype const t, realtype const* p,
         realtype const* k, realtype const* y
     );
 
@@ -619,7 +619,7 @@ class AbstractModel {
      * @param ip sensitivity index
      */
     virtual void fdsigmaydp(
-        realtype* dsigmaydp, const realtype t, realtype const* p,
+        realtype* dsigmaydp, realtype const t, realtype const* p,
         realtype const* k, realtype const* y, int ip
     );
     /**
@@ -632,7 +632,7 @@ class AbstractModel {
      * @param y model output at timepoint t
      */
     virtual void fdsigmaydy(
-        realtype* dsigmaydy, const realtype t, realtype const* p,
+        realtype* dsigmaydy, realtype const t, realtype const* p,
         realtype const* k, realtype const* y
     );
 
@@ -644,7 +644,7 @@ class AbstractModel {
      * @param k constant vector
      */
     virtual void fsigmaz(
-        realtype* sigmaz, const realtype t, realtype const* p, realtype const* k
+        realtype* sigmaz, realtype const t, realtype const* p, realtype const* k
     );
 
     /**
@@ -657,7 +657,7 @@ class AbstractModel {
      * @param ip sensitivity index
      */
     virtual void fdsigmazdp(
-        realtype* dsigmazdp, const realtype t, realtype const* p,
+        realtype* dsigmazdp, realtype const t, realtype const* p,
         realtype const* k, int ip
     );
 
@@ -822,7 +822,7 @@ class AbstractModel {
      * @param spl spline value vector
      */
     virtual void
-    fw(realtype* w, const realtype t, realtype const* x, realtype const* p,
+    fw(realtype* w, realtype const t, realtype const* x, realtype const* p,
        realtype const* k, realtype const* h, realtype const* tcl,
        realtype const* spl);
 
@@ -842,7 +842,7 @@ class AbstractModel {
      * \f$
      */
     virtual void fdwdp(
-        realtype* dwdp, const realtype t, realtype const* x, realtype const* p,
+        realtype* dwdp, realtype const t, realtype const* x, realtype const* p,
         realtype const* k, realtype const* h, realtype const* w,
         realtype const* tcl, realtype const* stcl, realtype const* spl,
         realtype const* sspl
@@ -876,7 +876,7 @@ class AbstractModel {
      * @param ip sensitivity parameter index
      */
     virtual void fdwdp(
-        realtype* dwdp, const realtype t, realtype const* x, realtype const* p,
+        realtype* dwdp, realtype const t, realtype const* x, realtype const* p,
         realtype const* k, realtype const* h, realtype const* w,
         realtype const* tcl, realtype const* stcl, realtype const* spl,
         realtype const* sspl, int ip
@@ -895,7 +895,7 @@ class AbstractModel {
      * @param spl spline value vector
      */
     virtual void fdwdx(
-        realtype* dwdx, const realtype t, realtype const* x, realtype const* p,
+        realtype* dwdx, realtype const t, realtype const* x, realtype const* p,
         realtype const* k, realtype const* h, realtype const* w,
         realtype const* tcl, realtype const* spl
     );
diff --git a/include/amici/backwardproblem.h b/include/amici/backwardproblem.h
index 1c26186c17..a5d32e3c47 100644
--- a/include/amici/backwardproblem.h
+++ b/include/amici/backwardproblem.h
@@ -134,7 +134,7 @@ class BackwardProblem {
     /** state derivative of data likelihood */
     std::vector<realtype> dJydx_;
     /** state derivative of event likelihood */
-    const std::vector<realtype> dJzdx_;
+    std::vector<realtype> const dJzdx_;
 };
 
 } // namespace amici
diff --git a/include/amici/misc.h b/include/amici/misc.h
index 6dc3294240..6874c19185 100644
--- a/include/amici/misc.h
+++ b/include/amici/misc.h
@@ -85,7 +85,7 @@ void checkBufferSize(
  * @param buffer buffer to which values are to be written
  */
 template <class T>
-void writeSlice(const gsl::span<T const> slice, gsl::span<T> buffer) {
+void writeSlice(gsl::span<T const> const slice, gsl::span<T> buffer) {
     checkBufferSize(buffer, slice.size());
     std::copy(slice.begin(), slice.end(), buffer.data());
 };
@@ -97,7 +97,7 @@ void writeSlice(const gsl::span<T const> slice, gsl::span<T> buffer) {
  * @param buffer buffer to which values are to be added
  */
 template <class T>
-void addSlice(const gsl::span<T const> slice, gsl::span<T> buffer) {
+void addSlice(gsl::span<T const> const slice, gsl::span<T> buffer) {
     checkBufferSize(buffer, slice.size());
     std::transform(
         slice.begin(), slice.end(), buffer.begin(), buffer.begin(),
diff --git a/include/amici/model.h b/include/amici/model.h
index 481164afdf..23f6ac63c8 100644
--- a/include/amici/model.h
+++ b/include/amici/model.h
@@ -92,7 +92,7 @@ enum class ModelQuantity {
     drzdx,
 };
 
-extern const std::map<ModelQuantity, std::string> model_quantity_to_str;
+extern std::map<ModelQuantity, std::string> const model_quantity_to_str;
 
 /**
  * @brief The Model class represents an AMICI ODE/DAE model.
@@ -918,7 +918,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param x Current state
      */
     void
-    getExpression(gsl::span<realtype> w, const realtype t, AmiVector const& x);
+    getExpression(gsl::span<realtype> w, realtype const t, AmiVector const& x);
 
     /**
      * @brief Get time-resolved observables.
@@ -927,7 +927,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param x Current state
      */
     void
-    getObservable(gsl::span<realtype> y, const realtype t, AmiVector const& x);
+    getObservable(gsl::span<realtype> y, realtype const t, AmiVector const& x);
 
     /**
      * @brief Get scaling type for observable
@@ -947,7 +947,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param sx State sensitivities
      */
     void getObservableSensitivity(
-        gsl::span<realtype> sy, const realtype t, AmiVector const& x,
+        gsl::span<realtype> sy, realtype const t, AmiVector const& x,
         AmiVectorArray const& sx
     );
 
@@ -1045,7 +1045,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param x State variables
      */
     void getEvent(
-        gsl::span<realtype> z, int const ie, const realtype t,
+        gsl::span<realtype> z, int const ie, realtype const t,
         AmiVector const& x
     );
     /**
@@ -1060,7 +1060,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param sx State sensitivities
      */
     void getEventSensitivity(
-        gsl::span<realtype> sz, int const ie, const realtype t,
+        gsl::span<realtype> sz, int const ie, realtype const t,
         AmiVector const& x, AmiVectorArray const& sx
     );
 
@@ -1082,7 +1082,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param x State variables
      */
     void getEventRegularization(
-        gsl::span<realtype> rz, int const ie, const realtype t,
+        gsl::span<realtype> rz, int const ie, realtype const t,
         AmiVector const& x
     );
 
@@ -1099,7 +1099,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param sx State sensitivities
      */
     void getEventRegularizationSensitivity(
-        gsl::span<realtype> srz, int const ie, const realtype t,
+        gsl::span<realtype> srz, int const ie, realtype const t,
         AmiVector const& x, AmiVectorArray const& sx
     );
     /**
@@ -1113,7 +1113,7 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void getEventSigma(
         gsl::span<realtype> sigmaz, int const ie, int const nroots,
-        const realtype t, ExpData const* edata
+        realtype const t, ExpData const* edata
     );
 
     /**
@@ -1131,7 +1131,7 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void getEventSigmaSensitivity(
         gsl::span<realtype> ssigmaz, int const ie, int const nroots,
-        const realtype t, ExpData const* edata
+        realtype const t, ExpData const* edata
     );
 
     /**
@@ -1144,7 +1144,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param edata Experimental data
      */
     void addEventObjective(
-        realtype& Jz, int const ie, int const nroots, const realtype t,
+        realtype& Jz, int const ie, int const nroots, realtype const t,
         AmiVector const& x, ExpData const& edata
     );
 
@@ -1158,7 +1158,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param edata Experimental data
      */
     void addEventObjectiveRegularization(
-        realtype& Jrz, int const ie, int const nroots, const realtype t,
+        realtype& Jrz, int const ie, int const nroots, realtype const t,
         AmiVector const& x, ExpData const& edata
     );
 
@@ -1180,7 +1180,7 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void addEventObjectiveSensitivity(
         std::vector<realtype>& sllh, std::vector<realtype>& s2llh, int const ie,
-        int const nroots, const realtype t, AmiVector const& x,
+        int const nroots, realtype const t, AmiVector const& x,
         AmiVectorArray const& sx, ExpData const& edata
     );
 
@@ -1201,7 +1201,7 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void addPartialEventObjectiveSensitivity(
         std::vector<realtype>& sllh, std::vector<realtype>& s2llh, int const ie,
-        int const nroots, const realtype t, AmiVector const& x,
+        int const nroots, realtype const t, AmiVector const& x,
         ExpData const& edata
     );
 
@@ -1219,7 +1219,7 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void getAdjointStateEventUpdate(
         gsl::span<realtype> dJzdx, int const ie, int const nroots,
-        const realtype t, AmiVector const& x, ExpData const& edata
+        realtype const t, AmiVector const& x, ExpData const& edata
     );
 
     /**
@@ -1234,7 +1234,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param sx State sensitivities
      */
     void getEventTimeSensitivity(
-        std::vector<realtype>& stau, const realtype t, int const ie,
+        std::vector<realtype>& stau, realtype const t, int const ie,
         AmiVector const& x, AmiVectorArray const& sx
     );
 
@@ -1247,7 +1247,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param xdot_old Value of residual function before event
      */
     void addStateEventUpdate(
-        AmiVector& x, int const ie, const realtype t, AmiVector const& xdot,
+        AmiVector& x, int const ie, realtype const t, AmiVector const& xdot,
         AmiVector const& xdot_old
     );
 
@@ -1263,7 +1263,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * `Model::getEventTimeSensitivity`
      */
     void addStateSensitivityEventUpdate(
-        AmiVectorArray& sx, int const ie, const realtype t,
+        AmiVectorArray& sx, int const ie, realtype const t,
         AmiVector const& x_old, AmiVector const& xdot,
         AmiVector const& xdot_old, std::vector<realtype> const& stau
     );
@@ -1278,7 +1278,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param xdot_old Value of residual function before event
      */
     void addAdjointStateEventUpdate(
-        AmiVector& xB, int const ie, const realtype t, AmiVector const& x,
+        AmiVector& xB, int const ie, realtype const t, AmiVector const& x,
         AmiVector const& xdot, AmiVector const& xdot_old
     );
 
@@ -1293,7 +1293,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param xdot_old Value of residual function before event
      */
     void addAdjointQuadratureEventUpdate(
-        AmiVector xQB, int const ie, const realtype t, AmiVector const& x,
+        AmiVector xQB, int const ie, realtype const t, AmiVector const& x,
         AmiVector const& xB, AmiVector const& xdot, AmiVector const& xdot_old
     );
 
@@ -1693,7 +1693,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param edata Experimental data
      */
     void fsigmaz(
-        int const ie, int const nroots, const realtype t, ExpData const* edata
+        int const ie, int const nroots, realtype const t, ExpData const* edata
     );
 
     /**
@@ -1727,7 +1727,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param edata Experimental data
      */
     void fdJzdz(
-        int const ie, int const nroots, const realtype t, AmiVector const& x,
+        int const ie, int const nroots, realtype const t, AmiVector const& x,
         ExpData const& edata
     );
 
@@ -1741,7 +1741,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param edata Pointer to experimental data instance
      */
     void fdJzdsigma(
-        int const ie, int const nroots, const realtype t, AmiVector const& x,
+        int const ie, int const nroots, realtype const t, AmiVector const& x,
         ExpData const& edata
     );
 
@@ -1794,7 +1794,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param edata Experimental data
      */
     void fdJrzdz(
-        int const ie, int const nroots, const realtype t, AmiVector const& x,
+        int const ie, int const nroots, realtype const t, AmiVector const& x,
         ExpData const& edata
     );
 
@@ -1808,7 +1808,7 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param edata pointer to experimental data instance
      */
     void fdJrzdsigma(
-        int const ie, int const nroots, const realtype t, AmiVector const& x,
+        int const ie, int const nroots, realtype const t, AmiVector const& x,
         ExpData const& edata
     );
 
@@ -2016,11 +2016,13 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /** method for steady-state computation */
     SteadyStateComputationMode steadystate_computation_mode_{
-        SteadyStateComputationMode::integrateIfNewtonFails};
+        SteadyStateComputationMode::integrateIfNewtonFails
+    };
 
     /** method for steadystate sensitivities computation */
     SteadyStateSensitivityMode steadystate_sensitivity_mode_{
-        SteadyStateSensitivityMode::integrateIfNewtonFails};
+        SteadyStateSensitivityMode::integrateIfNewtonFails
+    };
 
     /**
      * Indicates whether the result of every call to `Model::f*` should be
diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h
index b35cfa1d70..232bfc867d 100644
--- a/include/amici/model_dae.h
+++ b/include/amici/model_dae.h
@@ -46,7 +46,7 @@ class Model_DAE : public Model {
     Model_DAE(
         ModelDimensions const& model_dimensions,
         SimulationParameters simulation_parameters,
-        const SecondOrderMode o2mode, std::vector<realtype> const& idlist,
+        SecondOrderMode const o2mode, std::vector<realtype> const& idlist,
         std::vector<int> const& z2event, bool const pythonGenerated = false,
         int const ndxdotdp_explicit = 0, int const ndxdotdx_explicit = 0,
         int const w_recursion_depth = 0,
@@ -85,7 +85,7 @@ class Model_DAE : public Model {
        const_N_Vector xdot, SUNMatrix J);
 
     void
-    fJB(const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+    fJB(realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot,
         SUNMatrix JB) override;
 
@@ -122,7 +122,7 @@ class Model_DAE : public Model {
     );
 
     void fJSparseB(
-        const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+        realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot,
         SUNMatrix JB
     ) override;
@@ -253,7 +253,7 @@ class Model_DAE : public Model {
     );
 
     void fxBdot_ss(
-        const realtype t, AmiVector const& xB, AmiVector const& dxB,
+        realtype const t, AmiVector const& xB, AmiVector const& dxB,
         AmiVector& xBdot
     ) override;
 
@@ -298,7 +298,7 @@ class Model_DAE : public Model {
      * @param xBdot Vector with the adjoint state right hand side
      */
     void writeSteadystateJB(
-        const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+        realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot
     ) override;
 
@@ -308,8 +308,8 @@ class Model_DAE : public Model {
      * @param x Vector with the states
      * @param dx Vector with the derivative states
      */
-    void fdxdotdp(realtype t, const const_N_Vector x, const const_N_Vector dx);
-    void fdxdotdp(const realtype t, AmiVector const& x, AmiVector const& dx)
+    void fdxdotdp(realtype t, const_N_Vector const x, const_N_Vector const dx);
+    void fdxdotdp(realtype const t, AmiVector const& x, AmiVector const& dx)
         override {
         fdxdotdp(t, x.getNVector(), dx.getNVector());
     };
@@ -526,7 +526,7 @@ class Model_DAE : public Model {
      * @param k constants vector
      */
     virtual void
-    fM(realtype* M, const realtype t, realtype const* x, realtype const* p,
+    fM(realtype* M, realtype const t, realtype const* x, realtype const* p,
        realtype const* k);
 };
 } // namespace amici
diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h
index e03e1867d1..3632c2198f 100644
--- a/include/amici/model_ode.h
+++ b/include/amici/model_ode.h
@@ -45,7 +45,7 @@ class Model_ODE : public Model {
     Model_ODE(
         ModelDimensions const& model_dimensions,
         SimulationParameters simulation_parameters,
-        const SecondOrderMode o2mode, std::vector<realtype> const& idlist,
+        SecondOrderMode const o2mode, std::vector<realtype> const& idlist,
         std::vector<int> const& z2event, bool const pythonGenerated = false,
         int const ndxdotdp_explicit = 0, int const ndxdotdx_explicit = 0,
         int const w_recursion_depth = 0,
@@ -75,7 +75,7 @@ class Model_ODE : public Model {
     void fJ(realtype t, const_N_Vector x, const_N_Vector xdot, SUNMatrix J);
 
     void
-    fJB(const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+    fJB(realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot,
         SUNMatrix JB) override;
 
@@ -109,7 +109,7 @@ class Model_ODE : public Model {
     void fJSparse(realtype t, const_N_Vector x, SUNMatrix J);
 
     void fJSparseB(
-        const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+        realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot,
         SUNMatrix JB
     ) override;
@@ -229,7 +229,7 @@ class Model_ODE : public Model {
     fqBdot(realtype t, const_N_Vector x, const_N_Vector xB, N_Vector qBdot);
 
     void fxBdot_ss(
-        const realtype t, AmiVector const& xB, AmiVector const& /*dxB*/,
+        realtype const t, AmiVector const& xB, AmiVector const& /*dxB*/,
         AmiVector& xBdot
     ) override;
 
@@ -268,7 +268,7 @@ class Model_ODE : public Model {
      * @param xBdot Vector with the adjoint state right hand side
      */
     void writeSteadystateJB(
-        const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+        realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
         AmiVector const& xB, AmiVector const& dxB, AmiVector const& xBdot
     ) override;
 
diff --git a/include/amici/rdata.h b/include/amici/rdata.h
index 3b9ea01b9f..1de02c99db 100644
--- a/include/amici/rdata.h
+++ b/include/amici/rdata.h
@@ -688,7 +688,7 @@ class ReturnData : public ModelDimensions {
      * @param edata ExpData instance carrying experimental data
      */
     void getEventOutput(
-        realtype t, const std::vector<int> rootidx, Model& model,
+        realtype t, std::vector<int> const rootidx, Model& model,
         ExpData const* edata
     );
 
diff --git a/include/amici/returndata_matlab.h b/include/amici/returndata_matlab.h
index 6a51db5264..26b7d55d9e 100644
--- a/include/amici/returndata_matlab.h
+++ b/include/amici/returndata_matlab.h
@@ -55,7 +55,7 @@ void writeMatlabField0(
 template <typename T>
 void writeMatlabField1(
     mxArray* matlabStruct, char const* fieldName,
-    gsl::span<T const> const& fieldData, const mwSize dim0
+    gsl::span<T const> const& fieldData, mwSize const dim0
 );
 
 /**
diff --git a/include/amici/serialization.h b/include/amici/serialization.h
index c2a69d4b3a..6cab4f0353 100644
--- a/include/amici/serialization.h
+++ b/include/amici/serialization.h
@@ -34,11 +34,11 @@ void archiveVector(Archive& ar, T** p, int size) {
     if (Archive::is_loading::value) {
         if (*p != nullptr)
             delete[] *p;
-        ar& size;
+        ar & size;
         *p = size ? new T[size] : nullptr;
     } else {
         size = *p == nullptr ? 0 : size;
-        ar& size;
+        ar & size;
     }
     ar& make_array<T>(*p, size);
 }
@@ -51,40 +51,40 @@ void archiveVector(Archive& ar, T** p, int size) {
  */
 template <class Archive>
 void serialize(Archive& ar, amici::Solver& s, unsigned int const /*version*/) {
-    ar& s.sensi_;
-    ar& s.atol_;
-    ar& s.rtol_;
-    ar& s.atolB_;
-    ar& s.rtolB_;
-    ar& s.atol_fsa_;
-    ar& s.rtol_fsa_;
-    ar& s.quad_atol_;
-    ar& s.quad_rtol_;
-    ar& s.ss_tol_factor_;
-    ar& s.ss_atol_;
-    ar& s.ss_rtol_;
-    ar& s.ss_tol_sensi_factor_;
-    ar& s.ss_atol_sensi_;
-    ar& s.ss_rtol_sensi_;
-    ar& s.maxsteps_;
-    ar& s.maxstepsB_;
-    ar& s.newton_maxsteps_;
-    ar& s.newton_damping_factor_mode_;
-    ar& s.newton_damping_factor_lower_bound_;
-    ar& s.ism_;
-    ar& s.sensi_meth_;
-    ar& s.linsol_;
-    ar& s.interp_type_;
-    ar& s.lmm_;
-    ar& s.iter_;
-    ar& s.stldet_;
-    ar& s.ordering_;
-    ar& s.cpu_time_;
-    ar& s.cpu_timeB_;
-    ar& s.newton_step_steadystate_conv_;
-    ar& s.check_sensi_steadystate_conv_;
-    ar& s.rdata_mode_;
-    ar& s.maxtime_;
+    ar & s.sensi_;
+    ar & s.atol_;
+    ar & s.rtol_;
+    ar & s.atolB_;
+    ar & s.rtolB_;
+    ar & s.atol_fsa_;
+    ar & s.rtol_fsa_;
+    ar & s.quad_atol_;
+    ar & s.quad_rtol_;
+    ar & s.ss_tol_factor_;
+    ar & s.ss_atol_;
+    ar & s.ss_rtol_;
+    ar & s.ss_tol_sensi_factor_;
+    ar & s.ss_atol_sensi_;
+    ar & s.ss_rtol_sensi_;
+    ar & s.maxsteps_;
+    ar & s.maxstepsB_;
+    ar & s.newton_maxsteps_;
+    ar & s.newton_damping_factor_mode_;
+    ar & s.newton_damping_factor_lower_bound_;
+    ar & s.ism_;
+    ar & s.sensi_meth_;
+    ar & s.linsol_;
+    ar & s.interp_type_;
+    ar & s.lmm_;
+    ar & s.iter_;
+    ar & s.stldet_;
+    ar & s.ordering_;
+    ar & s.cpu_time_;
+    ar & s.cpu_timeB_;
+    ar & s.newton_step_steadystate_conv_;
+    ar & s.check_sensi_steadystate_conv_;
+    ar & s.rdata_mode_;
+    ar & s.maxtime_;
 }
 
 /**
@@ -99,11 +99,11 @@ void serialize(
 ) {
     Period tmp_period;
     if (Archive::is_loading::value) {
-        ar& tmp_period;
+        ar & tmp_period;
         d = std::chrono::duration<Period, Rep>(tmp_period);
     } else {
         tmp_period = d.count();
-        ar& tmp_period;
+        ar & tmp_period;
     }
 }
 
@@ -127,24 +127,24 @@ void serialize(
 template <class Archive>
 void serialize(Archive& ar, amici::Model& m, unsigned int const /*version*/) {
     ar& dynamic_cast<amici::ModelDimensions&>(m);
-    ar& m.simulation_parameters_;
-    ar& m.o2mode;
-    ar& m.z2event_;
-    ar& m.idlist;
-    ar& m.state_.h;
-    ar& m.state_.unscaledParameters;
-    ar& m.state_.fixedParameters;
-    ar& m.state_.plist;
-    ar& m.x0data_;
-    ar& m.sx0data_;
-    ar& m.nmaxevent_;
-    ar& m.state_is_non_negative_;
-    ar& m.pythonGenerated;
-    ar& m.min_sigma_;
-    ar& m.sigma_res_;
-    ar& m.steadystate_computation_mode_;
-    ar& m.steadystate_sensitivity_mode_;
-    ar& m.state_independent_events_;
+    ar & m.simulation_parameters_;
+    ar & m.o2mode;
+    ar & m.z2event_;
+    ar & m.idlist;
+    ar & m.state_.h;
+    ar & m.state_.unscaledParameters;
+    ar & m.state_.fixedParameters;
+    ar & m.state_.plist;
+    ar & m.x0data_;
+    ar & m.sx0data_;
+    ar & m.nmaxevent_;
+    ar & m.state_is_non_negative_;
+    ar & m.pythonGenerated;
+    ar & m.min_sigma_;
+    ar & m.sigma_res_;
+    ar & m.steadystate_computation_mode_;
+    ar & m.steadystate_sensitivity_mode_;
+    ar & m.state_independent_events_;
 }
 
 /**
@@ -156,18 +156,18 @@ template <class Archive>
 void serialize(
     Archive& ar, amici::SimulationParameters& s, unsigned int const /*version*/
 ) {
-    ar& s.fixedParameters;
-    ar& s.fixedParametersPreequilibration;
-    ar& s.fixedParametersPresimulation;
-    ar& s.parameters;
-    ar& s.x0;
-    ar& s.sx0;
-    ar& s.pscale;
-    ar& s.plist;
-    ar& s.ts_;
-    ar& s.tstart_;
-    ar& s.t_presim;
-    ar& s.reinitializeFixedParameterInitialStates;
+    ar & s.fixedParameters;
+    ar & s.fixedParametersPreequilibration;
+    ar & s.fixedParametersPresimulation;
+    ar & s.parameters;
+    ar & s.x0;
+    ar & s.sx0;
+    ar & s.pscale;
+    ar & s.plist;
+    ar & s.ts_;
+    ar & s.tstart_;
+    ar & s.t_presim;
+    ar & s.reinitializeFixedParameterInitialStates;
 }
 
 /**
@@ -181,63 +181,63 @@ void serialize(
     Archive& ar, amici::ReturnData& r, unsigned int const /*version*/
 ) {
     ar& dynamic_cast<amici::ModelDimensions&>(r);
-    ar& r.id;
-    ar& r.nx;
-    ar& r.nxtrue;
-    ar& r.nplist;
-    ar& r.nmaxevent;
-    ar& r.nt;
-    ar& r.newton_maxsteps;
-    ar& r.pscale;
-    ar& r.o2mode;
-    ar& r.sensi;
-    ar& r.sensi_meth;
-
-    ar& r.ts;
-    ar& r.xdot;
-    ar& r.J;
-    ar& r.w;
-    ar& r.z& r.sigmaz;
-    ar& r.sz& r.ssigmaz;
-    ar& r.rz;
-    ar& r.srz;
-    ar& r.s2rz;
-    ar& r.x;
-    ar& r.sx;
-    ar& r.y& r.sigmay;
-    ar& r.sy& r.ssigmay;
-
-    ar& r.numsteps;
-    ar& r.numstepsB;
-    ar& r.numrhsevals;
-    ar& r.numrhsevalsB;
-    ar& r.numerrtestfails;
-    ar& r.numerrtestfailsB;
-    ar& r.numnonlinsolvconvfails;
-    ar& r.numnonlinsolvconvfailsB;
-    ar& r.order;
-    ar& r.cpu_time;
-    ar& r.cpu_timeB;
-    ar& r.cpu_time_total;
-    ar& r.preeq_cpu_time;
-    ar& r.preeq_cpu_timeB;
-    ar& r.preeq_status;
-    ar& r.preeq_numsteps;
-    ar& r.preeq_wrms;
-    ar& r.preeq_t;
-    ar& r.posteq_cpu_time;
-    ar& r.posteq_cpu_timeB;
-    ar& r.posteq_status;
-    ar& r.posteq_numsteps;
-    ar& r.posteq_wrms;
-    ar& r.posteq_t;
-    ar& r.x0;
-    ar& r.sx0;
-    ar& r.llh;
-    ar& r.chi2;
-    ar& r.sllh;
-    ar& r.s2llh;
-    ar& r.status;
+    ar & r.id;
+    ar & r.nx;
+    ar & r.nxtrue;
+    ar & r.nplist;
+    ar & r.nmaxevent;
+    ar & r.nt;
+    ar & r.newton_maxsteps;
+    ar & r.pscale;
+    ar & r.o2mode;
+    ar & r.sensi;
+    ar & r.sensi_meth;
+
+    ar & r.ts;
+    ar & r.xdot;
+    ar & r.J;
+    ar & r.w;
+    ar & r.z & r.sigmaz;
+    ar & r.sz & r.ssigmaz;
+    ar & r.rz;
+    ar & r.srz;
+    ar & r.s2rz;
+    ar & r.x;
+    ar & r.sx;
+    ar & r.y & r.sigmay;
+    ar & r.sy & r.ssigmay;
+
+    ar & r.numsteps;
+    ar & r.numstepsB;
+    ar & r.numrhsevals;
+    ar & r.numrhsevalsB;
+    ar & r.numerrtestfails;
+    ar & r.numerrtestfailsB;
+    ar & r.numnonlinsolvconvfails;
+    ar & r.numnonlinsolvconvfailsB;
+    ar & r.order;
+    ar & r.cpu_time;
+    ar & r.cpu_timeB;
+    ar & r.cpu_time_total;
+    ar & r.preeq_cpu_time;
+    ar & r.preeq_cpu_timeB;
+    ar & r.preeq_status;
+    ar & r.preeq_numsteps;
+    ar & r.preeq_wrms;
+    ar & r.preeq_t;
+    ar & r.posteq_cpu_time;
+    ar & r.posteq_cpu_timeB;
+    ar & r.posteq_status;
+    ar & r.posteq_numsteps;
+    ar & r.posteq_wrms;
+    ar & r.posteq_t;
+    ar & r.x0;
+    ar & r.sx0;
+    ar & r.llh;
+    ar & r.chi2;
+    ar & r.sllh;
+    ar & r.s2llh;
+    ar & r.status;
 }
 
 /**
@@ -250,30 +250,30 @@ template <class Archive>
 void serialize(
     Archive& ar, amici::ModelDimensions& m, unsigned int const /*version*/
 ) {
-    ar& m.nx_rdata;
-    ar& m.nxtrue_rdata;
-    ar& m.nx_solver;
-    ar& m.nxtrue_solver;
-    ar& m.nx_solver_reinit;
-    ar& m.np;
-    ar& m.nk;
-    ar& m.ny;
-    ar& m.nytrue;
-    ar& m.nz;
-    ar& m.nztrue;
-    ar& m.ne;
-    ar& m.ne_solver;
-    ar& m.nspl;
-    ar& m.nw;
-    ar& m.ndwdx;
-    ar& m.ndwdp;
-    ar& m.ndwdw;
-    ar& m.ndxdotdw;
-    ar& m.ndJydy;
-    ar& m.nnz;
-    ar& m.nJ;
-    ar& m.ubw;
-    ar& m.lbw;
+    ar & m.nx_rdata;
+    ar & m.nxtrue_rdata;
+    ar & m.nx_solver;
+    ar & m.nxtrue_solver;
+    ar & m.nx_solver_reinit;
+    ar & m.np;
+    ar & m.nk;
+    ar & m.ny;
+    ar & m.nytrue;
+    ar & m.nz;
+    ar & m.nztrue;
+    ar & m.ne;
+    ar & m.ne_solver;
+    ar & m.nspl;
+    ar & m.nw;
+    ar & m.ndwdx;
+    ar & m.ndwdp;
+    ar & m.ndwdw;
+    ar & m.ndxdotdw;
+    ar & m.ndJydy;
+    ar & m.nnz;
+    ar & m.nJ;
+    ar & m.ubw;
+    ar & m.lbw;
 }
 #endif
 } // namespace serialization
diff --git a/include/amici/simulation_parameters.h b/include/amici/simulation_parameters.h
index 8a1c1ec23d..110b8fc03e 100644
--- a/include/amici/simulation_parameters.h
+++ b/include/amici/simulation_parameters.h
@@ -37,8 +37,16 @@ class SimulationParameters {
 
 #ifndef SWIGPYTHON
     /*
-     * include/amici/simulation_parameters.h:71: Warning 509: Overloaded method amici::SimulationParameters::SimulationParameters(std::vector< amici::realtype,std::allocator< amici::realtype > >,std::vector< amici::realtype,std::allocator< amici::realtype > >,std::vector< amici::realtype,std::allocator< amici::realtype > >) effectively ignored,
-     * include/amici/simulation_parameters.h:54: Warning 509: as it is shadowed by amici::SimulationParameters::SimulationParameters(std::vector< amici::realtype,std::allocator< amici::realtype > >,std::vector< amici::realtype,std::allocator< amici::realtype > >,std::vector< int,std::allocator< int > >).
+     * include/amici/simulation_parameters.h:71: Warning 509: Overloaded method
+     * amici::SimulationParameters::SimulationParameters(std::vector<
+     * amici::realtype,std::allocator< amici::realtype > >,std::vector<
+     * amici::realtype,std::allocator< amici::realtype > >,std::vector<
+     * amici::realtype,std::allocator< amici::realtype > >) effectively ignored,
+     * include/amici/simulation_parameters.h:54: Warning 509: as it is shadowed
+     * by amici::SimulationParameters::SimulationParameters(std::vector<
+     * amici::realtype,std::allocator< amici::realtype > >,std::vector<
+     * amici::realtype,std::allocator< amici::realtype > >,std::vector<
+     * int,std::allocator< int > >).
      */
     /**
      * @brief Constructor
diff --git a/include/amici/solver.h b/include/amici/solver.h
index 4a1c95b96a..92ae4c47fb 100644
--- a/include/amici/solver.h
+++ b/include/amici/solver.h
@@ -134,7 +134,7 @@ class Solver {
      */
 
     void setupSteadystate(
-        const realtype t0, Model* model, AmiVector const& x0,
+        realtype const t0, Model* model, AmiVector const& x0,
         AmiVector const& dx0, AmiVector const& xB0, AmiVector const& dxB0,
         AmiVector const& xQ0
     ) const;
@@ -1612,7 +1612,7 @@ class Solver {
     mutable std::unique_ptr<void, free_solver_ptr> solver_memory_;
 
     /** pointer to solver memory block */
-    mutable std::vector<std::unique_ptr<void,free_solver_ptr>>
+    mutable std::vector<std::unique_ptr<void, free_solver_ptr>>
         solver_memory_B_;
 
     /** Sundials user_data */
@@ -1709,7 +1709,7 @@ class Solver {
      * @param preequilibration flag indicating preequilibration or simulation
      */
     void checkSensitivityMethod(
-        const SensitivityMethod sensi_meth, bool preequilibration
+        SensitivityMethod const sensi_meth, bool preequilibration
     ) const;
 
     /** state (dimension: nx_solver) */
@@ -1784,7 +1784,8 @@ class Solver {
 
     /** Damping factor state used int the Newton method */
     NewtonDampingFactorMode newton_damping_factor_mode_{
-        NewtonDampingFactorMode::on};
+        NewtonDampingFactorMode::on
+    };
 
     /** Lower bound of the damping factor. */
     realtype newton_damping_factor_lower_bound_{1e-8};
diff --git a/include/amici/solver_cvodes.h b/include/amici/solver_cvodes.h
index d6d1dcea24..f98fb34c1a 100644
--- a/include/amici/solver_cvodes.h
+++ b/include/amici/solver_cvodes.h
@@ -218,7 +218,7 @@ class CVodeSolver : public Solver {
     init(realtype t0, AmiVector const& x0, AmiVector const& dx0) const override;
 
     void initSteadystate(
-        const realtype t0, AmiVector const& x0, AmiVector const& dx0
+        realtype const t0, AmiVector const& x0, AmiVector const& dx0
     ) const override;
 
     void sensInit1(AmiVectorArray const& sx0, AmiVectorArray const& sdx0)
diff --git a/include/amici/solver_idas.h b/include/amici/solver_idas.h
index 0dba1a9504..b985a090af 100644
--- a/include/amici/solver_idas.h
+++ b/include/amici/solver_idas.h
@@ -198,7 +198,7 @@ class IDASolver : public Solver {
     init(realtype t0, AmiVector const& x0, AmiVector const& dx0) const override;
 
     void initSteadystate(
-        const realtype t0, AmiVector const& x0, AmiVector const& dx0
+        realtype const t0, AmiVector const& x0, AmiVector const& dx0
     ) const override;
 
     void sensInit1(AmiVectorArray const& sx0, AmiVectorArray const& sdx0)
diff --git a/include/amici/splinefunctions.h b/include/amici/splinefunctions.h
index db4410de91..ed9afb9ffc 100644
--- a/include/amici/splinefunctions.h
+++ b/include/amici/splinefunctions.h
@@ -72,7 +72,7 @@ class AbstractSpline {
      * @param t point at which the spline is to be evaluated
      * @return value of the spline at `t`
      */
-    realtype get_value(const realtype t) const;
+    realtype get_value(realtype const t) const;
 
     /**
      * @brief Get the value of this spline at a given point
@@ -80,7 +80,7 @@ class AbstractSpline {
      * @param t point at which the spline is to be evaluated
      * @return scaled value of the spline at `t`
      */
-    virtual realtype get_value_scaled(const realtype t) const = 0;
+    virtual realtype get_value_scaled(realtype const t) const = 0;
 
     /**
      * @brief Get the value of this spline at a given node
@@ -105,7 +105,7 @@ class AbstractSpline {
      * @return sensitivity of the spline with respect to the `ip`th parameter
      * at `t`
      */
-    realtype get_sensitivity(const realtype t, int const ip) const;
+    realtype get_sensitivity(realtype const t, int const ip) const;
 
     /**
      * @brief Get the derivative of this spline with respect to a given
@@ -119,7 +119,7 @@ class AbstractSpline {
      * at `t`
      */
     realtype
-    get_sensitivity(const realtype t, int const ip, const realtype value) const;
+    get_sensitivity(realtype const t, int const ip, realtype const value) const;
 
     /**
      * @brief Get the derivative of this spline with respect to a given
@@ -131,7 +131,7 @@ class AbstractSpline {
      * parameter at `t`
      */
     virtual realtype
-    get_sensitivity_scaled(const realtype t, int const ip) const
+    get_sensitivity_scaled(realtype const t, int const ip) const
         = 0;
 
     /**
@@ -329,7 +329,7 @@ class HermiteSpline : public AbstractSpline {
         gsl::span<realtype> dspline_slopesdp
     ) override;
 
-    realtype get_value_scaled(const realtype t) const override;
+    realtype get_value_scaled(realtype const t) const override;
 
     /**
      * @brief Get the derivative of the spline at a given node
@@ -347,7 +347,7 @@ class HermiteSpline : public AbstractSpline {
     realtype get_node_derivative_scaled(int const i) const;
 
     realtype
-    get_sensitivity_scaled(const realtype t, int const ip) const override;
+    get_sensitivity_scaled(realtype const t, int const ip) const override;
 
     /**
      * @brief Whether derivatives of this spline are computed
diff --git a/include/amici/sundials_matrix_wrapper.h b/include/amici/sundials_matrix_wrapper.h
index 8d63eca5ea..ee2516f78d 100644
--- a/include/amici/sundials_matrix_wrapper.h
+++ b/include/amici/sundials_matrix_wrapper.h
@@ -266,7 +266,7 @@ class SUNMatrixWrapper {
      * @brief Set the index values of a sparse matrix
      * @param vals rows (CSC) or columns (CSR) for data entries
      */
-    void set_indexvals(const gsl::span<sunindextype const> vals) {
+    void set_indexvals(gsl::span<sunindextype const> const vals) {
         assert(matrix_);
         assert(matrix_id() == SUNMATRIX_SPARSE);
         assert(gsl::narrow<sunindextype>(vals.size()) == capacity());
@@ -309,7 +309,7 @@ class SUNMatrixWrapper {
      * @param ptrs starting data-indices where the columns (CSC) or rows (CSR)
      * start
      */
-    void set_indexptrs(const gsl::span<sunindextype const> ptrs) {
+    void set_indexptrs(gsl::span<sunindextype const> const ptrs) {
         assert(matrix_);
         assert(matrix_id() == SUNMATRIX_SPARSE);
         assert(gsl::narrow<sunindextype>(ptrs.size()) == num_indexptrs() + 1);
@@ -357,7 +357,7 @@ class SUNMatrixWrapper {
      */
     void multiply(
         gsl::span<realtype> c, gsl::span<realtype const> b,
-        const realtype alpha = 1.0
+        realtype const alpha = 1.0
     ) const;
 
     /**
@@ -440,8 +440,8 @@ class SUNMatrixWrapper {
      * @return updated number of nonzeros in C
      */
     sunindextype scatter(
-        const sunindextype k, const realtype beta, sunindextype* w,
-        gsl::span<realtype> x, const sunindextype mark, SUNMatrixWrapper* C,
+        sunindextype const k, realtype const beta, sunindextype* w,
+        gsl::span<realtype> x, sunindextype const mark, SUNMatrixWrapper* C,
         sunindextype nnz
     ) const;
 
@@ -455,7 +455,7 @@ class SUNMatrixWrapper {
      * set to ncols/nrows
      */
     void transpose(
-        SUNMatrixWrapper& C, const realtype alpha, sunindextype blocksize
+        SUNMatrixWrapper& C, realtype const alpha, sunindextype blocksize
     ) const;
 
     /**
diff --git a/python/sdist/amici/__init__.py b/python/sdist/amici/__init__.py
index b06ea6d0f3..204933565f 100644
--- a/python/sdist/amici/__init__.py
+++ b/python/sdist/amici/__init__.py
@@ -13,7 +13,7 @@
 import sys
 from pathlib import Path
 from types import ModuleType as ModelModule
-from typing import Any, Callable, Optional, Union
+from typing import Any, Callable, Union
 
 
 def _get_amici_path():
@@ -153,7 +153,7 @@ def __exit__(self, exc_type, exc_value, traceback):
 
 
 def import_model_module(
-    module_name: str, module_path: Optional[Union[Path, str]] = None
+    module_name: str, module_path: Union[Path, str]
 ) -> ModelModule:
     """
     Import Python module of an AMICI model
diff --git a/python/sdist/amici/antimony_import.py b/python/sdist/amici/antimony_import.py
index cc23079787..545a2654bd 100644
--- a/python/sdist/amici/antimony_import.py
+++ b/python/sdist/amici/antimony_import.py
@@ -28,11 +28,13 @@ def antimony2sbml(ant_model: Union[str, Path]) -> str:
         is_file = False
 
     if is_file:
-        status = ant.loadAntimonyFile(ant_model)
+        status = ant.loadAntimonyFile(str(ant_model))
     else:
         status = ant.loadAntimonyString(ant_model)
     if status < 0:
-        raise RuntimeError("Antimony model could not be loaded.")
+        raise RuntimeError(
+            f"Antimony model could not be loaded: {ant.getLastError()}"
+        )
 
     if (main_module_name := ant.getMainModuleName()) is None:
         raise AssertionError("There is no Antimony module.")
diff --git a/python/sdist/amici/swig_wrappers.py b/python/sdist/amici/swig_wrappers.py
index e16d4c2580..e942cc76dd 100644
--- a/python/sdist/amici/swig_wrappers.py
+++ b/python/sdist/amici/swig_wrappers.py
@@ -136,17 +136,22 @@ def ExpData(*args) -> "amici_swig.ExpData":
 
     :returns: ExpData Instance
     """
+    if not args:
+        return amici_swig.ExpData()
+
     if isinstance(args[0], numpy.ReturnDataView):
         return amici_swig.ExpData(_get_ptr(args[0]["ptr"]), *args[1:])
-    elif isinstance(args[0], (amici_swig.ExpData, amici_swig.ExpDataPtr)):
+
+    if isinstance(args[0], (amici_swig.ExpData, amici_swig.ExpDataPtr)):
         # the *args[:1] should be empty, but by the time you read this,
         # the constructor signature may have changed, and you are glad this
         # wrapper did not break.
         return amici_swig.ExpData(_get_ptr(args[0]), *args[1:])
-    elif isinstance(args[0], (amici_swig.Model, amici_swig.ModelPtr)):
+
+    if isinstance(args[0], (amici_swig.Model, amici_swig.ModelPtr)):
         return amici_swig.ExpData(_get_ptr(args[0]))
-    else:
-        return amici_swig.ExpData(*args)
+
+    return amici_swig.ExpData(*args)
 
 
 def runAmiciSimulations(
diff --git a/python/sdist/setup.cfg b/python/sdist/setup.cfg
index cc82ef08de..009d622a6b 100644
--- a/python/sdist/setup.cfg
+++ b/python/sdist/setup.cfg
@@ -35,6 +35,7 @@ install_requires =
     numpy; python_version>='3.12'
     python-libsbml
     pandas>=2.0.2
+    pyarrow
     wurlitzer
     toposort
     setuptools>=48
diff --git a/python/tests/valgrind-python.supp b/python/tests/valgrind-python.supp
index 9cd3f5de3b..b4e0177033 100644
--- a/python/tests/valgrind-python.supp
+++ b/python/tests/valgrind-python.supp
@@ -547,6 +547,14 @@
    ...
 }
 
+{
+   Pandas 2.2.0
+   Memcheck:Leak
+   match-leak-kinds: definite
+   fun:realloc
+   obj:*site-packages/pandas/*
+}
+
 {
    Scipy extensions
    Memcheck:Leak
@@ -728,6 +736,16 @@
    fun:_PyImport_FindSharedFuncptr
 }
 
+{
+   _dl_catch_exception
+   Memcheck:Addr8
+   fun:strncmp
+   fun:is_dst
+   ...
+   fun:_dl_catch_exception
+   ...
+}
+
 {
    Python dictkeys_get_index
    Memcheck:Value8
diff --git a/scripts/downloadAndBuildSwig.sh b/scripts/downloadAndBuildSwig.sh
index 5fa0896f6c..29ce32890b 100755
--- a/scripts/downloadAndBuildSwig.sh
+++ b/scripts/downloadAndBuildSwig.sh
@@ -7,7 +7,7 @@ set -euo pipefail
 SCRIPT_PATH=$(dirname "$BASH_SOURCE")
 AMICI_PATH=$(cd "$SCRIPT_PATH/.." && pwd)
 
-swig_version="${1:-"4.1.1"}"
+swig_version="${1:-"4.2.0"}"
 SWIG_ARCHIVE="swig-${swig_version}.tar.gz"
 SWIG_URL="http://downloads.sourceforge.net/project/swig/swig/swig-${swig_version}/${SWIG_ARCHIVE}"
 SWIG_DIR="swig-${swig_version}"
@@ -18,13 +18,13 @@ cd "${AMICI_PATH}/ThirdParty/"
 
 if [[ ! -d ${SWIG_DIR} ]]; then
     if [[ ! -f ${SWIG_ARCHIVE} ]]
-        then wget ${SWIG_URL}
+        then wget "${SWIG_URL}"
     fi
 
-    tar -xzf ${SWIG_ARCHIVE}
+    tar -xzf "${SWIG_ARCHIVE}"
 fi
 
-cd ${SWIG_DIR}
+cd "${SWIG_DIR}"
 ./configure \
   --prefix="${PREFIX}" \
   --without-alllang \
@@ -38,6 +38,6 @@ echo "================"
 echo "SWIG installation successful"
 echo
 echo "To use this version of SWIG, add directory ${SWIG_BIN_DIR} to your PATH,"
-echo "e.g. adding the following line to your .bashrc:"
+echo "e.g., adding the following line to your ~/.bashrc:"
 echo "    export PATH=${SWIG_BIN_DIR}:\$PATH"
 echo "================"
diff --git a/src/abstract_model.cpp b/src/abstract_model.cpp
index dc2c469173..041d06b361 100644
--- a/src/abstract_model.cpp
+++ b/src/abstract_model.cpp
@@ -11,7 +11,7 @@ std::string AbstractModel::getAmiciCommit() const {
 }
 
 void AbstractModel::
-    fx0(realtype* /*x0*/, const realtype /*t*/, realtype const* /*p*/,
+    fx0(realtype* /*x0*/, realtype const /*t*/, realtype const* /*p*/,
         realtype const* /*k*/) {
     throw AmiException(
         "Requested functionality is not supported as %s is "
@@ -25,14 +25,14 @@ bool AbstractModel::isFixedParameterStateReinitializationAllowed() const {
 }
 
 void AbstractModel::fx0_fixedParameters(
-    realtype* /*x0*/, const realtype /*t*/, realtype const* /*p*/,
+    realtype* /*x0*/, realtype const /*t*/, realtype const* /*p*/,
     realtype const* /*k*/, gsl::span<int const> /*reinitialization_state_idxs*/
 ) {
     // no-op default implementation
 }
 
 void AbstractModel::fsx0_fixedParameters(
-    realtype* /*sx0*/, const realtype /*t*/, realtype const* /*x0*/,
+    realtype* /*sx0*/, realtype const /*t*/, realtype const* /*x0*/,
     realtype const* /*p*/, realtype const* /*k*/, int const /*ip*/,
     gsl::span<int const> /*reinitialization_state_idxs*/
 ) {
@@ -40,7 +40,7 @@ void AbstractModel::fsx0_fixedParameters(
 }
 
 void AbstractModel::fsx0(
-    realtype* /*sx0*/, const realtype /*t*/, realtype const* /*x0*/,
+    realtype* /*sx0*/, realtype const /*t*/, realtype const* /*x0*/,
     realtype const* /*p*/, realtype const* /*k*/, int const /*ip*/
 ) {
     throw AmiException(
@@ -55,7 +55,7 @@ void AbstractModel::fdx0(AmiVector& /*x0*/, AmiVector& /*dx0*/) {
 }
 
 void AbstractModel::fstau(
-    realtype* /*stau*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*stau*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*tcl*/, realtype const* /*sx*/, int const /*ip*/,
     int const /*ie*/
@@ -68,7 +68,7 @@ void AbstractModel::fstau(
 }
 
 void AbstractModel::
-    fy(realtype* /*y*/, const realtype /*t*/, realtype const* /*x*/,
+    fy(realtype* /*y*/, realtype const /*t*/, realtype const* /*x*/,
        realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
        realtype const* /*w*/) {
     throw AmiException(
@@ -79,7 +79,7 @@ void AbstractModel::
 }
 
 void AbstractModel::fdydp(
-    realtype* /*dydp*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dydp*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     int const /*ip*/, realtype const* /*w*/, realtype const* /*dwdp*/
 ) {
@@ -91,7 +91,7 @@ void AbstractModel::fdydp(
 }
 
 void AbstractModel::fdydp(
-    realtype* /*dydp*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dydp*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     int /*ip*/, realtype const* /*w*/, realtype const* /*tcl*/,
     realtype const* /*dtcldp*/, realtype const* /*spl*/,
@@ -105,7 +105,7 @@ void AbstractModel::fdydp(
 }
 
 void AbstractModel::fdydx(
-    realtype* /*dydx*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dydx*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/, realtype const* /*dwdx*/
 ) {
@@ -117,7 +117,7 @@ void AbstractModel::fdydx(
 }
 
 void AbstractModel::
-    fz(realtype* /*z*/, int const /*ie*/, const realtype /*t*/,
+    fz(realtype* /*z*/, int const /*ie*/, realtype const /*t*/,
        realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
        realtype const* /*h*/) {
     throw AmiException(
@@ -128,7 +128,7 @@ void AbstractModel::
 }
 
 void AbstractModel::
-    fsz(realtype* /*sz*/, int const /*ie*/, const realtype /*t*/,
+    fsz(realtype* /*sz*/, int const /*ie*/, realtype const /*t*/,
         realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
         realtype const* /*h*/, realtype const* /*sx*/, int const /*ip*/) {
     throw AmiException(
@@ -139,7 +139,7 @@ void AbstractModel::
 }
 
 void AbstractModel::
-    frz(realtype* /*rz*/, int const /*ie*/, const realtype /*t*/,
+    frz(realtype* /*rz*/, int const /*ie*/, realtype const /*t*/,
         realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
         realtype const* /*h*/) {
     throw AmiException(
@@ -150,7 +150,7 @@ void AbstractModel::
 }
 
 void AbstractModel::fsrz(
-    realtype* /*srz*/, int const /*ie*/, const realtype /*t*/,
+    realtype* /*srz*/, int const /*ie*/, realtype const /*t*/,
     realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
     realtype const* /*h*/, realtype const* /*sx*/, int const /*ip*/
 ) {
@@ -162,7 +162,7 @@ void AbstractModel::fsrz(
 }
 
 void AbstractModel::fdzdp(
-    realtype* /*dzdp*/, int const /*ie*/, const realtype /*t*/,
+    realtype* /*dzdp*/, int const /*ie*/, realtype const /*t*/,
     realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
     realtype const* /*h*/, int const /*ip*/
 ) {
@@ -174,7 +174,7 @@ void AbstractModel::fdzdp(
 }
 
 void AbstractModel::fdzdx(
-    realtype* /*dzdx*/, int const /*ie*/, const realtype /*t*/,
+    realtype* /*dzdx*/, int const /*ie*/, realtype const /*t*/,
     realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
     realtype const* /*h*/
 ) {
@@ -186,7 +186,7 @@ void AbstractModel::fdzdx(
 }
 
 void AbstractModel::fdrzdp(
-    realtype* /*drzdp*/, int const /*ie*/, const realtype /*t*/,
+    realtype* /*drzdp*/, int const /*ie*/, realtype const /*t*/,
     realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
     realtype const* /*h*/, int const /*ip*/
 ) {
@@ -198,7 +198,7 @@ void AbstractModel::fdrzdp(
 }
 
 void AbstractModel::fdrzdx(
-    realtype* /*drzdx*/, int const /*ie*/, const realtype /*t*/,
+    realtype* /*drzdx*/, int const /*ie*/, realtype const /*t*/,
     realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
     realtype const* /*h*/
 ) {
@@ -210,7 +210,7 @@ void AbstractModel::fdrzdx(
 }
 
 void AbstractModel::fdeltax(
-    realtype* /*deltax*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*deltax*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     int const /*ie*/, realtype const* /*xdot*/, realtype const* /*xdot_old*/
 ) {
@@ -222,7 +222,7 @@ void AbstractModel::fdeltax(
 }
 
 void AbstractModel::fdeltasx(
-    realtype* /*deltasx*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*deltasx*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/, int const /*ip*/, int const /*ie*/,
     realtype const* /*xdot*/, realtype const* /*xdot_old*/,
@@ -236,7 +236,7 @@ void AbstractModel::fdeltasx(
 }
 
 void AbstractModel::fdeltaxB(
-    realtype* /*deltaxB*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*deltaxB*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     int const /*ie*/, realtype const* /*xdot*/, realtype const* /*xdot_old*/,
     realtype const* /*xB*/
@@ -249,7 +249,7 @@ void AbstractModel::fdeltaxB(
 }
 
 void AbstractModel::fdeltaqB(
-    realtype* /*deltaqB*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*deltaqB*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     int const /*ip*/, int const /*ie*/, realtype const* /*xdot*/,
     realtype const* /*xdot_old*/, realtype const* /*xB*/
@@ -262,7 +262,7 @@ void AbstractModel::fdeltaqB(
 }
 
 void AbstractModel::fsigmay(
-    realtype* /*sigmay*/, const realtype /*t*/, realtype const* /*p*/,
+    realtype* /*sigmay*/, realtype const /*t*/, realtype const* /*p*/,
     realtype const* /*k*/, realtype const* /*y*/
 ) {
     throw AmiException(
@@ -273,7 +273,7 @@ void AbstractModel::fsigmay(
 }
 
 void AbstractModel::fdsigmaydp(
-    realtype* /*dsigmaydp*/, const realtype /*t*/, realtype const* /*p*/,
+    realtype* /*dsigmaydp*/, realtype const /*t*/, realtype const* /*p*/,
     realtype const* /*k*/, realtype const* /*y*/, int const /*ip*/
 ) {
     throw AmiException(
@@ -284,7 +284,7 @@ void AbstractModel::fdsigmaydp(
 }
 
 void AbstractModel::fdsigmaydy(
-    realtype* /*dsigmaydy*/, const realtype /*t*/, realtype const* /*p*/,
+    realtype* /*dsigmaydy*/, realtype const /*t*/, realtype const* /*p*/,
     realtype const* /*k*/, realtype const* /*y*/
 ) {
     throw AmiException(
@@ -295,7 +295,7 @@ void AbstractModel::fdsigmaydy(
 }
 
 void AbstractModel::fsigmaz(
-    realtype* /*sigmaz*/, const realtype /*t*/, realtype const* /*p*/,
+    realtype* /*sigmaz*/, realtype const /*t*/, realtype const* /*p*/,
     realtype const* /*k*/
 ) {
     throw AmiException(
@@ -306,7 +306,7 @@ void AbstractModel::fsigmaz(
 }
 
 void AbstractModel::fdsigmazdp(
-    realtype* /*dsigmazdp*/, const realtype /*t*/, realtype const* /*p*/,
+    realtype* /*dsigmazdp*/, realtype const /*t*/, realtype const* /*p*/,
     realtype const* /*k*/, int const /*ip*/
 ) {
     throw AmiException(
@@ -438,7 +438,7 @@ void AbstractModel::fdJrzdsigma(
 }
 
 void AbstractModel::
-    fw(realtype* /*w*/, const realtype /*t*/, realtype const* /*x*/,
+    fw(realtype* /*w*/, realtype const /*t*/, realtype const* /*x*/,
        realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
        realtype const* /*tcl*/, realtype const* /*spl*/) {
     throw AmiException(
@@ -449,7 +449,7 @@ void AbstractModel::
 }
 
 void AbstractModel::fdwdp(
-    realtype* /*dwdp*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dwdp*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/, realtype const* /*tcl*/, realtype const* /*stcl*/,
     realtype const* /*spl*/, realtype const* /*sspl*/
@@ -478,7 +478,7 @@ void AbstractModel::fdwdp_rowvals(SUNMatrixWrapper& /*dwdp*/) {
 }
 
 void AbstractModel::fdwdp(
-    realtype* /*dwdp*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dwdp*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/, realtype const* /*tcl*/, realtype const* /*stcl*/,
     realtype const* /*spl*/, realtype const* /*sspl*/, int const /*ip*/
@@ -491,7 +491,7 @@ void AbstractModel::fdwdp(
 }
 
 void AbstractModel::fdwdx(
-    realtype* /*dwdx*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dwdx*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/, realtype const* /*tcl*/, realtype const* /*spl*/
 ) {
diff --git a/src/edata.cpp b/src/edata.cpp
index 05499bb9cf..2408e4becd 100644
--- a/src/edata.cpp
+++ b/src/edata.cpp
@@ -194,7 +194,7 @@ void ExpData::setObservedDataStdDev(
         observed_data_std_dev_.clear();
 }
 
-void ExpData::setObservedDataStdDev(const realtype stdDev) {
+void ExpData::setObservedDataStdDev(realtype const stdDev) {
     checkSigmaPositivity(stdDev, "stdDev");
     std::fill(
         observed_data_std_dev_.begin(), observed_data_std_dev_.end(), stdDev
@@ -216,7 +216,7 @@ void ExpData::setObservedDataStdDev(
             = observedDataStdDev.at(it);
 }
 
-void ExpData::setObservedDataStdDev(const realtype stdDev, int iy) {
+void ExpData::setObservedDataStdDev(realtype const stdDev, int iy) {
     checkSigmaPositivity(stdDev, "stdDev");
     for (int it = 0; it < nt(); ++it)
         observed_data_std_dev_.at(iy + it * nytrue_) = stdDev;
@@ -290,7 +290,7 @@ void ExpData::setObservedEventsStdDev(
         observed_events_std_dev_.clear();
 }
 
-void ExpData::setObservedEventsStdDev(const realtype stdDev) {
+void ExpData::setObservedEventsStdDev(realtype const stdDev) {
     checkSigmaPositivity(stdDev, "stdDev");
     std::fill(
         observed_events_std_dev_.begin(), observed_events_std_dev_.end(), stdDev
@@ -313,7 +313,7 @@ void ExpData::setObservedEventsStdDev(
             = observedEventsStdDev.at(ie);
 }
 
-void ExpData::setObservedEventsStdDev(const realtype stdDev, int iz) {
+void ExpData::setObservedEventsStdDev(realtype const stdDev, int iz) {
     checkSigmaPositivity(stdDev, "stdDev");
 
     for (int ie = 0; ie < nmaxevent_; ++ie)
@@ -392,7 +392,7 @@ void checkSigmaPositivity(
         checkSigmaPositivity(sigma, vectorName);
 }
 
-void checkSigmaPositivity(const realtype sigma, char const* sigmaName) {
+void checkSigmaPositivity(realtype const sigma, char const* sigmaName) {
     if (sigma <= 0.0)
         throw AmiException(
             "Encountered sigma <= 0 in %s! value: %f", sigmaName, sigma
diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp
index c090f51280..b4bcd9740f 100644
--- a/src/forwardproblem.cpp
+++ b/src/forwardproblem.cpp
@@ -22,11 +22,12 @@ namespace amici {
  */
 bool is_next_t_too_close(realtype cur_t, realtype t_next) {
     auto tdiff = t_next - cur_t;
-    if(tdiff == 0.0)
+    if (tdiff == 0.0)
         return true;
 
     auto tdist = std::fabs(tdiff);
-    auto tround = std::numeric_limits<realtype>::epsilon() * std::max(std::fabs(cur_t), std::fabs(t_next));
+    auto tround = std::numeric_limits<realtype>::epsilon()
+                  * std::max(std::fabs(cur_t), std::fabs(t_next));
     if (tdist < 2.0 * tround)
         return true;
 
diff --git a/src/hdf5.cpp b/src/hdf5.cpp
index d9d875cdf7..e459eee138 100644
--- a/src/hdf5.cpp
+++ b/src/hdf5.cpp
@@ -655,7 +655,7 @@ void createAndWriteDouble2DDataset(
     const H5::H5File& file, std::string const& datasetName,
     gsl::span<double const> buffer, hsize_t m, hsize_t n
 ) {
-    const hsize_t adims[]{m, n};
+    hsize_t const adims[]{m, n};
     H5::DataSpace dataspace(2, adims);
     auto dataset = file.createDataSet(
         datasetName.c_str(), H5::PredType::NATIVE_DOUBLE, dataspace
@@ -667,7 +667,7 @@ void createAndWriteInt2DDataset(
     H5::H5File const& file, std::string const& datasetName,
     gsl::span<int const> buffer, hsize_t m, hsize_t n
 ) {
-    const hsize_t adims[]{m, n};
+    hsize_t const adims[]{m, n};
     H5::DataSpace dataspace(2, adims);
     auto dataset = file.createDataSet(
         datasetName.c_str(), H5::PredType::NATIVE_INT, dataspace
@@ -679,7 +679,7 @@ void createAndWriteDouble3DDataset(
     H5::H5File const& file, std::string const& datasetName,
     gsl::span<double const> buffer, hsize_t m, hsize_t n, hsize_t o
 ) {
-    const hsize_t adims[]{m, n, o};
+    hsize_t const adims[]{m, n, o};
     H5::DataSpace dataspace(3, adims);
     auto dataset = file.createDataSet(
         datasetName.c_str(), H5::PredType::NATIVE_DOUBLE, dataspace
diff --git a/src/interface_matlab.cpp b/src/interface_matlab.cpp
index 3caae66a96..e0f4702f34 100644
--- a/src/interface_matlab.cpp
+++ b/src/interface_matlab.cpp
@@ -71,12 +71,12 @@ void amici_dgemm(
 ) {
     assert(layout == BLASLayout::colMajor);
 
-    const ptrdiff_t M_ = M;
-    const ptrdiff_t N_ = N;
-    const ptrdiff_t K_ = K;
-    const ptrdiff_t lda_ = lda;
-    const ptrdiff_t ldb_ = ldb;
-    const ptrdiff_t ldc_ = ldc;
+    ptrdiff_t const M_ = M;
+    ptrdiff_t const N_ = N;
+    ptrdiff_t const K_ = K;
+    ptrdiff_t const lda_ = lda;
+    ptrdiff_t const ldb_ = ldb;
+    ptrdiff_t const ldc_ = ldc;
     char const transA = amici_blasCBlasTransToBlasTrans(TransA);
     char const transB = amici_blasCBlasTransToBlasTrans(TransB);
 
@@ -92,11 +92,11 @@ void amici_dgemv(
 ) {
     assert(layout == BLASLayout::colMajor);
 
-    const ptrdiff_t M_ = M;
-    const ptrdiff_t N_ = N;
-    const ptrdiff_t lda_ = lda;
-    const ptrdiff_t incX_ = incX;
-    const ptrdiff_t incY_ = incY;
+    ptrdiff_t const M_ = M;
+    ptrdiff_t const N_ = N;
+    ptrdiff_t const lda_ = lda;
+    ptrdiff_t const incX_ = incX;
+    ptrdiff_t const incY_ = incY;
     char const transA = amici_blasCBlasTransToBlasTrans(TransA);
 
     FORTRAN_WRAPPER(dgemv)
@@ -107,9 +107,9 @@ void amici_daxpy(
     int n, double alpha, double const* x, int const incx, double* y, int incy
 ) {
 
-    const ptrdiff_t n_ = n;
-    const ptrdiff_t incx_ = incx;
-    const ptrdiff_t incy_ = incy;
+    ptrdiff_t const n_ = n;
+    ptrdiff_t const incx_ = incx;
+    ptrdiff_t const incy_ = incy;
 
     FORTRAN_WRAPPER(daxpy)(&n_, &alpha, x, &incx_, y, &incy_);
 }
diff --git a/src/model.cpp b/src/model.cpp
index 3c78802731..2485867a4c 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -19,7 +19,7 @@ namespace amici {
 /**
  * @brief Maps ModelQuantity items to their string value
  */
-const std::map<ModelQuantity, std::string> model_quantity_to_str{
+std::map<ModelQuantity, std::string> const model_quantity_to_str{
     {ModelQuantity::J, "J"},
     {ModelQuantity::JB, "JB"},
     {ModelQuantity::Jv, "Jv"},
@@ -996,7 +996,7 @@ void Model::setUnscaledInitialStateSensitivities(
     sx0data_ = sx0;
 }
 
-void Model::setSteadyStateComputationMode(const SteadyStateComputationMode mode
+void Model::setSteadyStateComputationMode(SteadyStateComputationMode const mode
 ) {
     steadystate_computation_mode_ = mode;
 }
@@ -1005,7 +1005,7 @@ SteadyStateComputationMode Model::getSteadyStateComputationMode() const {
     return steadystate_computation_mode_;
 }
 
-void Model::setSteadyStateSensitivityMode(const SteadyStateSensitivityMode mode
+void Model::setSteadyStateSensitivityMode(SteadyStateSensitivityMode const mode
 ) {
     steadystate_sensitivity_mode_ = mode;
 }
@@ -1046,14 +1046,14 @@ void Model::requireSensitivitiesForAllParameters() {
 }
 
 void Model::getExpression(
-    gsl::span<realtype> w, const realtype t, AmiVector const& x
+    gsl::span<realtype> w, realtype const t, AmiVector const& x
 ) {
     fw(t, computeX_pos(x));
     writeSlice(derived_state_.w_, w);
 }
 
 void Model::getObservable(
-    gsl::span<realtype> y, const realtype t, AmiVector const& x
+    gsl::span<realtype> y, realtype const t, AmiVector const& x
 ) {
     fy(t, x);
     writeSlice(derived_state_.y_, y);
@@ -1064,7 +1064,7 @@ ObservableScaling Model::getObservableScaling(int /*iy*/) const {
 }
 
 void Model::getObservableSensitivity(
-    gsl::span<realtype> sy, const realtype t, AmiVector const& x,
+    gsl::span<realtype> sy, realtype const t, AmiVector const& x,
     AmiVectorArray const& sx
 ) {
     if (!ny)
@@ -1205,14 +1205,14 @@ void Model::getAdjointStateObservableUpdate(
 }
 
 void Model::getEvent(
-    gsl::span<realtype> z, int const ie, const realtype t, AmiVector const& x
+    gsl::span<realtype> z, int const ie, realtype const t, AmiVector const& x
 ) {
     fz(ie, t, x);
     writeSliceEvent(derived_state_.z_, z, ie);
 }
 
 void Model::getEventSensitivity(
-    gsl::span<realtype> sz, int const ie, const realtype t, AmiVector const& x,
+    gsl::span<realtype> sz, int const ie, realtype const t, AmiVector const& x,
     AmiVectorArray const& sx
 ) {
     if (pythonGenerated) {
@@ -1263,14 +1263,14 @@ void Model::getUnobservedEventSensitivity(
 }
 
 void Model::getEventRegularization(
-    gsl::span<realtype> rz, int const ie, const realtype t, AmiVector const& x
+    gsl::span<realtype> rz, int const ie, realtype const t, AmiVector const& x
 ) {
     frz(ie, t, x);
     writeSliceEvent(derived_state_.rz_, rz, ie);
 }
 
 void Model::getEventRegularizationSensitivity(
-    gsl::span<realtype> srz, int const ie, const realtype t, AmiVector const& x,
+    gsl::span<realtype> srz, int const ie, realtype const t, AmiVector const& x,
     AmiVectorArray const& sx
 ) {
     if (pythonGenerated) {
@@ -1313,7 +1313,7 @@ void Model::getEventRegularizationSensitivity(
 
 void Model::getEventSigma(
     gsl::span<realtype> sigmaz, int const ie, int const nroots,
-    const realtype t, ExpData const* edata
+    realtype const t, ExpData const* edata
 ) {
     fsigmaz(ie, nroots, t, edata);
     writeSliceEvent(derived_state_.sigmaz_, sigmaz, ie);
@@ -1321,14 +1321,14 @@ void Model::getEventSigma(
 
 void Model::getEventSigmaSensitivity(
     gsl::span<realtype> ssigmaz, int const ie, int const nroots,
-    const realtype t, ExpData const* edata
+    realtype const t, ExpData const* edata
 ) {
     fdsigmazdp(ie, nroots, t, edata);
     writeSensitivitySliceEvent(derived_state_.dsigmazdp_, ssigmaz, ie);
 }
 
 void Model::addEventObjective(
-    realtype& Jz, int const ie, int const nroots, const realtype t,
+    realtype& Jz, int const ie, int const nroots, realtype const t,
     AmiVector const& x, ExpData const& edata
 ) {
     fz(ie, t, x);
@@ -1348,7 +1348,7 @@ void Model::addEventObjective(
 }
 
 void Model::addEventObjectiveRegularization(
-    realtype& Jrz, int const ie, int const nroots, const realtype t,
+    realtype& Jrz, int const ie, int const nroots, realtype const t,
     AmiVector const& x, ExpData const& edata
 ) {
     frz(ie, t, x);
@@ -1370,7 +1370,7 @@ void Model::addEventObjectiveRegularization(
 
 void Model::addEventObjectiveSensitivity(
     std::vector<realtype>& sllh, std::vector<realtype>& s2llh, int const ie,
-    int const nroots, const realtype t, AmiVector const& x,
+    int const nroots, realtype const t, AmiVector const& x,
     AmiVectorArray const& sx, ExpData const& edata
 ) {
 
@@ -1405,7 +1405,7 @@ void Model::addEventObjectiveSensitivity(
 }
 
 void Model::getAdjointStateEventUpdate(
-    gsl::span<realtype> dJzdx, int const ie, int const nroots, const realtype t,
+    gsl::span<realtype> dJzdx, int const ie, int const nroots, realtype const t,
     AmiVector const& x, ExpData const& edata
 ) {
     fdJzdx(ie, nroots, t, x, edata);
@@ -1414,7 +1414,7 @@ void Model::getAdjointStateEventUpdate(
 
 void Model::addPartialEventObjectiveSensitivity(
     std::vector<realtype>& sllh, std::vector<realtype>& s2llh, int const ie,
-    int const nroots, const realtype t, AmiVector const& x, ExpData const& edata
+    int const nroots, realtype const t, AmiVector const& x, ExpData const& edata
 ) {
     if (!nz)
         return;
@@ -1425,7 +1425,7 @@ void Model::addPartialEventObjectiveSensitivity(
 }
 
 void Model::getEventTimeSensitivity(
-    std::vector<realtype>& stau, const realtype t, int const ie,
+    std::vector<realtype>& stau, realtype const t, int const ie,
     AmiVector const& x, AmiVectorArray const& sx
 ) {
 
@@ -1441,7 +1441,7 @@ void Model::getEventTimeSensitivity(
 }
 
 void Model::addStateEventUpdate(
-    AmiVector& x, int const ie, const realtype t, AmiVector const& xdot,
+    AmiVector& x, int const ie, realtype const t, AmiVector const& xdot,
     AmiVector const& xdot_old
 ) {
 
@@ -1465,7 +1465,7 @@ void Model::addStateEventUpdate(
 }
 
 void Model::addStateSensitivityEventUpdate(
-    AmiVectorArray& sx, int const ie, const realtype t, AmiVector const& x_old,
+    AmiVectorArray& sx, int const ie, realtype const t, AmiVector const& x_old,
     AmiVector const& xdot, AmiVector const& xdot_old,
     std::vector<realtype> const& stau
 ) {
@@ -1497,7 +1497,7 @@ void Model::addStateSensitivityEventUpdate(
 }
 
 void Model::addAdjointStateEventUpdate(
-    AmiVector& xB, int const ie, const realtype t, AmiVector const& x,
+    AmiVector& xB, int const ie, realtype const t, AmiVector const& x,
     AmiVector const& xdot, AmiVector const& xdot_old
 ) {
 
@@ -1522,7 +1522,7 @@ void Model::addAdjointStateEventUpdate(
 }
 
 void Model::addAdjointQuadratureEventUpdate(
-    AmiVector xQB, int const ie, const realtype t, AmiVector const& x,
+    AmiVector xQB, int const ie, realtype const t, AmiVector const& x,
     AmiVector const& xB, AmiVector const& xdot, AmiVector const& xdot_old
 ) {
     for (int ip = 0; ip < nplist(); ip++) {
@@ -2032,7 +2032,7 @@ void Model::initializeVectors() {
         derived_state_.dxdotdp = AmiVectorArray(nx_solver, nplist());
 }
 
-void Model::fy(const realtype t, AmiVector const& x) {
+void Model::fy(realtype const t, AmiVector const& x) {
     if (!ny)
         return;
 
@@ -2052,7 +2052,7 @@ void Model::fy(const realtype t, AmiVector const& x) {
     }
 }
 
-void Model::fdydp(const realtype t, AmiVector const& x) {
+void Model::fdydp(realtype const t, AmiVector const& x) {
     if (!ny)
         return;
 
@@ -2086,7 +2086,7 @@ void Model::fdydp(const realtype t, AmiVector const& x) {
     }
 }
 
-void Model::fdydx(const realtype t, AmiVector const& x) {
+void Model::fdydx(realtype const t, AmiVector const& x) {
     if (!ny)
         return;
 
@@ -2420,7 +2420,7 @@ void Model::fdJydx(int const it, AmiVector const& x, ExpData const& edata) {
     }
 }
 
-void Model::fz(int const ie, const realtype t, AmiVector const& x) {
+void Model::fz(int const ie, realtype const t, AmiVector const& x) {
 
     derived_state_.z_.assign(nz, 0.0);
 
@@ -2429,7 +2429,7 @@ void Model::fz(int const ie, const realtype t, AmiVector const& x) {
        state_.h.data());
 }
 
-void Model::fdzdp(int const ie, const realtype t, AmiVector const& x) {
+void Model::fdzdp(int const ie, realtype const t, AmiVector const& x) {
     if (!nz)
         return;
 
@@ -2448,7 +2448,7 @@ void Model::fdzdp(int const ie, const realtype t, AmiVector const& x) {
     }
 }
 
-void Model::fdzdx(int const ie, const realtype t, AmiVector const& x) {
+void Model::fdzdx(int const ie, realtype const t, AmiVector const& x) {
     if (!nz)
         return;
 
@@ -2465,7 +2465,7 @@ void Model::fdzdx(int const ie, const realtype t, AmiVector const& x) {
     }
 }
 
-void Model::frz(int const ie, const realtype t, AmiVector const& x) {
+void Model::frz(int const ie, realtype const t, AmiVector const& x) {
 
     derived_state_.rz_.assign(nz, 0.0);
 
@@ -2474,7 +2474,7 @@ void Model::frz(int const ie, const realtype t, AmiVector const& x) {
         state_.h.data());
 }
 
-void Model::fdrzdp(int const ie, const realtype t, AmiVector const& x) {
+void Model::fdrzdp(int const ie, realtype const t, AmiVector const& x) {
     if (!nz)
         return;
 
@@ -2493,7 +2493,7 @@ void Model::fdrzdp(int const ie, const realtype t, AmiVector const& x) {
     }
 }
 
-void Model::fdrzdx(int const ie, const realtype t, AmiVector const& x) {
+void Model::fdrzdx(int const ie, realtype const t, AmiVector const& x) {
     if (!nz)
         return;
 
@@ -2511,7 +2511,7 @@ void Model::fdrzdx(int const ie, const realtype t, AmiVector const& x) {
 }
 
 void Model::fsigmaz(
-    int const ie, int const nroots, const realtype t, ExpData const* edata
+    int const ie, int const nroots, realtype const t, ExpData const* edata
 ) {
     if (!nz)
         return;
@@ -2547,7 +2547,7 @@ void Model::fsigmaz(
 }
 
 void Model::fdsigmazdp(
-    int const ie, int const nroots, const realtype t, ExpData const* edata
+    int const ie, int const nroots, realtype const t, ExpData const* edata
 ) {
     if (!nz)
         return;
@@ -2583,7 +2583,7 @@ void Model::fdsigmazdp(
 }
 
 void Model::fdJzdz(
-    int const ie, int const nroots, const realtype t, AmiVector const& x,
+    int const ie, int const nroots, realtype const t, AmiVector const& x,
     ExpData const& edata
 ) {
     if (!nz)
@@ -2615,7 +2615,7 @@ void Model::fdJzdz(
 }
 
 void Model::fdJzdsigma(
-    int const ie, int const nroots, const realtype t, AmiVector const& x,
+    int const ie, int const nroots, realtype const t, AmiVector const& x,
     ExpData const& edata
 ) {
     if (!nz)
@@ -2715,7 +2715,7 @@ void Model::fdJzdp(
 }
 
 void Model::fdJzdx(
-    int const ie, int const nroots, const realtype t, AmiVector const& x,
+    int const ie, int const nroots, realtype const t, AmiVector const& x,
     ExpData const& edata
 ) {
     // dJzdz         nJ x nz        x nztrue
@@ -2763,7 +2763,7 @@ void Model::fdJzdx(
 }
 
 void Model::fdJrzdz(
-    int const ie, int const nroots, const realtype t, AmiVector const& x,
+    int const ie, int const nroots, realtype const t, AmiVector const& x,
     ExpData const& edata
 ) {
     if (!nz)
@@ -2794,7 +2794,7 @@ void Model::fdJrzdz(
 }
 
 void Model::fdJrzdsigma(
-    int const ie, int const nroots, const realtype t, AmiVector const& x,
+    int const ie, int const nroots, realtype const t, AmiVector const& x,
     ExpData const& edata
 ) {
     if (!nz)
@@ -2825,12 +2825,12 @@ void Model::fdJrzdsigma(
     }
 }
 
-void Model::fspl(const realtype t) {
+void Model::fspl(realtype const t) {
     for (int ispl = 0; ispl < nspl; ispl++)
         state_.spl_[ispl] = splines_[ispl].get_value(t);
 }
 
-void Model::fsspl(const realtype t) {
+void Model::fsspl(realtype const t) {
     derived_state_.sspl_.zero();
     realtype* sspl_data = derived_state_.sspl_.data();
     for (int ip = 0; ip < nplist(); ip++) {
@@ -2840,7 +2840,7 @@ void Model::fsspl(const realtype t) {
     }
 }
 
-void Model::fw(const realtype t, realtype const* x) {
+void Model::fw(realtype const t, realtype const* x) {
     std::fill(derived_state_.w_.begin(), derived_state_.w_.end(), 0.0);
     fspl(t);
     fw(derived_state_.w_.data(), t, x, state_.unscaledParameters.data(),
@@ -2852,7 +2852,7 @@ void Model::fw(const realtype t, realtype const* x) {
     }
 }
 
-void Model::fdwdp(const realtype t, realtype const* x) {
+void Model::fdwdp(realtype const t, realtype const* x) {
     if (!nw)
         return;
 
@@ -2901,7 +2901,7 @@ void Model::fdwdp(const realtype t, realtype const* x) {
     }
 }
 
-void Model::fdwdx(const realtype t, realtype const* x) {
+void Model::fdwdx(realtype const t, realtype const* x) {
     if (!nw)
         return;
 
@@ -2947,7 +2947,7 @@ void Model::fdwdx(const realtype t, realtype const* x) {
     }
 }
 
-void Model::fdwdw(const realtype t, realtype const* x) {
+void Model::fdwdw(realtype const t, realtype const* x) {
     if (!nw || !dwdw_.capacity())
         return;
     dwdw_.zero();
diff --git a/src/model_dae.cpp b/src/model_dae.cpp
index 43a8d81313..3b2e74e0e1 100644
--- a/src/model_dae.cpp
+++ b/src/model_dae.cpp
@@ -4,7 +4,7 @@
 namespace amici {
 
 void Model_DAE::fJ(
-    const realtype t, const realtype cj, AmiVector const& x,
+    realtype const t, realtype const cj, AmiVector const& x,
     AmiVector const& dx, AmiVector const& xdot, SUNMatrix J
 ) {
     fJ(t, cj, x.getNVector(), dx.getNVector(), xdot.getNVector(), J);
@@ -21,7 +21,7 @@ void Model_DAE::fJ(
 }
 
 void Model_DAE::fJSparse(
-    const realtype t, const realtype cj, AmiVector const& x,
+    realtype const t, realtype const cj, AmiVector const& x,
     AmiVector const& dx, AmiVector const& /*xdot*/, SUNMatrix J
 ) {
     fJSparse(t, cj, x.getNVector(), dx.getNVector(), J);
@@ -75,9 +75,9 @@ void Model_DAE::fJSparse(
 }
 
 void Model_DAE::fJv(
-    const realtype t, AmiVector const& x, AmiVector const& dx,
+    realtype const t, AmiVector const& x, AmiVector const& dx,
     AmiVector const& /*xdot*/, AmiVector const& v, AmiVector& Jv,
-    const realtype cj
+    realtype const cj
 ) {
     fJv(t, x.getNVector(), dx.getNVector(), v.getNVector(), Jv.getNVector(),
         cj);
@@ -94,7 +94,7 @@ void Model_DAE::fJv(
 }
 
 void Model_DAE::froot(
-    const realtype t, AmiVector const& x, AmiVector const& dx,
+    realtype const t, AmiVector const& x, AmiVector const& dx,
     gsl::span<realtype> root
 ) {
     froot(t, x.getNVector(), dx.getNVector(), root);
@@ -113,7 +113,7 @@ void Model_DAE::froot(
 }
 
 void Model_DAE::fxdot(
-    const realtype t, AmiVector const& x, AmiVector const& dx, AmiVector& xdot
+    realtype const t, AmiVector const& x, AmiVector const& dx, AmiVector& xdot
 ) {
     fxdot(t, x.getNVector(), dx.getNVector(), xdot.getNVector());
 }
@@ -132,7 +132,7 @@ void Model_DAE::fxdot(
 }
 
 void Model_DAE::fJDiag(
-    const realtype t, AmiVector& JDiag, const realtype /*cj*/,
+    realtype const t, AmiVector& JDiag, realtype const /*cj*/,
     AmiVector const& x, AmiVector const& dx
 ) {
     fJSparse(t, 0.0, x.getNVector(), dx.getNVector(), derived_state_.J_.get());
@@ -143,7 +143,7 @@ void Model_DAE::fJDiag(
 }
 
 void Model_DAE::fdxdotdw(
-    const realtype t, const_N_Vector x, const const_N_Vector dx
+    realtype const t, const_N_Vector x, const_N_Vector const dx
 ) {
     derived_state_.dxdotdw_.zero();
     if (nw > 0 && derived_state_.dxdotdw_.capacity()) {
@@ -161,7 +161,7 @@ void Model_DAE::fdxdotdw(
 }
 
 void Model_DAE::fdxdotdp(
-    const realtype t, const const_N_Vector x, const const_N_Vector dx
+    realtype const t, const_N_Vector const x, const_N_Vector const dx
 ) {
     auto x_pos = computeX_pos(x);
 
@@ -230,7 +230,7 @@ void Model_DAE::fJSparse(
 }
 
 void Model_DAE::froot(
-    realtype* /*root*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*root*/, realtype const /*t*/, realtype const* /*x*/,
     double const* /*p*/, double const* /*k*/, realtype const* /*h*/,
     realtype const* /*dx*/
 ) {
@@ -242,7 +242,7 @@ void Model_DAE::froot(
 }
 
 void Model_DAE::fdxdotdp(
-    realtype* /*dxdotdp*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dxdotdp*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     int const /*ip*/, realtype const* /*dx*/, realtype const* /*w*/,
     realtype const* /*dwdp*/
@@ -255,7 +255,7 @@ void Model_DAE::fdxdotdp(
 }
 
 void Model_DAE::fdxdotdp_explicit(
-    realtype* /*dxdotdp_explicit*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dxdotdp_explicit*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*dx*/, realtype const* /*w*/
 ) {
@@ -283,7 +283,7 @@ void Model_DAE::fdxdotdp_explicit_rowvals(SUNMatrixWrapper& /*dxdotdp*/) {
 }
 
 void Model_DAE::fdxdotdx_explicit(
-    realtype* /*dxdotdx_explicit*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dxdotdx_explicit*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*dx*/, realtype const* /*w*/
 ) {
@@ -311,7 +311,7 @@ void Model_DAE::fdxdotdx_explicit_rowvals(SUNMatrixWrapper& /*dxdotdx*/) {
 }
 
 void Model_DAE::fdxdotdw(
-    realtype* /*dxdotdw*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dxdotdw*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*dx*/, realtype const* /*w*/
 ) {
@@ -339,11 +339,11 @@ void Model_DAE::fdxdotdw_rowvals(SUNMatrixWrapper& /*dxdotdw*/) {
 }
 
 void Model_DAE::
-    fM(realtype* /*M*/, const realtype /*t*/, realtype const* /*x*/,
+    fM(realtype* /*M*/, realtype const /*t*/, realtype const* /*x*/,
        realtype const* /*p*/, realtype const* /*k*/) {}
 
 void Model_DAE::fJB(
-    const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+    realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
     AmiVector const& xB, AmiVector const& /*dxB*/, AmiVector const& /*xBdot*/,
     SUNMatrix JB
 ) {
@@ -362,7 +362,7 @@ void Model_DAE::fJB(
 }
 
 void Model_DAE::fJSparseB(
-    const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+    realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
     AmiVector const& xB, AmiVector const& dxB, AmiVector const& /*xBdot*/,
     SUNMatrix JB
 ) {
@@ -427,7 +427,7 @@ void Model_DAE::fqBdot(
 }
 
 void Model_DAE::fxBdot_ss(
-    const realtype t, AmiVector const& xB, AmiVector const& dxB,
+    realtype const t, AmiVector const& xB, AmiVector const& dxB,
     AmiVector& xBdot
 ) {
     fxBdot_ss(t, xB.getNVector(), dxB.getNVector(), xBdot.getNVector());
@@ -459,7 +459,7 @@ void Model_DAE::fJSparseB_ss(SUNMatrix JB) {
 }
 
 void Model_DAE::writeSteadystateJB(
-    const realtype t, realtype cj, AmiVector const& x, AmiVector const& dx,
+    realtype const t, realtype cj, AmiVector const& x, AmiVector const& dx,
     AmiVector const& xB, AmiVector const& dxB, AmiVector const& /*xBdot*/
 ) {
     /* Get backward Jacobian */
@@ -473,7 +473,7 @@ void Model_DAE::writeSteadystateJB(
 }
 
 void Model_DAE::fsxdot(
-    const realtype t, AmiVector const& x, AmiVector const& dx, int const ip,
+    realtype const t, AmiVector const& x, AmiVector const& dx, int const ip,
     AmiVector const& sx, AmiVector const& sdx, AmiVector& sxdot
 ) {
     fsxdot(
diff --git a/src/model_ode.cpp b/src/model_ode.cpp
index 24787df8ad..b3d517b36e 100644
--- a/src/model_ode.cpp
+++ b/src/model_ode.cpp
@@ -5,7 +5,7 @@
 namespace amici {
 
 void Model_ODE::fJ(
-    const realtype t, const realtype /*cj*/, AmiVector const& x,
+    realtype const t, realtype const /*cj*/, AmiVector const& x,
     AmiVector const& /*dx*/, AmiVector const& xdot, SUNMatrix J
 ) {
     fJ(t, x.getNVector(), xdot.getNVector(), J);
@@ -23,7 +23,7 @@ void Model_ODE::fJ(
 }
 
 void Model_ODE::fJSparse(
-    const realtype t, const realtype /*cj*/, AmiVector const& x,
+    realtype const t, realtype const /*cj*/, AmiVector const& x,
     AmiVector const& /*dx*/, AmiVector const& /*xdot*/, SUNMatrix J
 ) {
     fJSparse(t, x.getNVector(), J);
@@ -69,9 +69,9 @@ void Model_ODE::fJSparse(realtype t, const_N_Vector x, SUNMatrix J) {
 }
 
 void Model_ODE::
-    fJv(const realtype t, AmiVector const& x, AmiVector const& /*dx*/,
+    fJv(realtype const t, AmiVector const& x, AmiVector const& /*dx*/,
         AmiVector const& /*xdot*/, AmiVector const& v, AmiVector& Jv,
-        const realtype /*cj*/) {
+        realtype const /*cj*/) {
     fJv(v.getNVector(), Jv.getNVector(), t, x.getNVector());
 }
 
@@ -85,7 +85,7 @@ void Model_ODE::fJv(
 }
 
 void Model_ODE::froot(
-    const realtype t, AmiVector const& x, AmiVector const& /*dx*/,
+    realtype const t, AmiVector const& x, AmiVector const& /*dx*/,
     gsl::span<realtype> root
 ) {
     froot(t, x.getNVector(), root);
@@ -102,7 +102,7 @@ void Model_ODE::froot(realtype t, const_N_Vector x, gsl::span<realtype> root) {
 }
 
 void Model_ODE::fxdot(
-    const realtype t, AmiVector const& x, AmiVector const& /*dx*/,
+    realtype const t, AmiVector const& x, AmiVector const& /*dx*/,
     AmiVector& xdot
 ) {
     fxdot(t, x.getNVector(), xdot.getNVector());
@@ -120,7 +120,7 @@ void Model_ODE::fxdot(realtype t, const_N_Vector x, N_Vector xdot) {
 }
 
 void Model_ODE::fJDiag(
-    const realtype t, AmiVector& JDiag, const realtype /*cj*/,
+    realtype const t, AmiVector& JDiag, realtype const /*cj*/,
     AmiVector const& x, AmiVector const& /*dx*/
 ) {
     fJDiag(t, JDiag.getNVector(), x.getNVector());
@@ -128,7 +128,7 @@ void Model_ODE::fJDiag(
         throw AmiException("Evaluation of fJDiag failed!");
 }
 
-void Model_ODE::fdxdotdw(const realtype t, const_N_Vector x) {
+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);
@@ -143,7 +143,7 @@ void Model_ODE::fdxdotdw(const realtype t, const_N_Vector x) {
     }
 }
 
-void Model_ODE::fdxdotdp(const realtype t, const_N_Vector x) {
+void Model_ODE::fdxdotdp(realtype const t, const_N_Vector x) {
     auto x_pos = computeX_pos(x);
     fdwdp(t, N_VGetArrayPointerConst(x_pos));
 
@@ -189,7 +189,7 @@ void Model_ODE::fdxdotdp(const realtype t, const_N_Vector x) {
 }
 
 void Model_ODE::
-    fdxdotdp(const realtype t, AmiVector const& x, AmiVector const& /*dx*/) {
+    fdxdotdp(realtype const t, AmiVector const& x, AmiVector const& /*dx*/) {
     fdxdotdp(t, x.getNVector());
 }
 
@@ -198,7 +198,7 @@ std::unique_ptr<Solver> Model_ODE::getSolver() {
 }
 
 void Model_ODE::fJSparse(
-    SUNMatrixContent_Sparse /*JSparse*/, const realtype /*t*/,
+    SUNMatrixContent_Sparse /*JSparse*/, realtype const /*t*/,
     realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/,
     realtype const* /*h*/, realtype const* /*w*/, realtype const* /*dwdx*/
 ) {
@@ -210,7 +210,7 @@ void Model_ODE::fJSparse(
 }
 
 void Model_ODE::fJSparse(
-    realtype* /*JSparse*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*JSparse*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/, realtype const* /*dwdx*/
 ) {
@@ -238,7 +238,7 @@ void Model_ODE::fJSparse_rowvals(SUNMatrixWrapper& /*JSparse*/) {
 }
 
 void Model_ODE::froot(
-    realtype* /*root*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*root*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*tcl*/
 ) {
@@ -250,7 +250,7 @@ void Model_ODE::froot(
 }
 
 void Model_ODE::fdxdotdp(
-    realtype* /*dxdotdp*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dxdotdp*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     int const /*ip*/, realtype const* /*w*/, realtype const* /*dwdp*/
 ) {
@@ -262,7 +262,7 @@ void Model_ODE::fdxdotdp(
 }
 
 void Model_ODE::fdxdotdp_explicit(
-    realtype* /*dxdotdp_explicit*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dxdotdp_explicit*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/
 ) {
@@ -290,7 +290,7 @@ void Model_ODE::fdxdotdp_explicit_rowvals(SUNMatrixWrapper& /*dxdotdp*/) {
 }
 
 void Model_ODE::fdxdotdx_explicit(
-    realtype* /*dxdotdx_explicit*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dxdotdx_explicit*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/
 ) {
@@ -318,7 +318,7 @@ void Model_ODE::fdxdotdx_explicit_rowvals(SUNMatrixWrapper& /*dxdotdx*/) {
 }
 
 void Model_ODE::fdxdotdw(
-    realtype* /*dxdotdw*/, const realtype /*t*/, realtype const* /*x*/,
+    realtype* /*dxdotdw*/, realtype const /*t*/, realtype const* /*x*/,
     realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/,
     realtype const* /*w*/
 ) {
@@ -346,7 +346,7 @@ void Model_ODE::fdxdotdw_rowvals(SUNMatrixWrapper& /*dxdotdw*/) {
 }
 
 void Model_ODE::fJB(
-    const realtype t, realtype /*cj*/, AmiVector const& x,
+    realtype const t, realtype /*cj*/, AmiVector const& x,
     AmiVector const& /*dx*/, AmiVector const& xB, AmiVector const& /*dxB*/,
     AmiVector const& xBdot, SUNMatrix JB
 ) {
@@ -364,7 +364,7 @@ void Model_ODE::fJB(
 }
 
 void Model_ODE::fJSparseB(
-    const realtype t, realtype /*cj*/, AmiVector const& x,
+    realtype const t, realtype /*cj*/, AmiVector const& x,
     AmiVector const& /*dx*/, AmiVector const& xB, AmiVector const& /*dxB*/,
     AmiVector const& xBdot, SUNMatrix JB
 ) {
@@ -436,7 +436,7 @@ void Model_ODE::fqBdot(
 }
 
 void Model_ODE::fxBdot_ss(
-    const realtype t, AmiVector const& xB, AmiVector const& /*dx*/,
+    realtype const t, AmiVector const& xB, AmiVector const& /*dx*/,
     AmiVector& xBdot
 ) {
     fxBdot_ss(t, xB.getNVector(), xBdot.getNVector());
@@ -463,7 +463,7 @@ void Model_ODE::fJSparseB_ss(SUNMatrix JB) {
 }
 
 void Model_ODE::writeSteadystateJB(
-    const realtype t, realtype /*cj*/, AmiVector const& x,
+    realtype const t, realtype /*cj*/, AmiVector const& x,
     AmiVector const& /*dx*/, AmiVector const& xB, AmiVector const& /*dxB*/,
     AmiVector const& xBdot
 ) {
@@ -478,7 +478,7 @@ void Model_ODE::writeSteadystateJB(
 }
 
 void Model_ODE::fsxdot(
-    const realtype t, AmiVector const& x, AmiVector const& /*dx*/, int const ip,
+    realtype const t, AmiVector const& x, AmiVector const& /*dx*/, int const ip,
     AmiVector const& sx, AmiVector const& /*sdx*/, AmiVector& sxdot
 ) {
     fsxdot(t, x.getNVector(), ip, sx.getNVector(), sxdot.getNVector());
diff --git a/src/rdata.cpp b/src/rdata.cpp
index e96f295d23..d7b99f3c8a 100644
--- a/src/rdata.cpp
+++ b/src/rdata.cpp
@@ -707,7 +707,7 @@ void ReturnData::applyChainRuleFactorToSimulationResults(Model const& model) {
         for (int IND1 = 0; (IND1) < (N1T); ++(IND1))                           \
             for (int ip = 0; ip < nplist; ++ip)                                \
                 for (int IND2 = 0; (IND2) < (N2); ++(IND2)) {                  \
-                    s##QUANT.at(((IND2)*nplist + ip) * (N1) + (IND1))          \
+                    s##QUANT.at(((IND2) * nplist + ip) * (N1) + (IND1))        \
                         *= pcoefficient.at(ip);                                \
                 }
 
diff --git a/src/solver.cpp b/src/solver.cpp
index 22e1723640..1a74e8b44b 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -81,7 +81,7 @@ void Solver::apply_max_num_steps_B() const {
     }
 }
 
-int Solver::run(const realtype tout) const {
+int Solver::run(realtype const tout) const {
     setStopTime(tout);
     CpuTimer cpu_timer;
     int status = AMICI_SUCCESS;
@@ -100,7 +100,7 @@ int Solver::run(const realtype tout) const {
     return status;
 }
 
-int Solver::step(const realtype tout) const {
+int Solver::step(realtype const tout) const {
     int status = AMICI_SUCCESS;
 
     apply_max_num_steps();
@@ -116,7 +116,7 @@ int Solver::step(const realtype tout) const {
     return status;
 }
 
-void Solver::runB(const realtype tout) const {
+void Solver::runB(realtype const tout) const {
     CpuTimer cpu_timer;
 
     apply_max_num_steps_B();
@@ -128,7 +128,7 @@ void Solver::runB(const realtype tout) const {
 }
 
 void Solver::setup(
-    const realtype t0, Model* model, AmiVector const& x0, AmiVector const& dx0,
+    realtype const t0, Model* model, AmiVector const& x0, AmiVector const& dx0,
     AmiVectorArray const& sx0, AmiVectorArray const& sdx0
 ) const {
     if (nx() != model->nx_solver || nplist() != model->nplist()
@@ -196,7 +196,7 @@ void Solver::setup(
 }
 
 void Solver::setupB(
-    int* which, const realtype tf, Model* model, AmiVector const& xB0,
+    int* which, realtype const tf, Model* model, AmiVector const& xB0,
     AmiVector const& dxB0, AmiVector const& xQB0
 ) const {
     if (!solver_memory_)
@@ -228,7 +228,7 @@ void Solver::setupB(
 }
 
 void Solver::setupSteadystate(
-    const realtype t0, Model* model, AmiVector const& x0, AmiVector const& dx0,
+    realtype const t0, Model* model, AmiVector const& x0, AmiVector const& dx0,
     AmiVector const& xB0, AmiVector const& dxB0, AmiVector const& xQ0
 ) const {
     /* Initialize CVodes/IDAs solver with steadystate RHS function */
@@ -642,20 +642,20 @@ SensitivityMethod Solver::getSensitivityMethodPreequilibration() const {
     return sensi_meth_preeq_;
 }
 
-void Solver::setSensitivityMethod(const SensitivityMethod sensi_meth) {
+void Solver::setSensitivityMethod(SensitivityMethod const sensi_meth) {
     checkSensitivityMethod(sensi_meth, false);
     this->sensi_meth_ = sensi_meth;
 }
 
 void Solver::setSensitivityMethodPreequilibration(
-    const SensitivityMethod sensi_meth_preeq
+    SensitivityMethod const sensi_meth_preeq
 ) {
     checkSensitivityMethod(sensi_meth_preeq, true);
     sensi_meth_preeq_ = sensi_meth_preeq;
 }
 
 void Solver::checkSensitivityMethod(
-    const SensitivityMethod sensi_meth, bool preequilibration
+    SensitivityMethod const sensi_meth, bool preequilibration
 ) const {
     if (rdata_mode_ == RDataReporting::residuals
         && sensi_meth == SensitivityMethod::adjoint)
@@ -693,7 +693,7 @@ void Solver::setNewtonDampingFactorLowerBound(double dampingFactorLowerBound) {
 
 SensitivityOrder Solver::getSensitivityOrder() const { return sensi_; }
 
-void Solver::setSensitivityOrder(const SensitivityOrder sensi) {
+void Solver::setSensitivityOrder(SensitivityOrder const sensi) {
     if (sensi_ != sensi)
         resetMutableMemory(nx(), nplist(), nquad());
     sensi_ = sensi;
@@ -950,7 +950,7 @@ void Solver::setMaxStepsBackwardProblem(long int const maxsteps) {
 
 LinearMultistepMethod Solver::getLinearMultistepMethod() const { return lmm_; }
 
-void Solver::setLinearMultistepMethod(const LinearMultistepMethod lmm) {
+void Solver::setLinearMultistepMethod(LinearMultistepMethod const lmm) {
     if (solver_memory_)
         resetMutableMemory(nx(), nplist(), nquad());
     lmm_ = lmm;
@@ -960,7 +960,7 @@ NonlinearSolverIteration Solver::getNonlinearSolverIteration() const {
     return iter_;
 }
 
-void Solver::setNonlinearSolverIteration(const NonlinearSolverIteration iter) {
+void Solver::setNonlinearSolverIteration(NonlinearSolverIteration const iter) {
     if (solver_memory_)
         resetMutableMemory(nx(), nplist(), nquad());
     iter_ = iter;
@@ -968,7 +968,7 @@ void Solver::setNonlinearSolverIteration(const NonlinearSolverIteration iter) {
 
 InterpolationType Solver::getInterpolationType() const { return interp_type_; }
 
-void Solver::setInterpolationType(const InterpolationType interpType) {
+void Solver::setInterpolationType(InterpolationType const interpType) {
     if (!solver_memory_B_.empty())
         resetMutableMemory(nx(), nplist(), nquad());
     interp_type_ = interpType;
@@ -1020,7 +1020,7 @@ InternalSensitivityMethod Solver::getInternalSensitivityMethod() const {
     return ism_;
 }
 
-void Solver::setInternalSensitivityMethod(const InternalSensitivityMethod ism) {
+void Solver::setInternalSensitivityMethod(InternalSensitivityMethod const ism) {
     if (solver_memory_)
         resetMutableMemory(nx(), nplist(), nquad());
     ism_ = ism;
@@ -1181,7 +1181,7 @@ void Solver::writeSolutionB(
     xQB.copy(getAdjointQuadrature(which, *t));
 }
 
-AmiVector const& Solver::getState(const realtype t) const {
+AmiVector const& Solver::getState(realtype const t) const {
     if (t == t_)
         return x_;
 
@@ -1191,7 +1191,7 @@ AmiVector const& Solver::getState(const realtype t) const {
     return dky_;
 }
 
-AmiVector const& Solver::getDerivativeState(const realtype t) const {
+AmiVector const& Solver::getDerivativeState(realtype const t) const {
     if (t == t_)
         return dx_;
 
@@ -1201,7 +1201,7 @@ AmiVector const& Solver::getDerivativeState(const realtype t) const {
     return dky_;
 }
 
-AmiVectorArray const& Solver::getStateSensitivity(const realtype t) const {
+AmiVectorArray const& Solver::getStateSensitivity(realtype const t) const {
     if (sens_initialized_ && solver_was_called_F_) {
         if (t == t_) {
             getSens();
@@ -1213,7 +1213,7 @@ AmiVectorArray const& Solver::getStateSensitivity(const realtype t) const {
 }
 
 AmiVector const&
-Solver::getAdjointState(int const which, const realtype t) const {
+Solver::getAdjointState(int const which, realtype const t) const {
     if (adj_initialized_) {
         if (solver_was_called_B_) {
             if (t == t_) {
@@ -1229,7 +1229,7 @@ Solver::getAdjointState(int const which, const realtype t) const {
 }
 
 AmiVector const&
-Solver::getAdjointDerivativeState(int const which, const realtype t) const {
+Solver::getAdjointDerivativeState(int const which, realtype const t) const {
     if (adj_initialized_) {
         if (solver_was_called_B_) {
             if (t == t_) {
@@ -1245,7 +1245,7 @@ Solver::getAdjointDerivativeState(int const which, const realtype t) const {
 }
 
 AmiVector const&
-Solver::getAdjointQuadrature(int const which, const realtype t) const {
+Solver::getAdjointQuadrature(int const which, realtype const t) const {
     if (adj_initialized_) {
         if (solver_was_called_B_) {
             if (t == t_) {
diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp
index b53be68c2e..f5150a87c7 100644
--- a/src/solver_cvodes.cpp
+++ b/src/solver_cvodes.cpp
@@ -104,7 +104,7 @@ static int fsxdot(
 /* Function implementations */
 
 void CVodeSolver::
-    init(const realtype t0, AmiVector const& x0, AmiVector const& /*dx0*/)
+    init(realtype const t0, AmiVector const& x0, AmiVector const& /*dx0*/)
         const {
     solver_was_called_F_ = false;
     force_reinit_postprocess_F_ = false;
@@ -122,7 +122,7 @@ void CVodeSolver::
 }
 
 void CVodeSolver::initSteadystate(
-    const realtype /*t0*/, AmiVector const& /*x0*/, AmiVector const& /*dx0*/
+    realtype const /*t0*/, AmiVector const& /*x0*/, AmiVector const& /*dx0*/
 ) const {
     // We need to set the steadystate rhs function. Sundials doesn't have this
     // in its public API, so we have to change it in the solver memory,
@@ -161,7 +161,7 @@ void CVodeSolver::
 }
 
 void CVodeSolver::binit(
-    int const which, const realtype tf, AmiVector const& xB0,
+    int const which, realtype const tf, AmiVector const& xB0,
     AmiVector const& /*dxB0*/
 ) const {
     solver_was_called_B_ = false;
@@ -452,12 +452,12 @@ void CVodeSolver::resetState(void* ami_mem, const_N_Vector y0) const {
     N_VScale(ONE, const_cast<N_Vector>(y0), cv_mem->cv_zn[0]);
 }
 
-void CVodeSolver::reInitPostProcessF(const realtype tnext) const {
+void CVodeSolver::reInitPostProcessF(realtype const tnext) const {
     reInitPostProcess(solver_memory_.get(), &t_, &x_, tnext);
     force_reinit_postprocess_F_ = false;
 }
 
-void CVodeSolver::reInitPostProcessB(const realtype tnext) const {
+void CVodeSolver::reInitPostProcessB(realtype const tnext) const {
     realtype tBret;
     auto cv_mem = static_cast<CVodeMem>(solver_memory_.get());
     auto ca_mem = cv_mem->cv_adj_mem;
@@ -477,7 +477,7 @@ void CVodeSolver::reInitPostProcessB(const realtype tnext) const {
 }
 
 void CVodeSolver::reInitPostProcess(
-    void* ami_mem, realtype* t, AmiVector* yout, const realtype tout
+    void* ami_mem, realtype* t, AmiVector* yout, realtype const tout
 ) const {
     auto cv_mem = static_cast<CVodeMem>(ami_mem);
     auto nst_tmp = cv_mem->cv_nst;
@@ -499,7 +499,7 @@ void CVodeSolver::reInitPostProcess(
         );
     if (status != CV_SUCCESS) {
         std::stringstream msg;
-        msg<<"tout: "<<tout<<", t: "<<*t<<".";
+        msg << "tout: " << tout << ", t: " << *t << ".";
         throw CvodeException(status, "reInitPostProcess", msg.str().c_str());
     }
 
@@ -531,7 +531,7 @@ void CVodeSolver::reInitPostProcess(
 }
 
 void CVodeSolver::
-    reInit(const realtype t0, AmiVector const& yy0, AmiVector const& /*yp0*/)
+    reInit(realtype const t0, AmiVector const& yy0, AmiVector const& /*yp0*/)
         const {
     auto cv_mem = static_cast<CVodeMem>(solver_memory_.get());
     cv_mem->cv_tn = t0;
@@ -559,7 +559,7 @@ void CVodeSolver::
 }
 
 void CVodeSolver::reInitB(
-    int const which, const realtype tB0, AmiVector const& yyB0,
+    int const which, realtype const tB0, AmiVector const& yyB0,
     AmiVector const& /*ypB0*/
 ) const {
     auto cv_memB = static_cast<CVodeMem>(
@@ -616,14 +616,14 @@ void CVodeSolver::getSens() const {
         throw CvodeException(status, "CVodeGetSens");
 }
 
-void CVodeSolver::getSensDky(const realtype t, int const k) const {
+void CVodeSolver::getSensDky(realtype const t, int const k) const {
     int status
         = CVodeGetSensDky(solver_memory_.get(), t, k, sx_.getNVectorArray());
     if (status != CV_SUCCESS)
         throw CvodeException(status, "CVodeGetSens");
 }
 
-void CVodeSolver::getDkyB(const realtype t, int const k, int const which)
+void CVodeSolver::getDkyB(realtype const t, int const k, int const which)
     const {
     int status = CVodeGetDky(
         CVodeGetAdjCVodeBmem(solver_memory_.get(), which), t, k,
@@ -648,7 +648,7 @@ void CVodeSolver::getQuad(realtype& t) const {
         throw CvodeException(status, "CVodeGetQuad");
 }
 
-void CVodeSolver::getQuadDkyB(const realtype t, int const k, int which) const {
+void CVodeSolver::getQuadDkyB(realtype const t, int const k, int which) const {
     int status = CVodeGetQuadDky(
         CVodeGetAdjCVodeBmem(solver_memory_.get(), which), t, k,
         xQB_.getNVector()
@@ -657,7 +657,7 @@ void CVodeSolver::getQuadDkyB(const realtype t, int const k, int which) const {
         throw CvodeException(status, "CVodeGetQuadDkyB");
 }
 
-void CVodeSolver::getQuadDky(const realtype t, int const k) const {
+void CVodeSolver::getQuadDky(realtype const t, int const k) const {
     int status = CVodeGetQuadDky(solver_memory_.get(), t, k, xQ_.getNVector());
     if (status != CV_SUCCESS)
         throw CvodeException(status, "CVodeGetQuadDky");
@@ -712,7 +712,7 @@ void CVodeSolver::allocateSolverB(int* which) const {
 }
 
 void CVodeSolver::setSStolerancesB(
-    int const which, const realtype relTolB, const realtype absTolB
+    int const which, realtype const relTolB, realtype const absTolB
 ) const {
     int status
         = CVodeSStolerancesB(solver_memory_.get(), which, relTolB, absTolB);
@@ -721,7 +721,7 @@ void CVodeSolver::setSStolerancesB(
 }
 
 void CVodeSolver::quadSStolerancesB(
-    int const which, const realtype reltolQB, const realtype abstolQB
+    int const which, realtype const reltolQB, realtype const abstolQB
 ) const {
     int status = CVodeQuadSStolerancesB(
         solver_memory_.get(), which, reltolQB, abstolQB
@@ -731,7 +731,7 @@ void CVodeSolver::quadSStolerancesB(
 }
 
 void CVodeSolver::quadSStolerances(
-    const realtype reltolQB, const realtype abstolQB
+    realtype const reltolQB, realtype const abstolQB
 ) const {
     int status
         = CVodeQuadSStolerances(solver_memory_.get(), reltolQB, abstolQB);
@@ -747,7 +747,7 @@ void CVodeSolver::getB(int const which) const {
         throw CvodeException(status, "CVodeGetB");
 }
 
-int CVodeSolver::solve(const realtype tout, int const itask) const {
+int CVodeSolver::solve(realtype const tout, int const itask) const {
     if (force_reinit_postprocess_F_)
         reInitPostProcessF(tout);
     int status = CVode(solver_memory_.get(), tout, x_.getNVector(), &t_, itask);
@@ -757,7 +757,7 @@ int CVodeSolver::solve(const realtype tout, int const itask) const {
     return status;
 }
 
-int CVodeSolver::solveF(const realtype tout, int const itask, int* ncheckPtr)
+int CVodeSolver::solveF(realtype const tout, int const itask, int* ncheckPtr)
     const {
     if (force_reinit_postprocess_F_)
         reInitPostProcessF(tout);
@@ -770,7 +770,7 @@ int CVodeSolver::solveF(const realtype tout, int const itask, int* ncheckPtr)
     return status;
 }
 
-void CVodeSolver::solveB(const realtype tBout, int const itaskB) const {
+void CVodeSolver::solveB(realtype const tBout, int const itaskB) const {
     if (force_reinit_postprocess_B_)
         reInitPostProcessB(tBout);
     int status = CVodeB(solver_memory_.get(), tBout, itaskB);
@@ -839,12 +839,12 @@ void* CVodeSolver::getAdjBmem(void* ami_mem, int which) const {
     return CVodeGetAdjCVodeBmem(ami_mem, which);
 }
 
-void CVodeSolver::calcIC(const realtype /*tout1*/) const {};
+void CVodeSolver::calcIC(realtype const /*tout1*/) const {};
 
-void CVodeSolver::calcICB(int const /*which*/, const realtype /*tout1*/)
-    const {};
+void CVodeSolver::calcICB(int const /*which*/, realtype const /*tout1*/) const {
+};
 
-void CVodeSolver::setStopTime(const realtype tstop) const {
+void CVodeSolver::setStopTime(realtype const tstop) const {
     int status = CVodeSetStopTime(solver_memory_.get(), tstop);
     if (status != CV_SUCCESS)
         throw CvodeException(status, "CVodeSetStopTime");
diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp
index 63f65b184d..dac085e98d 100644
--- a/src/solver_idas.cpp
+++ b/src/solver_idas.cpp
@@ -100,7 +100,7 @@ static int fsxdot(
 /* Function implementations */
 
 void IDASolver::init(
-    const realtype t0, AmiVector const& x0, AmiVector const& dx0
+    realtype const t0, AmiVector const& x0, AmiVector const& dx0
 ) const {
     int status;
     solver_was_called_F_ = false;
@@ -122,7 +122,7 @@ void IDASolver::init(
 }
 
 void IDASolver::initSteadystate(
-    const realtype /*t0*/, AmiVector const& /*x0*/, AmiVector const& /*dx0*/
+    realtype const /*t0*/, AmiVector const& /*x0*/, AmiVector const& /*dx0*/
 ) const {
     /* We need to set the steadystate rhs function. SUndials doesn't have this
        in its public api, so we have to change it in the solver memory,
@@ -157,7 +157,7 @@ void IDASolver::sensInit1(AmiVectorArray const& sx0, AmiVectorArray const& sdx0)
 }
 
 void IDASolver::binit(
-    int const which, const realtype tf, AmiVector const& xB0,
+    int const which, realtype const tf, AmiVector const& xB0,
     AmiVector const& dxB0
 ) const {
     int status;
@@ -263,13 +263,13 @@ void IDASolver::allocateSolver() const {
         );
 }
 
-void IDASolver::setSStolerances(const realtype rtol, const realtype atol)
+void IDASolver::setSStolerances(realtype const rtol, realtype const atol)
     const {
     int status = IDASStolerances(solver_memory_.get(), rtol, atol);
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDASStolerances");
 }
-void IDASolver::setSensSStolerances(const realtype rtol, realtype const* atol)
+void IDASolver::setSensSStolerances(realtype const rtol, realtype const* atol)
     const {
     int status = IDASensSStolerances(
         solver_memory_.get(), rtol, const_cast<realtype*>(atol)
@@ -374,11 +374,11 @@ void IDASolver::resetState(
     ida_mem->ida_kk = 0;
 }
 
-void IDASolver::reInitPostProcessF(const realtype tnext) const {
+void IDASolver::reInitPostProcessF(realtype const tnext) const {
     reInitPostProcess(solver_memory_.get(), &t_, &x_, &dx_, tnext);
 }
 
-void IDASolver::reInitPostProcessB(const realtype tnext) const {
+void IDASolver::reInitPostProcessB(realtype const tnext) const {
     realtype tBret;
     auto ida_mem = static_cast<IDAMem>(solver_memory_.get());
     auto idaadj_mem = ida_mem->ida_adj_mem;
@@ -445,7 +445,7 @@ void IDASolver::reInitPostProcess(
 }
 
 void IDASolver::reInit(
-    const realtype t0, AmiVector const& yy0, AmiVector const& yp0
+    realtype const t0, AmiVector const& yy0, AmiVector const& yp0
 ) const {
 
     auto ida_mem = static_cast<IDAMem>(solver_memory_.get());
@@ -492,7 +492,7 @@ void IDASolver::sensToggleOff() const {
 }
 
 void IDASolver::reInitB(
-    int const which, const realtype tB0, AmiVector const& yyB0,
+    int const which, realtype const tB0, AmiVector const& yyB0,
     AmiVector const& ypB0
 ) const {
 
@@ -526,7 +526,7 @@ void IDASolver::setSensParams(
         throw IDAException(status, "IDASetSensParams");
 }
 
-void IDASolver::getDky(const realtype t, int const k) const {
+void IDASolver::getDky(realtype const t, int const k) const {
     int status = IDAGetDky(solver_memory_.get(), t, k, dky_.getNVector());
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDAGetDky");
@@ -540,7 +540,7 @@ void IDASolver::getSens() const {
         throw IDAException(status, "IDAGetSens");
 }
 
-void IDASolver::getSensDky(const realtype t, int const k) const {
+void IDASolver::getSensDky(realtype const t, int const k) const {
     int status
         = IDAGetSensDky(solver_memory_.get(), t, k, sx_.getNVectorArray());
     if (status != IDA_SUCCESS)
@@ -557,7 +557,7 @@ void IDASolver::getB(int const which) const {
         throw IDAException(status, "IDAGetB");
 }
 
-void IDASolver::getDkyB(const realtype t, int k, int const which) const {
+void IDASolver::getDkyB(realtype const t, int k, int const which) const {
     int status = IDAGetDky(
         IDAGetAdjIDABmem(solver_memory_.get(), which), t, k, dky_.getNVector()
     );
@@ -578,7 +578,7 @@ void IDASolver::getQuad(realtype& t) const {
         throw IDAException(status, "IDAGetQuad");
 }
 
-void IDASolver::getQuadDkyB(const realtype t, int k, int const which) const {
+void IDASolver::getQuadDkyB(realtype const t, int k, int const which) const {
     int status = IDAGetQuadDky(
         IDAGetAdjIDABmem(solver_memory_.get(), which), t, k, xQB_.getNVector()
     );
@@ -586,7 +586,7 @@ void IDASolver::getQuadDkyB(const realtype t, int k, int const which) const {
         throw IDAException(status, "IDAGetB");
 }
 
-void IDASolver::getQuadDky(const realtype t, int const k) const {
+void IDASolver::getQuadDky(realtype const t, int const k) const {
     int status = IDAGetQuadDky(solver_memory_.get(), t, k, xQ_.getNVector());
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDAGetQuadDky");
@@ -639,7 +639,7 @@ void IDASolver::allocateSolverB(int* which) const {
 }
 
 void IDASolver::setSStolerancesB(
-    int const which, const realtype relTolB, const realtype absTolB
+    int const which, realtype const relTolB, realtype const absTolB
 ) const {
     int status
         = IDASStolerancesB(solver_memory_.get(), which, relTolB, absTolB);
@@ -648,7 +648,7 @@ void IDASolver::setSStolerancesB(
 }
 
 void IDASolver::quadSStolerancesB(
-    int const which, const realtype reltolQB, const realtype abstolQB
+    int const which, realtype const reltolQB, realtype const abstolQB
 ) const {
     int status
         = IDAQuadSStolerancesB(solver_memory_.get(), which, reltolQB, abstolQB);
@@ -657,14 +657,14 @@ void IDASolver::quadSStolerancesB(
 }
 
 void IDASolver::quadSStolerances(
-    const realtype reltolQB, const realtype abstolQB
+    realtype const reltolQB, realtype const abstolQB
 ) const {
     int status = IDAQuadSStolerances(solver_memory_.get(), reltolQB, abstolQB);
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDAQuadSStolerances");
 }
 
-int IDASolver::solve(const realtype tout, int const itask) const {
+int IDASolver::solve(realtype const tout, int const itask) const {
     if (force_reinit_postprocess_F_)
         reInitPostProcessF(tout);
     int status = IDASolve(
@@ -677,7 +677,7 @@ int IDASolver::solve(const realtype tout, int const itask) const {
     return status;
 }
 
-int IDASolver::solveF(const realtype tout, int const itask, int* ncheckPtr)
+int IDASolver::solveF(realtype const tout, int const itask, int* ncheckPtr)
     const {
     if (force_reinit_postprocess_F_)
         reInitPostProcessF(tout);
@@ -691,7 +691,7 @@ int IDASolver::solveF(const realtype tout, int const itask, int* ncheckPtr)
     return status;
 }
 
-void IDASolver::solveB(const realtype tBout, int const itaskB) const {
+void IDASolver::solveB(realtype const tBout, int const itaskB) const {
     if (force_reinit_postprocess_B_)
         reInitPostProcessB(tBout);
     int status = IDASolveB(solver_memory_.get(), tBout, itaskB);
@@ -768,7 +768,7 @@ void IDASolver::calcIC(realtype tout1) const {
         throw IDAException(status, "IDACalcIC");
 }
 
-void IDASolver::calcICB(int const which, const realtype tout1) const {
+void IDASolver::calcICB(int const which, realtype const tout1) const {
     int status = IDACalcICB(
         solver_memory_.get(), which, tout1, xB_.getNVector(), dxB_.getNVector()
     );
@@ -776,7 +776,7 @@ void IDASolver::calcICB(int const which, const realtype tout1) const {
         throw IDAException(status, "IDACalcICB");
 }
 
-void IDASolver::setStopTime(const realtype tstop) const {
+void IDASolver::setStopTime(realtype const tstop) const {
     int status = IDASetStopTime(solver_memory_.get(), tstop);
     if (status != IDA_SUCCESS)
         throw IDAException(status, "IDASetStopTime");
diff --git a/src/splinefunctions.cpp b/src/splinefunctions.cpp
index 8c888b450a..f5d2599ea7 100644
--- a/src/splinefunctions.cpp
+++ b/src/splinefunctions.cpp
@@ -54,18 +54,18 @@ AbstractSpline::AbstractSpline(
     }
 }
 
-realtype AbstractSpline::get_value(const realtype t) const {
+realtype AbstractSpline::get_value(realtype const t) const {
     auto y = get_value_scaled(t);
     return logarithmic_parametrization_ ? std::exp(y) : y;
 }
 
-realtype AbstractSpline::get_sensitivity(const realtype t, int const ip) const {
+realtype AbstractSpline::get_sensitivity(realtype const t, int const ip) const {
     auto s = get_sensitivity_scaled(t, ip);
     return logarithmic_parametrization_ ? s * get_value(t) : s;
 }
 
 realtype AbstractSpline::get_sensitivity(
-    const realtype t, int const ip, const realtype value
+    realtype const t, int const ip, realtype const value
 ) const {
     auto s = get_sensitivity_scaled(t, ip);
     return logarithmic_parametrization_ ? s * value : s;
@@ -387,12 +387,12 @@ void HermiteSpline::compute_coefficients_extrapolation() {
 #ifdef DVALUESDP
 #error "Preprocessor macro DVALUESDP already defined?!"
 #else
-#define DVALUESDP(i_node) dvaluesdp[node_offset + (i_node)*nplist]
+#define DVALUESDP(i_node) dvaluesdp[node_offset + (i_node) * nplist]
 #endif
 #ifdef DSLOPESDP
 #error "Preprocessor macro DSLOPESDP already defined?!"
 #else
-#define DSLOPESDP(i_node) dslopesdp[node_offset + (i_node)*nplist]
+#define DSLOPESDP(i_node) dslopesdp[node_offset + (i_node) * nplist]
 #endif
 
 void HermiteSpline::compute_coefficients_sensi(
@@ -821,7 +821,7 @@ void HermiteSpline::compute_final_sensitivity(
     set_final_sensitivity_scaled(finalSensitivity);
 }
 
-realtype HermiteSpline::get_value_scaled(const realtype t) const {
+realtype HermiteSpline::get_value_scaled(realtype const t) const {
     /* Is this a steady state computation? */
     if (std::isinf(t))
         return get_final_value_scaled();
@@ -922,7 +922,7 @@ realtype HermiteSpline::get_value_scaled(const realtype t) const {
 }
 
 realtype
-HermiteSpline::get_sensitivity_scaled(const realtype t, int const ip) const {
+HermiteSpline::get_sensitivity_scaled(realtype const t, int const ip) const {
     /* Is this a steady state computation? */
     if (std::isinf(t))
         return get_final_sensitivity_scaled(ip);
diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp
index a86574cd82..c415c26942 100644
--- a/src/sundials_matrix_wrapper.cpp
+++ b/src/sundials_matrix_wrapper.cpp
@@ -214,7 +214,7 @@ void SUNMatrixWrapper::scale(realtype a) {
 }
 
 void SUNMatrixWrapper::multiply(
-    N_Vector c, const_N_Vector b, const realtype alpha
+    N_Vector c, const_N_Vector b, realtype const alpha
 ) const {
     multiply(
         gsl::make_span<realtype>(NV_DATA_S(c), NV_LENGTH_S(c)),
@@ -233,7 +233,7 @@ inline static void check_csc(SUNMatrixWrapper const* /*mat*/) {}
 #endif
 
 void SUNMatrixWrapper::multiply(
-    gsl::span<realtype> c, gsl::span<realtype const> b, const realtype alpha
+    gsl::span<realtype> c, gsl::span<realtype const> b, realtype const alpha
 ) const {
 
     if (!matrix_)
@@ -522,8 +522,8 @@ void SUNMatrixWrapper::sparse_sum(std::vector<SUNMatrixWrapper> const& mats) {
 }
 
 sunindextype SUNMatrixWrapper::scatter(
-    const sunindextype acol, const realtype beta, sunindextype* w,
-    gsl::span<realtype> x, const sunindextype mark, SUNMatrixWrapper* C,
+    sunindextype const acol, realtype const beta, sunindextype* w,
+    gsl::span<realtype> x, sunindextype const mark, SUNMatrixWrapper* C,
     sunindextype nnz
 ) const {
     if (!matrix_)
@@ -576,7 +576,7 @@ static void cumsum(gsl::span<sunindextype> p, std::vector<sunindextype>& c) {
 }
 
 void SUNMatrixWrapper::transpose(
-    SUNMatrixWrapper& C, const realtype alpha, sunindextype blocksize
+    SUNMatrixWrapper& C, realtype const alpha, sunindextype blocksize
 ) const {
     if (!matrix_ || !C.matrix_)
         return;
@@ -804,7 +804,7 @@ unravel_index(sunindextype i, SUNMatrix m) {
             ++col;
 
         // This can happen if indexvals / indexptrs haven't been set.
-        if(col == ncols)
+        if (col == ncols)
             return {-1, -1};
 
         gsl_EnsuresDebug(row >= 0);
diff --git a/swig/CMakeLists.txt b/swig/CMakeLists.txt
index 87d607d33a..73faccda7e 100644
--- a/swig/CMakeLists.txt
+++ b/swig/CMakeLists.txt
@@ -112,10 +112,9 @@ if(${SWIG_VERSION} VERSION_LESS 4.1.0)
     PROPERTY SWIG_COMPILE_OPTIONS -py3)
 endif()
 
-# NOTE: No public definitions of any dependency are forwarded to swig,
-# they are only used for compiling the swig-generated source file.
-# Any definition that are relevant for swig-code generation, need to be
-# forwarded manually.
+# NOTE: No public definitions of any dependency are forwarded to swig, they are
+# only used for compiling the swig-generated source file. Any definition that
+# are relevant for swig-code generation, need to be forwarded manually.
 target_link_libraries(_amici amici Python3::Python)
 if(WIN32)
   add_custom_command(
diff --git a/swig/solver_cvodes.i b/swig/solver_cvodes.i
index 52feb1f4dd..7e2c28a17e 100644
--- a/swig/solver_cvodes.i
+++ b/swig/solver_cvodes.i
@@ -9,5 +9,16 @@ using namespace amici;
 %newobject amici::CVodeSolver::clone;
 %feature("notabstract") amici::CVodeSolver;
 
+// Required with SWIG 4.2.0 https://github.com/AMICI-dev/AMICI/issues/2275
+%extend amici::CVodeSolver {
+    CVodeSolver() {
+        return new CVodeSolver();
+    }
+
+    CVodeSolver(Solver const& solver) {
+        return new CVodeSolver(dynamic_cast<CVodeSolver const&>(solver));
+    }
+}
+
 // Process symbols in header
 %include "amici/solver_cvodes.h"
diff --git a/swig/solver_idas.i b/swig/solver_idas.i
index b4a49a836f..7ec9f663bb 100644
--- a/swig/solver_idas.i
+++ b/swig/solver_idas.i
@@ -9,5 +9,16 @@ using namespace amici;
 %newobject amici::IDASolver::clone;
 %feature("notabstract") amici::IDASolver;
 
+// Required for SWIG 4.2.0 https://github.com/AMICI-dev/AMICI/issues/2275
+%extend amici::IDASolver {
+    IDASolver() {
+        return new IDASolver();
+    }
+
+    IDASolver(Solver const& solver) {
+        return new IDASolver(dynamic_cast<IDASolver const&>(solver));
+    }
+}
+
 // Process symbols in header
 %include "amici/solver_idas.h"
diff --git a/tests/cpp/unittests/testMisc.cpp b/tests/cpp/unittests/testMisc.cpp
index f18de96b79..1f464b3433 100644
--- a/tests/cpp/unittests/testMisc.cpp
+++ b/tests/cpp/unittests/testMisc.cpp
@@ -261,6 +261,13 @@ TEST(SolverIdasTest, DefaultConstructableAndNotLeaky)
     IDASolver solver;
 }
 
+TEST(SolverIdasTest, CopyCtor)
+{
+    IDASolver solver1;
+    IDASolver solver2(solver1);
+}
+
+
 
 class SolverTest : public ::testing::Test {
   protected:
diff --git a/tox.ini b/tox.ini
index a61ee8a6e2..9f060a7cf0 100644
--- a/tox.ini
+++ b/tox.ini
@@ -13,7 +13,6 @@ description =
     Build documentation
 deps =
     -r documentation/rtd_requirements.txt
-    -r documentation/rtd_requirements2.txt
 # don't install the package, this is already handled by `deps` above
 skip_install = true
 change_dir = documentation/
diff --git a/version.txt b/version.txt
index a67cebaf7f..59dad104b0 100644
--- a/version.txt
+++ b/version.txt
@@ -1 +1 @@
-0.21.1
+0.21.2