diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml
index 018f0925..b01f7d66 100644
--- a/.github/workflows/ci_tests.yml
+++ b/.github/workflows/ci_tests.yml
@@ -11,7 +11,7 @@ jobs:
     strategy:
       matrix:
         platform: [windows-latest, macos-latest, ubuntu-latest]
-        python-version: ["3.7", "3.10"]
+        python-version: ["3.8", "3.10"]
     runs-on: ${{ matrix.platform }}
 
     steps:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4b6f44d1..f9b25ec9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,22 @@
 
 ## 0.1 series
 
+### 0.1.28
+
+* Fixed validation for output parameters columns in the condition table
+  by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/161
+* Added Python support policy
+  by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/162
+* Fixed typehints and deprecation warning
+  by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/165
+* Fixed SBML validation
+  by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/168
+* Fixed deprecation warning from `get_model_for_condition`
+  by @dweindl in https://github.com/PEtab-dev/libpetab-python/pull/169
+
+**Full Changelog**:
+https://github.com/PEtab-dev/libpetab-python/compare/v0.1.27...v0.1.28
+
 ### 0.1.27
 
 Features:
diff --git a/README.md b/README.md
index 2880878c..58fdc2bd 100644
--- a/README.md
+++ b/README.md
@@ -22,7 +22,8 @@ and the easiest way to install it is running
 
     pip3 install petab
     
-It will require Python>=3.7.1 to run.
+It will require Python>=3.8 to run. (We are following the
+[numpy Python support policy](https://numpy.org/neps/nep-0029-deprecation_policy.html)).
 
 Development versions of the PEtab library can be installed using
 
diff --git a/petab/conditions.py b/petab/conditions.py
index 59b07270..3e206463 100644
--- a/petab/conditions.py
+++ b/petab/conditions.py
@@ -93,9 +93,8 @@ def get_parametric_overrides(condition_df: pd.DataFrame) -> List[str]:
     Returns:
         List of parameter IDs that are mapped in a condition-specific way
     """
-    constant_parameters = list(
-        set(condition_df.columns.values.tolist()) - {CONDITION_ID,
-                                                     CONDITION_NAME})
+    constant_parameters = (set(condition_df.columns.values.tolist())
+                           - {CONDITION_ID, CONDITION_NAME})
     result = []
 
     for column in constant_parameters:
@@ -104,7 +103,5 @@ def get_parametric_overrides(condition_df: pd.DataFrame) -> List[str]:
 
         floatified = condition_df.loc[:, column].apply(core.to_float_if_float)
 
-        for x in floatified:
-            if not isinstance(x, float):
-                result.append(x)
+        result.extend(x for x in floatified if not isinstance(x, float))
     return result
diff --git a/petab/core.py b/petab/core.py
index 30256ca6..19f24faf 100644
--- a/petab/core.py
+++ b/petab/core.py
@@ -20,7 +20,7 @@
            'create_combine_archive', 'unique_preserve_order']
 
 
-def get_simulation_df(simulation_file: str) -> pd.DataFrame:
+def get_simulation_df(simulation_file: Union[str, Path]) -> pd.DataFrame:
     """Read PEtab simulation table
 
     Arguments:
@@ -33,7 +33,7 @@ def get_simulation_df(simulation_file: str) -> pd.DataFrame:
                        float_precision='round_trip')
 
 
-def write_simulation_df(df: pd.DataFrame, filename: str) -> None:
+def write_simulation_df(df: pd.DataFrame, filename: Union[str, Path]) -> None:
     """Write PEtab simulation table
 
     Arguments:
@@ -91,7 +91,8 @@ def get_notnull_columns(df: pd.DataFrame, candidates: Iterable):
 
 
 def flatten_timepoint_specific_output_overrides(
-        petab_problem: 'petab.problem.Problem') -> None:
+        petab_problem: 'petab.problem.Problem',
+) -> None:
     """Flatten timepoint-specific output parameter overrides.
 
     If the PEtab problem definition has timepoint-specific
diff --git a/petab/lint.py b/petab/lint.py
index d05d02de..50a4781b 100644
--- a/petab/lint.py
+++ b/petab/lint.py
@@ -60,9 +60,7 @@ def _check_df(df: pd.DataFrame, req_cols: Iterable, name: str) -> None:
     Raises:
           AssertionError: if a column is missing
     """
-    cols_set = df.columns.values
-    missing_cols = set(req_cols) - set(cols_set)
-    if missing_cols:
+    if missing_cols := set(req_cols) - set(df.columns.values):
         raise AssertionError(
             f"DataFrame {name} requires the columns {missing_cols}.")
 
@@ -85,12 +83,16 @@ def assert_no_leading_trailing_whitespace(
 
 
 def check_condition_df(
-        df: pd.DataFrame, model: Optional[Model] = None) -> None:
+        df: pd.DataFrame,
+        model: Optional[Model] = None,
+        observable_df: Optional[pd.DataFrame] = None
+) -> None:
     """Run sanity checks on PEtab condition table
 
     Arguments:
         df: PEtab condition DataFrame
         model: Model for additional checking of parameter IDs
+        observable_df: PEtab observables DataFrame
 
     Raises:
         AssertionError: in case of problems
@@ -101,7 +103,7 @@ def check_condition_df(
     _check_df(df, req_cols, "condition")
 
     # Check for correct index
-    if not df.index.name == CONDITION_ID:
+    if df.index.name != CONDITION_ID:
         raise AssertionError(
             f"Condition table has wrong index {df.index.name}."
             f"expected {CONDITION_ID}.")
@@ -119,6 +121,9 @@ def check_condition_df(
 
     if model is not None:
         allowed_cols = set(model.get_valid_ids_for_condition_table())
+        if observable_df is not None:
+            allowed_cols |= set(petab.get_output_parameters(
+                model=model, observable_df=observable_df))
         for column_name in df.columns:
             if column_name != CONDITION_NAME \
                     and column_name not in allowed_cols:
@@ -154,14 +159,9 @@ def check_measurement_df(df: pd.DataFrame,
                 df[column_name].values, column_name)
 
     if observable_df is not None:
-        # Check all observables are defined
-        observables_defined = set(observable_df.index.values)
-        observables_used = set(df[OBSERVABLE_ID])
-        observables_undefined = observables_used - observables_defined
-        if observables_undefined:
-            raise ValueError(f"Observables {observables_undefined} used in "
-                             "measurement table but not defined in "
-                             "observables table.")
+        assert_measured_observables_defined(df, observable_df)
+        measurements.assert_overrides_match_parameter_count(
+            df, observable_df)
 
         if OBSERVABLE_TRANSFORMATION in observable_df:
             # Check for positivity of measurements in case of
@@ -176,11 +176,6 @@ def check_measurement_df(df: pd.DataFrame,
                                      f'transformation {trafo} must be '
                                      f'positive, but {measurement} <= 0.')
 
-    if observable_df is not None:
-        assert_measured_observables_defined(df, observable_df)
-        measurements.assert_overrides_match_parameter_count(
-            df, observable_df)
-
     assert_measurements_not_null(df)
     assert_measurements_numeric(df)
 
@@ -206,7 +201,7 @@ def check_parameter_df(
 
     _check_df(df, PARAMETER_DF_REQUIRED_COLS[1:], "parameter")
 
-    if not df.index.name == PARAMETER_ID:
+    if df.index.name != PARAMETER_ID:
         raise AssertionError(
             f"Parameter table has wrong index {df.index.name}."
             f"expected {PARAMETER_ID}.")
@@ -232,10 +227,11 @@ def check_parameter_df(
                                  f"but column {NOMINAL_VALUE} is missing.")
         try:
             df.loc[non_estimated_par_ids, NOMINAL_VALUE].apply(float)
-        except ValueError:
-            raise AssertionError("Expected numeric values for "
-                                 f"`{NOMINAL_VALUE}` in parameter table for "
-                                 "all non-estimated parameters.")
+        except ValueError as e:
+            raise AssertionError(
+                f"Expected numeric values for `{NOMINAL_VALUE}` in parameter "
+                "table for all non-estimated parameters."
+            ) from e
 
     assert_parameter_id_is_string(df)
     assert_parameter_scale_is_valid(df)
@@ -285,8 +281,9 @@ def check_observable_df(observable_df: pd.DataFrame) -> None:
         try:
             sp.sympify(obs)
         except sp.SympifyError as e:
-            raise AssertionError(f"Cannot parse expression '{obs}' "
-                                 f"for observable {row.Index}: {e}")
+            raise AssertionError(
+                f"Cannot parse expression '{obs}' "
+                f"for observable {row.Index}: {e}") from e
 
         noise = getattr(row, NOISE_FORMULA)
         try:
@@ -297,9 +294,10 @@ def check_observable_df(observable_df: pd.DataFrame) -> None:
                 raise AssertionError(f"No or non-finite {NOISE_FORMULA} "
                                      f"given for observable {row.Index}.")
         except sp.SympifyError as e:
-            raise AssertionError(f"Cannot parse expression '{noise}' "
-                                 f"for noise model for observable "
-                                 f"{row.Index}: {e}")
+            raise AssertionError(
+                f"Cannot parse expression '{noise}' "
+                f"for noise model for observable " f"{row.Index}: {e}"
+            ) from e
 
 
 def assert_all_parameters_present_in_parameter_df(
@@ -346,7 +344,8 @@ def assert_all_parameters_present_in_parameter_df(
 
 def assert_measured_observables_defined(
         measurement_df: pd.DataFrame,
-        observable_df: pd.DataFrame) -> None:
+        observable_df: pd.DataFrame
+) -> None:
     """Check if all observables in the measurement table have been defined in
     the observable table
 
@@ -360,12 +359,11 @@ def assert_measured_observables_defined(
 
     used_observables = set(measurement_df[OBSERVABLE_ID].values)
     defined_observables = set(observable_df.index.values)
-    undefined_observables = used_observables - defined_observables
-
-    if undefined_observables:
+    if undefined_observables := (used_observables - defined_observables):
         raise AssertionError(
-            "Undefined observables in measurement file: "
-            f"{undefined_observables}.")
+            f"Observables {undefined_observables} used in "
+            "measurement table but not defined in observables table."
+        )
 
 
 def condition_table_is_parameter_free(condition_df: pd.DataFrame) -> bool:
@@ -540,9 +538,10 @@ def assert_parameter_prior_parameters_are_valid(
                 pars = tuple(
                     float(val) for val in pars_str.split(PARAMETER_SEPARATOR)
                 )
-            except ValueError:
+            except ValueError as e:
                 raise AssertionError(
-                    f"Could not parse prior parameters '{pars_str}'.")
+                    f"Could not parse prior parameters '{pars_str}'.") from e
+
             # all distributions take 2 parameters
             if len(pars) != 2:
                 raise AssertionError(
@@ -795,7 +794,8 @@ def lint_problem(problem: 'petab.Problem') -> bool:
     if problem.condition_df is not None:
         logger.info("Checking condition table...")
         try:
-            check_condition_df(problem.condition_df, problem.model)
+            check_condition_df(problem.condition_df, problem.model,
+                               problem.observable_df)
         except AssertionError as e:
             logger.error(e)
             errors_occurred = True
@@ -925,9 +925,7 @@ def assert_measurement_conditions_present_in_condition_table(
         used_conditions |= \
             set(measurement_df[PREEQUILIBRATION_CONDITION_ID].dropna().values)
     available_conditions = set(condition_df.index.values)
-    missing_conditions = used_conditions - available_conditions
-
-    if missing_conditions:
+    if missing_conditions := (used_conditions - available_conditions):
         raise AssertionError("Measurement table references conditions that "
                              "are not specified in the condition table: "
                              + str(missing_conditions))
diff --git a/petab/measurements.py b/petab/measurements.py
index 72cc2ce0..1ddac966 100644
--- a/petab/measurements.py
+++ b/petab/measurements.py
@@ -2,6 +2,7 @@
 # noqa: F405
 
 import itertools
+import math
 import numbers
 from pathlib import Path
 from typing import Dict, List, Union
@@ -201,7 +202,7 @@ def create_measurement_df() -> pd.DataFrame:
         Created DataFrame
     """
 
-    df = pd.DataFrame(data={
+    return pd.DataFrame(data={
         OBSERVABLE_ID: [],
         PREEQUILIBRATION_CONDITION_ID: [],
         SIMULATION_CONDITION_ID: [],
@@ -213,8 +214,6 @@ def create_measurement_df() -> pd.DataFrame:
         REPLICATE_ID: []
     })
 
-    return df
-
 
 def measurements_have_replicates(measurement_df: pd.DataFrame) -> bool:
     """Tests whether the measurements come with replicates
@@ -235,7 +234,8 @@ def measurements_have_replicates(measurement_df: pd.DataFrame) -> bool:
 
 def assert_overrides_match_parameter_count(
         measurement_df: pd.DataFrame,
-        observable_df: pd.DataFrame) -> None:
+        observable_df: pd.DataFrame
+) -> None:
     """Ensure that number of parameters in the observable definition matches
     the number of overrides in ``measurement_df``
 
@@ -243,7 +243,6 @@ def assert_overrides_match_parameter_count(
         measurement_df: PEtab measurement table
         observable_df: PEtab observable table
     """
-
     # sympify only once and save number of parameters
     observable_parameters_count = {
         obs_id: len(observables.get_formula_placeholders(
@@ -260,10 +259,11 @@ def assert_overrides_match_parameter_count(
         # check observable parameters
         try:
             expected = observable_parameters_count[row[OBSERVABLE_ID]]
-        except KeyError:
+        except KeyError as e:
             raise ValueError(
                 f"Observable {row[OBSERVABLE_ID]} used in measurement table "
-                f"is not defined.")
+                f"is not defined.") from e
+
         actual = len(split_parameter_replacement_list(
             row.get(OBSERVABLE_PARAMETERS, None)))
         # No overrides are also allowed
@@ -289,7 +289,7 @@ def assert_overrides_match_parameter_count(
         except KeyError:
             # no overrides defined, but a numerical sigma can be provided
             # anyways
-            if not len(replacements) == 1 \
+            if len(replacements) != 1 \
                     or not isinstance(replacements[0], numbers.Number):
                 raise AssertionError(
                     f'No placeholders have been specified in the noise model '
diff --git a/petab/models/sbml_model.py b/petab/models/sbml_model.py
index 3dcabe2f..5f79c3c3 100644
--- a/petab/models/sbml_model.py
+++ b/petab/models/sbml_model.py
@@ -9,7 +9,7 @@
 from . import MODEL_TYPE_SBML
 from .model import Model
 from ..sbml import (get_sbml_model, is_sbml_consistent, load_sbml_from_string,
-                    log_sbml_errors, write_sbml)
+                    write_sbml)
 
 
 class SbmlModel(Model):
@@ -109,9 +109,7 @@ def symbol_allowed_in_observable_formula(self, id_: str) -> bool:
         return self.sbml_model.getElementBySId(id_) or id_ == 'time'
 
     def is_valid(self) -> bool:
-        valid = is_sbml_consistent(self.sbml_model.getSBMLDocument())
-        log_sbml_errors(self.sbml_model.getSBMLDocument())
-        return valid
+        return is_sbml_consistent(self.sbml_model.getSBMLDocument())
 
     def is_state_variable(self, id_: str) -> bool:
         return (self.sbml_model.getSpecies(id_) is not None
diff --git a/petab/observables.py b/petab/observables.py
index 05fa142f..b1ecb986 100644
--- a/petab/observables.py
+++ b/petab/observables.py
@@ -3,7 +3,7 @@
 import re
 from collections import OrderedDict
 from pathlib import Path
-from typing import List, Union
+from typing import List, Union, Literal
 
 import pandas as pd
 import sympy as sp
@@ -105,8 +105,11 @@ def get_output_parameters(
     return list(output_parameters.keys())
 
 
-def get_formula_placeholders(formula_string: str, observable_id: str,
-                             override_type: str) -> List[str]:
+def get_formula_placeholders(
+        formula_string: str,
+        observable_id: str,
+        override_type: Literal['observable', 'noise'],
+) -> List[str]:
     """
     Get placeholder variables in noise or observable definition for the
     given observable ID.
@@ -114,8 +117,8 @@ def get_formula_placeholders(formula_string: str, observable_id: str,
     Arguments:
         formula_string: observable formula
         observable_id: ID of current observable
-        override_type: 'observable' or 'noise', depending on whether `formula`
-            is for observable or for noise model
+        override_type: ``'observable'`` or ``'noise'``, depending on whether
+            ``formula`` is for observable or for noise model
 
     Returns:
         List of placeholder parameter IDs in the order expected in the
@@ -192,6 +195,4 @@ def create_observable_df() -> pd.DataFrame:
         Created DataFrame
     """
 
-    df = pd.DataFrame(data={col: [] for col in OBSERVABLE_DF_COLS})
-
-    return df
+    return pd.DataFrame(data={col: [] for col in OBSERVABLE_DF_COLS})
diff --git a/petab/parameter_mapping.py b/petab/parameter_mapping.py
index 6a309011..01f8760a 100644
--- a/petab/parameter_mapping.py
+++ b/petab/parameter_mapping.py
@@ -6,7 +6,7 @@
 import os
 import re
 import warnings
-from typing import Any, Dict, Iterable, List, Optional, Tuple, Union
+from typing import Any, Dict, Iterable, List, Optional, Tuple, Union, Literal
 
 import libsbml
 import numpy as np
@@ -362,7 +362,8 @@ def _output_parameters_to_nan(mapping: ParMappingDict) -> None:
 
 def _apply_output_parameter_overrides(
         mapping: ParMappingDict,
-        cur_measurement_df: pd.DataFrame) -> None:
+        cur_measurement_df: pd.DataFrame
+) -> None:
     """
     Apply output parameter overrides to the parameter mapping dict for a given
     condition as defined in the measurement table (``observableParameter``,
@@ -390,8 +391,9 @@ def _apply_output_parameter_overrides(
 def _apply_overrides_for_observable(
         mapping: ParMappingDict,
         observable_id: str,
-        override_type: str,
-        overrides: List[str]) -> None:
+        override_type: Literal['observable', 'noise'],
+        overrides: List[str],
+) -> None:
     """
     Apply parameter-overrides for observables and noises to mapping
     matrix.
@@ -399,7 +401,7 @@ def _apply_overrides_for_observable(
     Arguments:
         mapping: mapping dict to which to apply overrides
         observable_id: observable ID
-        override_type: 'observable' or 'noise'
+        override_type: ``'observable'`` or ``'noise'``
         overrides: list of overrides for noise or observable parameters
     """
     for i, override in enumerate(overrides):
@@ -407,11 +409,13 @@ def _apply_overrides_for_observable(
         mapping[overridee_id] = override
 
 
-def _apply_condition_parameters(par_mapping: ParMappingDict,
-                                scale_mapping: ScaleMappingDict,
-                                condition_id: str,
-                                condition_df: pd.DataFrame,
-                                model: Model) -> None:
+def _apply_condition_parameters(
+        par_mapping: ParMappingDict,
+        scale_mapping: ScaleMappingDict,
+        condition_id: str,
+        condition_df: pd.DataFrame,
+        model: Model,
+) -> None:
     """Replace parameter IDs in parameter mapping dictionary by condition
     table parameter values (in-place).
 
@@ -447,12 +451,13 @@ def _apply_condition_parameters(par_mapping: ParMappingDict,
         scale_mapping[overridee_id] = LIN
 
 
-def _apply_parameter_table(par_mapping: ParMappingDict,
-                           scale_mapping: ScaleMappingDict,
-                           parameter_df: Optional[pd.DataFrame] = None,
-                           scaled_parameters: bool = False,
-                           fill_fixed_parameters: bool = True,
-                           ) -> None:
+def _apply_parameter_table(
+        par_mapping: ParMappingDict,
+        scale_mapping: ScaleMappingDict,
+        parameter_df: Optional[pd.DataFrame] = None,
+        scaled_parameters: bool = False,
+        fill_fixed_parameters: bool = True,
+) -> None:
     """Replace parameters from parameter table in mapping list for a given
     condition and set the corresponding scale.
 
@@ -536,15 +541,17 @@ def _perform_mapping_checks(
             "Timepoint-specific parameter overrides currently unsupported.")
 
 
-def handle_missing_overrides(mapping_par_opt_to_par_sim: ParMappingDict,
-                             warn: bool = True,
-                             condition_id: str = None) -> None:
+def handle_missing_overrides(
+        mapping_par_opt_to_par_sim: ParMappingDict,
+        warn: bool = True,
+        condition_id: str = None,
+) -> None:
     """
     Find all observable parameters and noise parameters that were not mapped
     and set their mapping to np.nan.
 
-    Assumes that parameters matching "(noise|observable)Parameter[0-9]+_" were
-    all supposed to be overwritten.
+    Assumes that parameters matching the regular expression
+    ``(noise|observable)Parameter[0-9]+_`` were all supposed to be overwritten.
 
     Parameters:
         mapping_par_opt_to_par_sim:
@@ -578,14 +585,15 @@ def merge_preeq_and_sim_pars_condition(
         condition_map_sim: ParMappingDict,
         condition_scale_map_preeq: ScaleMappingDict,
         condition_scale_map_sim: ScaleMappingDict,
-        condition: Any) -> None:
+        condition: Any,
+) -> None:
     """Merge preequilibration and simulation parameters and scales for a single
     condition while checking for compatibility.
 
     This function is meant for the case where we cannot have different
     parameters (and scales) for preequilibration and simulation. Therefore,
     merge both and ensure matching scales and parameters.
-    ``condition_map_sim`` and ``condition_scale_map_sim`` will ne modified in
+    ``condition_map_sim`` and ``condition_scale_map_sim`` will be modified in
     place.
 
     Arguments:
diff --git a/petab/parameters.py b/petab/parameters.py
index 45f6febf..e90c6f90 100644
--- a/petab/parameters.py
+++ b/petab/parameters.py
@@ -4,7 +4,9 @@
 import warnings
 from collections import OrderedDict
 from pathlib import Path
-from typing import Dict, Iterable, List, Set, Tuple, Union, Optional
+from typing import (
+    Dict, Iterable, List, Set, Tuple, Union, Optional, Literal, Sequence
+)
 
 import libsbml
 import numpy as np
@@ -19,6 +21,7 @@
            'get_optimization_parameters',
            'get_parameter_df',
            'get_priors_from_df',
+           'get_valid_parameters_for_parameter_table',
            'map_scale',
            'map_unscale',
            'normalize_parameter_df',
@@ -26,6 +29,8 @@
            'unscale',
            'write_parameter_df']
 
+PARAMETER_SCALE_ARGS = Literal['', 'lin', 'log', 'log10']
+
 
 def get_parameter_df(
         parameter_file: Union[str, Path, pd.DataFrame,
@@ -147,13 +152,13 @@ def create_parameter_df(
     matching the number of parameters
 
     Arguments:
-        sbml_model: SBML Model
-        model: PEtab model
+        sbml_model: SBML Model (deprecated, mutually exclusive with ``model``)
+        model: PEtab model (mutually exclusive with ``sbml_model``)
         condition_df: PEtab condition DataFrame
         observable_df: PEtab observable DataFrame
         measurement_df: PEtab measurement DataFrame
         include_optional: By default this only returns parameters that are
-            required to be present in the parameter table. If set to True,
+            required to be present in the parameter table. If set to ``True``,
             this returns all parameters that are allowed to be present in the
             parameter table (i.e. also including parameters specified in the
             model).
@@ -265,6 +270,13 @@ def append_overrides(overrides):
         if not model.has_entity_with_id(p):
             parameter_ids[p] = None
 
+    # remove parameters that occur in the condition table and are overridden
+    #  for ALL conditions
+    for p in condition_df.columns[~condition_df.isnull().any()]:
+        try:
+            del parameter_ids[p]
+        except KeyError:
+            pass
     return parameter_ids.keys()
 
 
@@ -272,7 +284,8 @@ def get_valid_parameters_for_parameter_table(
         model: Model,
         condition_df: pd.DataFrame,
         observable_df: pd.DataFrame,
-        measurement_df: pd.DataFrame) -> Set[str]:
+        measurement_df: pd.DataFrame,
+) -> Set[str]:
     """
     Get set of parameters which may be present inside the parameter table
 
@@ -337,13 +350,15 @@ def append_overrides(overrides):
     return parameter_ids.keys()
 
 
-def get_priors_from_df(parameter_df: pd.DataFrame,
-                       mode: str) -> List[Tuple]:
+def get_priors_from_df(
+        parameter_df: pd.DataFrame,
+        mode: Literal['initialization', 'objective'],
+) -> List[Tuple]:
     """Create list with information about the parameter priors
 
     Arguments:
         parameter_df: PEtab parameter table
-        mode: 'initialization' or 'objective'
+        mode: ``'initialization'`` or ``'objective'``
 
     Returns:
         List with prior information.
@@ -384,55 +399,59 @@ def get_priors_from_df(parameter_df: pd.DataFrame,
     return prior_list
 
 
-def scale(parameter: numbers.Number, scale_str: 'str') -> numbers.Number:
-    """Scale parameter according to `scale_str`.
+def scale(
+        parameter: numbers.Number,
+        scale_str: PARAMETER_SCALE_ARGS,
+) -> numbers.Number:
+    """Scale parameter according to ``scale_str``.
 
     Arguments:
         parameter:
             Parameter to be scaled.
         scale_str:
-            One of 'lin' (synonymous with ''), 'log', 'log10'.
+            One of ``'lin'`` (synonymous with ``''``), ``'log'``, ``'log10'``.
 
     Returns:
         The scaled parameter.
     """
-
     if scale_str == LIN or not scale_str:
         return parameter
     if scale_str == LOG:
         return np.log(parameter)
     if scale_str == LOG10:
         return np.log10(parameter)
-    raise ValueError("Invalid parameter scaling: " + scale_str)
+    raise ValueError(f"Invalid parameter scaling: {scale_str}")
 
 
-def unscale(parameter: numbers.Number, scale_str: 'str') -> numbers.Number:
-    """Unscale parameter according to `scale_str`.
+def unscale(
+        parameter: numbers.Number,
+        scale_str: PARAMETER_SCALE_ARGS,
+) -> numbers.Number:
+    """Unscale parameter according to ``scale_str``.
 
     Arguments:
         parameter:
             Parameter to be unscaled.
         scale_str:
-            One of 'lin' (synonymous with ''), 'log', 'log10'.
+            One of ``'lin'`` (synonymous with ``''``), ``'log'``, ``'log10'``.
 
     Returns:
         The unscaled parameter.
     """
-
     if scale_str == LIN or not scale_str:
         return parameter
     if scale_str == LOG:
         return np.exp(parameter)
     if scale_str == LOG10:
         return 10**parameter
-    raise ValueError("Invalid parameter scaling: " + scale_str)
+    raise ValueError(f"Invalid parameter scaling: {scale_str}")
 
 
 def map_scale(
-    parameters: Iterable[numbers.Number],
-    scale_strs: Union[Iterable[str], str]
+    parameters: Sequence[numbers.Number],
+    scale_strs: Union[Iterable[PARAMETER_SCALE_ARGS], PARAMETER_SCALE_ARGS],
 ) -> Iterable[numbers.Number]:
-    """Scale the parameters, i.e. as `scale()`, but for Iterables.
+    """Scale the parameters, i.e. as :func:`scale`, but for Sequences.
 
     Arguments:
         parameters:
@@ -449,10 +468,10 @@ def map_scale(
 
 
 def map_unscale(
-    parameters: Iterable[numbers.Number],
-    scale_strs: Union[Iterable[str], str]
+    parameters: Sequence[numbers.Number],
+    scale_strs: Union[Iterable[PARAMETER_SCALE_ARGS], PARAMETER_SCALE_ARGS],
 ) -> Iterable[numbers.Number]:
-    """Unscale the parameters, i.e. as `unscale()`, but for Iterables.
+    """Unscale the parameters, i.e. as :func:`unscale`, but for Sequences.
 
     Arguments:
         parameters:
diff --git a/petab/petablint.py b/petab/petablint.py
index 2c2d8f64..8be50aa0 100755
--- a/petab/petablint.py
+++ b/petab/petablint.py
@@ -32,7 +32,6 @@ def format(self, record):
 
 def parse_cli_args():
     """Parse command line arguments"""
-
     parser = argparse.ArgumentParser(
         description='Check if a set of files adheres to the PEtab format.')
 
diff --git a/petab/problem.py b/petab/problem.py
index 4e274615..73349f7b 100644
--- a/petab/problem.py
+++ b/petab/problem.py
@@ -71,8 +71,8 @@ def __init__(
         if any((sbml_model, sbml_document, sbml_reader),):
             warn("Passing `sbml_model`, `sbml_document`, or `sbml_reader` "
                  "to petab.Problem is deprecated and will be removed in a "
-                 "future version. Use `model=petab.models.SbmlModel(...)` "
-                 "instead.", DeprecationWarning, stacklevel=2)
+                 "future version. Use `model=petab.models.sbml_model."
+                 "SbmlModel(...)` instead.", DeprecationWarning, stacklevel=2)
             if model:
                 raise ValueError("Must only provide one of (`sbml_model`, "
                                  "`sbml_document`, `sbml_reader`) or `model`.")
@@ -129,6 +129,7 @@ def from_files(
             parameter_file: PEtab parameter table
             visualization_files: PEtab visualization tables
             observable_files: PEtab observables tables
+            model_id: PEtab ID of the model
             extensions_config: Information on the extensions used
         """
         warn("petab.Problem.from_files is deprecated and will be removed in a "
@@ -327,8 +328,8 @@ def to_files_generic(
             prefix_path:
                 Specify a prefix to all paths, to avoid specifying the
                 prefix for all paths individually. NB: the prefix is added to
-                paths before `relative_paths` is handled downstream in
-                `petab.yaml.create_problem_yaml`.
+                paths before ``relative_paths`` is handled downstream in
+                :func:`petab.yaml.create_problem_yaml`.
 
         Returns:
             The path to the PEtab problem YAML file.
@@ -362,18 +363,19 @@ def to_files_generic(
             return filenames['yaml_file']
         return str(prefix_path / filenames['yaml_file'])
 
-    def to_files(self,
-                 sbml_file: Union[None, str, Path] = None,
-                 condition_file: Union[None, str, Path] = None,
-                 measurement_file: Union[None, str, Path] = None,
-                 parameter_file: Union[None, str, Path] = None,
-                 visualization_file: Union[None, str, Path] = None,
-                 observable_file: Union[None, str, Path] = None,
-                 yaml_file: Union[None, str, Path] = None,
-                 prefix_path: Union[None, str, Path] = None,
-                 relative_paths: bool = True,
-                 model_file: Union[None, str, Path] = None,
-                 ) -> None:
+    def to_files(
+            self,
+            sbml_file: Union[None, str, Path] = None,
+            condition_file: Union[None, str, Path] = None,
+            measurement_file: Union[None, str, Path] = None,
+            parameter_file: Union[None, str, Path] = None,
+            visualization_file: Union[None, str, Path] = None,
+            observable_file: Union[None, str, Path] = None,
+            yaml_file: Union[None, str, Path] = None,
+            prefix_path: Union[None, str, Path] = None,
+            relative_paths: bool = True,
+            model_file: Union[None, str, Path] = None,
+    ) -> None:
         """
         Write PEtab tables to files for this problem
 
@@ -395,10 +397,10 @@ def to_files(self,
             prefix_path:
                 Specify a prefix to all paths, to avoid specifying the
                 prefix for all paths individually. NB: the prefix is added to
-                paths before `relative_paths` is handled.
+                paths before ``relative_paths`` is handled.
             relative_paths:
                 whether all paths in the YAML file should be
-                relative to the location of the YAML file. If `False`, then
+                relative to the location of the YAML file. If ``False``, then
                 paths are left unchanged.
 
         Raises:
@@ -420,9 +422,7 @@ def to_files(self,
             prefix_path = Path(prefix_path)
 
             def add_prefix(path0: Union[None, str, Path]) -> str:
-                if path0 is None:
-                    return path0
-                return str(prefix_path / path0)
+                return path0 if path0 is None else str(prefix_path / path0)
 
             model_file = add_prefix(model_file)
             condition_file = add_prefix(condition_file)
diff --git a/petab/sampling.py b/petab/sampling.py
index 9c226cd0..51096b43 100644
--- a/petab/sampling.py
+++ b/petab/sampling.py
@@ -11,18 +11,20 @@
 __all__ = ['sample_from_prior', 'sample_parameter_startpoints']
 
 
-def sample_from_prior(prior: Tuple[str, list, str, list],
-                      n_starts: int) -> np.array:
+def sample_from_prior(
+        prior: Tuple[str, list, str, list],
+        n_starts: int
+) -> np.array:
     """Creates samples for one parameter based on prior
 
     Arguments:
-        prior: A tuple as obtained from ``petab.parameter.get_priors_from_df``
+        prior: A tuple as obtained from
+            :func:`petab.parameter.get_priors_from_df`
         n_starts: Number of samples
 
     Returns:
         Array with sampled values
     """
-
     # unpack info
     p_type, p_params, scaling, bounds = prior
 
@@ -40,8 +42,7 @@ def scale(x):
 
     def clip_to_bounds(x: np.array):
         """Clip values in array x to bounds"""
-        x = np.maximum(np.minimum(scale(bounds[1]), x), scale(bounds[0]))
-        return x
+        return np.maximum(np.minimum(scale(bounds[1]), x), scale(bounds[0]))
 
     # define lambda functions for each parameter
     if p_type == UNIFORM:
@@ -83,15 +84,17 @@ def clip_to_bounds(x: np.array):
     return clip_to_bounds(sp)
 
 
-def sample_parameter_startpoints(parameter_df: pd.DataFrame,
-                                 n_starts: int = 100,
-                                 seed: int = None) -> np.array:
-    """Create numpy.array with starting points for an optimization
+def sample_parameter_startpoints(
+        parameter_df: pd.DataFrame,
+        n_starts: int = 100,
+        seed: int = None,
+) -> np.array:
+    """Create :class:`numpy.array` with starting points for an optimization
 
     Arguments:
         parameter_df: PEtab parameter DataFrame
         n_starts: Number of points to be sampled
-        seed: Random number generator seed (see numpy.random.seed)
+        seed: Random number generator seed (see :func:`numpy.random.seed`)
 
     Returns:
         Array of sampled starting points with dimensions
diff --git a/petab/sbml.py b/petab/sbml.py
index 47c64587..f12fec76 100644
--- a/petab/sbml.py
+++ b/petab/sbml.py
@@ -24,8 +24,10 @@
 ]
 
 
-def is_sbml_consistent(sbml_document: libsbml.SBMLDocument,
-                       check_units: bool = False) -> bool:
+def is_sbml_consistent(
+        sbml_document: libsbml.SBMLDocument,
+        check_units: bool = False,
+) -> bool:
     """Check for SBML validity / consistency
 
     Arguments:
@@ -49,31 +51,35 @@ def is_sbml_consistent(sbml_document: libsbml.SBMLDocument,
     return not has_problems
 
 
-def log_sbml_errors(sbml_document: libsbml.SBMLDocument,
-                    minimum_severity=libsbml.LIBSBML_SEV_WARNING) -> None:
+def log_sbml_errors(
+        sbml_document: libsbml.SBMLDocument,
+        minimum_severity=libsbml.LIBSBML_SEV_WARNING,
+) -> None:
     """Log libsbml errors
 
     Arguments:
         sbml_document: SBML document to check
-        minimum_severity: Minimum severity level to report (see libsbml)
+        minimum_severity: Minimum severity level to report (see libsbml
+            documentation)
     """
-
+    severity_to_log_level = {
+        libsbml.LIBSBML_SEV_INFO: logging.INFO,
+        libsbml.LIBSBML_SEV_WARNING: logging.WARNING,
+    }
     for error_idx in range(sbml_document.getNumErrors()):
         error = sbml_document.getError(error_idx)
-        if error.getSeverity() >= minimum_severity:
+        if (severity := error.getSeverity()) >= minimum_severity:
             category = error.getCategoryAsString()
-            severity = error.getSeverityAsString()
+            severity_str = error.getSeverityAsString()
             message = error.getMessage()
-            if severity == libsbml.LIBSBML_SEV_INFO:
-                logger.info(f'libSBML {severity} ({category}): {message}')
-            elif severity == libsbml.LIBSBML_SEV_WARNING:
-                logger.warning(f'libSBML {severity} ({category}): {message}')
-            else:
-                logger.error(f'libSBML {severity} ({category}): {message}')
+            logger.log(severity_to_log_level.get(severity, logging.ERROR),
+                       f'libSBML {severity_str} ({category}): {message}')
 
 
-def globalize_parameters(sbml_model: libsbml.Model,
-                         prepend_reaction_id: bool = False) -> None:
+def globalize_parameters(
+        sbml_model: libsbml.Model,
+        prepend_reaction_id: bool = False,
+) -> None:
     """Turn all local parameters into global parameters with the same
     properties
 
@@ -81,8 +87,8 @@ def globalize_parameters(sbml_model: libsbml.Model,
     function to convert them to global parameters. There may exist local
     parameters with identical IDs within different kinetic laws. This is not
     checked here. If in doubt that local parameter IDs are unique, enable
-    `prepend_reaction_id` to create global parameters named
-    ${reaction_id}_${local_parameter_id}.
+    ``prepend_reaction_id`` to create global parameters named
+    ``${reaction_id}_${local_parameter_id}``.
 
     Arguments:
         sbml_model:
@@ -91,7 +97,6 @@ def globalize_parameters(sbml_model: libsbml.Model,
             Prepend reaction id of local parameter when
             creating global parameters
     """
-
     warn("This function will be removed in future releases.",
          DeprecationWarning)
 
@@ -118,15 +123,16 @@ def globalize_parameters(sbml_model: libsbml.Model,
             law.removeParameter(lp.getId())
 
 
-def get_model_parameters(sbml_model: libsbml.Model, with_values=False
-                         ) -> Union[List[str], Dict[str, float]]:
+def get_model_parameters(
+        sbml_model: libsbml.Model, with_values=False
+) -> Union[List[str], Dict[str, float]]:
     """Return SBML model parameters which are not Rule targets
 
     Arguments:
         sbml_model: SBML model
         with_values:
-            If False, returns list of SBML model parameter IDs which
-            are not Rule targets. If True, returns a dictionary with those
+            If ``False``, returns list of SBML model parameter IDs which
+            are not Rule targets. If ``True``, returns a dictionary with those
             parameter IDs as keys and parameter values from the SBML model as
             values.
     """
@@ -197,9 +203,8 @@ def load_sbml_from_file(
     """Load SBML model from file
 
     :param sbml_file: Filename of the SBML file
-    :return: The SBML document, model and reader
+    :return: The SBML reader, document, model
     """
-
     sbml_reader = libsbml.SBMLReader()
     sbml_document = sbml_reader.readSBML(sbml_file)
     sbml_model = sbml_document.getModel()
@@ -227,6 +232,9 @@ def get_model_for_condition(
         in the resulting model.
     :return: The generated SBML document, and SBML model
     """
+    from .models.sbml_model import SbmlModel
+    assert isinstance(petab_problem.model, SbmlModel)
+
     condition_dict = {petab.SIMULATION_CONDITION_ID: sim_condition_id}
     if preeq_condition_id:
         condition_dict[petab.PREEQUILIBRATION_CONDITION_ID] = \
@@ -240,7 +248,7 @@ def get_model_for_condition(
             condition_id=sim_condition_id,
             is_preeq=False,
             cur_measurement_df=cur_measurement_df,
-            sbml_model=petab_problem.sbml_model,
+            model=petab_problem.model,
             condition_df=petab_problem.condition_df,
             parameter_df=petab_problem.parameter_df,
             warn_unmapped=True,
@@ -251,7 +259,7 @@ def get_model_for_condition(
             allow_timepoint_specific_numeric_noise_parameters=True,
         )
     # create a copy of the model
-    sbml_doc = petab_problem.sbml_model.getSBMLDocument().clone()
+    sbml_doc = petab_problem.model.sbml_model.getSBMLDocument().clone()
     sbml_model = sbml_doc.getModel()
 
     # fill in parameters
diff --git a/petab/simulate.py b/petab/simulate.py
index 5d2a065f..139170d5 100644
--- a/petab/simulate.py
+++ b/petab/simulate.py
@@ -17,7 +17,7 @@ class Simulator(abc.ABC):
     """Base class that specific simulators should inherit.
 
     Specific simulators should minimally implement the
-    `simulate_without_noise` method.
+    :meth:`petab.simulate.Simulator.simulate_without_noise` method.
     Example (AMICI): https://bit.ly/33SUSG4
 
     Attributes:
@@ -29,7 +29,7 @@ class Simulator(abc.ABC):
         rng:
             A NumPy random generator, used to sample from noise distributions.
         temporary_working_dir:
-            Whether `working_dir` is a temporary directory, which can be
+            Whether ``working_dir`` is a temporary directory, which can be
             deleted without significant consequence.
         working_dir:
             All simulator-specific output files will be saved here. This
@@ -74,14 +74,15 @@ def __init__(
     def remove_working_dir(self, force: bool = False, **kwargs) -> None:
         """Remove the simulator working directory, and all files within.
 
-        See the `__init__` method arguments.
+        See the :meth:`petab.simulate.Simulator.__init__` method arguments.
 
         Arguments:
             force:
-                If True, the working directory is removed regardless of
+                If ``True``, the working directory is removed regardless of
                 whether it is a temporary directory.
             **kwargs:
-                Additional keyword arguments are passed to `shutil.rmtree`.
+                Additional keyword arguments are passed to
+                :func:`shutil.rmtree`.
         """
         if force or self.temporary_working_dir:
             shutil.rmtree(self.working_dir, **kwargs)
@@ -102,11 +103,13 @@ def simulate_without_noise(self) -> pd.DataFrame:
 
         Returns:
             Simulated data, as a PEtab measurements table, which should be
-            equivalent to replacing all values in the `petab.C.MEASUREMENT`
-            column of the measurements table (of the PEtab problem supplied to
-            the `__init__` method), with simulated values.
+            equivalent to replacing all values in the
+            :const:`petab.C.MEASUREMENT` column of the measurements table (of
+            the PEtab problem supplied to the
+            :meth:`petab.simulate.Simulator.__init__` method), with
+            simulated values.
         """
-        raise NotImplementedError
+        raise NotImplementedError()
 
     def simulate(
             self,
@@ -122,7 +125,7 @@ def simulate(
                 A multiplier of the scale of the noise distribution.
             **kwargs:
                 Additional keyword arguments are passed to
-                `simulate_without_noise`.
+                :meth:`petab.simulate.Simulator.simulate_without_noise`.
 
         Returns:
             Simulated data, as a PEtab measurements table.
@@ -146,7 +149,8 @@ def add_noise(
             noise_scaling_factor:
                 A multiplier of the scale of the noise distribution.
             **kwargs:
-                Additional keyword arguments are passed to `sample_noise`.
+                Additional keyword arguments are passed to
+                :func:`sample_noise`.
 
         Returns:
             Simulated data with noise, as a PEtab measurements table.
@@ -181,7 +185,7 @@ def sample_noise(
     Arguments:
         petab_problem:
             The PEtab problem used to generate the simulated value.
-            Instance of `petab.Problem`.
+            Instance of :class:`petab.Problem`.
         measurement_row:
             The row in the PEtab problem measurement table that corresponds
             to the simulated value.
@@ -189,16 +193,15 @@ def sample_noise(
             A simulated value without noise.
         noise_formulas:
             Processed noise formulas from the PEtab observables table, in the
-            form output by the `petab.calculate.get_symbolic_noise_formulas`
-            method.
+            form output by :func:`petab.calculate.get_symbolic_noise_formulas`.
         rng:
             A NumPy random generator.
         noise_scaling_factor:
             A multiplier of the scale of the noise distribution.
         zero_bounded:
-            Return zero if the sign of the return value and `simulated_value`
+            Return zero if the sign of the return value and ``simulated_value``
             differ. Can be used to ensure non-negative and non-positive values,
-            if the sign of `simulated_value` should not change.
+            if the sign of ``simulated_value`` should not change.
 
     Returns:
         The sample from the PEtab noise distribution.
diff --git a/petab/version.py b/petab/version.py
index 7f84ffb2..e7c1d81b 100644
--- a/petab/version.py
+++ b/petab/version.py
@@ -1,2 +1,2 @@
 """PEtab library version"""
-__version__ = '0.1.27'
+__version__ = '0.1.28'
diff --git a/petab/visualize/plotter.py b/petab/visualize/plotter.py
index 2e0caa5e..4e2ca82d 100644
--- a/petab/visualize/plotter.py
+++ b/petab/visualize/plotter.py
@@ -33,8 +33,10 @@ def __init__(self, figure: Figure, data_provider: DataProvider):
         self.data_provider = data_provider
 
     @abstractmethod
-    def generate_figure(self, subplot_dir: Optional[str] = None
-                        ) -> Optional[Dict[str, plt.Subplot]]:
+    def generate_figure(
+            self,
+            subplot_dir: Optional[str] = None
+    ) -> Optional[Dict[str, plt.Subplot]]:
         pass
 
 
@@ -65,9 +67,12 @@ def _error_column_for_plot_type_data(plot_type_data: str) -> Optional[str]:
             return 'noise_model'
         return None
 
-    def generate_lineplot(self, ax: 'matplotlib.pyplot.Axes',
-                          dataplot: DataPlot,
-                          plotTypeData: str) -> None:
+    def generate_lineplot(
+            self,
+            ax: 'matplotlib.pyplot.Axes',
+            dataplot: DataPlot,
+            plotTypeData: str
+    ) -> None:
         """
         Generate lineplot.
 
@@ -82,7 +87,6 @@ def generate_lineplot(self, ax: 'matplotlib.pyplot.Axes',
         plotTypeData:
             Specifies how replicates should be handled.
         """
-
         simu_color = None
         measurements_to_plot, simulations_to_plot = \
             self.data_provider.get_data_to_plot(dataplot,
@@ -152,9 +156,12 @@ def generate_lineplot(self, ax: 'matplotlib.pyplot.Axes',
                 label=label_base + " simulation", color=simu_color
             )
 
-    def generate_barplot(self, ax: 'matplotlib.pyplot.Axes',
-                         dataplot: DataPlot,
-                         plotTypeData: str) -> None:
+    def generate_barplot(
+            self,
+            ax: 'matplotlib.pyplot.Axes',
+            dataplot: DataPlot,
+            plotTypeData: str
+    ) -> None:
         """
         Generate barplot.
 
@@ -200,9 +207,12 @@ def generate_barplot(self, ax: 'matplotlib.pyplot.Axes',
                    color='white', edgecolor=color, **bar_kwargs,
                    label='simulation')
 
-    def generate_scatterplot(self, ax: 'matplotlib.pyplot.Axes',
-                             dataplot: DataPlot,
-                             plotTypeData: str) -> None:
+    def generate_scatterplot(
+            self,
+            ax: 'matplotlib.pyplot.Axes',
+            dataplot: DataPlot,
+            plotTypeData: str
+    ) -> None:
         """
         Generate scatterplot.
 
@@ -215,7 +225,6 @@ def generate_scatterplot(self, ax: 'matplotlib.pyplot.Axes',
         plotTypeData:
             Specifies how replicates should be handled.
         """
-
         measurements_to_plot, simulations_to_plot = \
             self.data_provider.get_data_to_plot(dataplot,
                                                 plotTypeData == PROVIDED)
@@ -228,9 +237,11 @@ def generate_scatterplot(self, ax: 'matplotlib.pyplot.Axes',
                    label=getattr(dataplot, LEGEND_ENTRY))
         self._square_plot_equal_ranges(ax)
 
-    def generate_subplot(self,
-                         ax,
-                         subplot: Subplot) -> None:
+    def generate_subplot(
+            self,
+            ax: plt.Axes,
+            subplot: Subplot
+    ) -> None:
         """
         Generate subplot based on markup provided by subplot.
 
@@ -241,7 +252,6 @@ def generate_subplot(self,
         subplot:
             Subplot visualization settings.
         """
-
         # set yScale
         if subplot.yScale == LIN:
             ax.set_yscale("linear")
@@ -270,7 +280,6 @@ def generate_subplot(self,
             for data_plot in subplot.data_plots:
                 self.generate_scatterplot(ax, data_plot, subplot.plotTypeData)
         else:
-
             # set xScale
             if subplot.xScale == LIN:
                 ax.set_xscale("linear")
@@ -345,7 +354,6 @@ def generate_figure(
         None:
             In case subplots are saved to file.
         """
-
         if subplot_dir is None:
             # compute, how many rows and columns we need for the subplots
             num_row = int(np.round(np.sqrt(self.figure.num_subplots)))
@@ -361,7 +369,7 @@ def generate_figure(
             axes = dict(zip([plot.plotId for plot in self.figure.subplots],
                             axes.flat))
 
-        for idx, subplot in enumerate(self.figure.subplots):
+        for subplot in self.figure.subplots:
             if subplot_dir is not None:
                 fig, ax = plt.subplots(figsize=self.figure.size)
                 fig.set_tight_layout(True)
@@ -419,6 +427,8 @@ class SeabornPlotter(Plotter):
     def __init__(self, figure: Figure, data_provider: DataProvider):
         super().__init__(figure, data_provider)
 
-    def generate_figure(self, subplot_dir: Optional[str] = None
-                        ) -> Optional[Dict[str, plt.Subplot]]:
+    def generate_figure(
+            self,
+            subplot_dir: Optional[str] = None
+    ) -> Optional[Dict[str, plt.Subplot]]:
         pass
diff --git a/petab/visualize/plotting.py b/petab/visualize/plotting.py
index 45956b17..041b3042 100644
--- a/petab/visualize/plotting.py
+++ b/petab/visualize/plotting.py
@@ -2,7 +2,7 @@
 import warnings
 from numbers import Number, Real
 from pathlib import Path
-from typing import Dict, List, Optional, Tuple, Union
+from typing import Dict, List, Optional, Tuple, Union, Literal
 
 import numpy as np
 import pandas as pd
@@ -130,6 +130,9 @@ def from_df(cls, plot_spec: pd.DataFrame):
 
         return cls(vis_spec_dict)
 
+    def __repr__(self):
+        return f"{self.__class__.__name__}({self.__dict__})"
+
 
 class Subplot:
     """
@@ -480,9 +483,13 @@ def _get_independent_var_values(self, data_df: pd.DataFrame,
 
         return uni_condition_id, col_name_unique, conditions_
 
-    def get_data_series(self, data_df: pd.DataFrame, data_col: str,
-                        dataplot: DataPlot,
-                        provided_noise: bool) -> DataSeries:
+    def get_data_series(
+            self,
+            data_df: pd.DataFrame,
+            data_col: Literal['measurement', 'simulation'],
+            dataplot: DataPlot,
+            provided_noise: bool
+    ) -> DataSeries:
         """
         Get data to plot from measurement or simulation DataFrame.
 
@@ -499,10 +506,8 @@ def get_data_series(self, data_df: pd.DataFrame, data_col: str,
         -------
         Data to plot
         """
-
         uni_condition_id, col_name_unique, conditions_ = \
-            self._get_independent_var_values(data_df,
-                                             dataplot)
+            self._get_independent_var_values(data_df, dataplot)
 
         dataset_id = getattr(dataplot, DATASET_ID)
 
@@ -643,8 +648,10 @@ def _data_df(self):
             None else self.simulations_data
 
     @staticmethod
-    def create_subplot(plot_id: str,
-                       subplot_vis_spec: pd.DataFrame) -> Subplot:
+    def create_subplot(
+            plot_id: str,
+            subplot_vis_spec: pd.DataFrame
+    ) -> Subplot:
         """
         Create subplot.
 
@@ -661,7 +668,6 @@ def create_subplot(plot_id: str,
 
                 Subplot
         """
-
         subplot_columns = [col for col in subplot_vis_spec.columns if col in
                            VISUALIZATION_DF_SUBPLOT_LEVEL_COLS]
         subplot = Subplot.from_df(plot_id,
@@ -677,9 +683,10 @@ def create_subplot(plot_id: str,
 
         return subplot
 
-    def parse_from_vis_spec(self,
-                            vis_spec: Optional[Union[str, Path, pd.DataFrame]],
-                            ) -> Tuple[Figure, DataProvider]:
+    def parse_from_vis_spec(
+            self,
+            vis_spec: Optional[Union[str, Path, pd.DataFrame]],
+    ) -> Tuple[Figure, DataProvider]:
         """
         Get visualization settings from a visualization specification.
 
@@ -694,7 +701,6 @@ def parse_from_vis_spec(self,
 
         A figure template with visualization settings and a data provider
         """
-
         # import visualization specification, if file was specified
         if isinstance(vis_spec, (str, Path)):
             vis_spec = core.get_visualization_df(vis_spec)
diff --git a/petab/yaml.py b/petab/yaml.py
index d5c35cba..00d456ef 100644
--- a/petab/yaml.py
+++ b/petab/yaml.py
@@ -28,7 +28,8 @@
 
 def validate(
         yaml_config: Union[Dict, str, Path],
-        path_prefix: Union[None, str, Path] = None):
+        path_prefix: Union[None, str, Path] = None,
+):
     """Validate syntax and semantics of PEtab config YAML
 
     Arguments:
@@ -57,7 +58,7 @@ def validate_yaml_syntax(
             Custom schema for validation
 
     Raises:
-        see jsonschema.validate
+        see :func:`jsonschema.validate`
     """
     yaml_config = load_yaml(yaml_config)
 
@@ -193,7 +194,6 @@ def write_yaml(
         yaml_config: Data to write
         filename: File to create
     """
-
     with open(filename, 'w') as outfile:
         yaml.dump(yaml_config, outfile, default_flow_style=False,
                   sort_keys=False)
@@ -210,8 +210,7 @@ def create_problem_yaml(
         Optional[Union[str, Path, List[Union[str, Path]]]] = None,
         relative_paths: bool = True,
 ) -> None:
-    """
-    Create and write default YAML file for a single PEtab problem
+    """Create and write default YAML file for a single PEtab problem
 
     Arguments:
         sbml_files: Path of SBML model file or list of such
diff --git a/setup.py b/setup.py
index aaa8de2e..0ea76eae 100644
--- a/setup.py
+++ b/setup.py
@@ -29,8 +29,8 @@ def absolute_links(txt):
 
 
 # Python version check. We need >= 3.6 due to e.g. f-strings
-if sys.version_info < (3, 7, 1):
-    sys.exit('PEtab requires at least Python version 3.7.1')
+if sys.version_info < (3, 8, 0):
+    sys.exit('PEtab requires at least Python version 3.8')
 
 # read version from file
 __version__ = ''
@@ -67,7 +67,7 @@ def absolute_links(txt):
                         'jsonschema',
                         ],
       include_package_data=True,
-      python_requires='>=3.7.1',
+      python_requires='>=3.8.0',
       entry_points=ENTRY_POINTS,
       extras_require={
           'tests': [
@@ -79,7 +79,7 @@ def absolute_links(txt):
           'reports': ['Jinja2'],
           'combine': ['python-libcombine>=0.2.6'],
           'doc': [
-              'sphinx>=3.5.3',
+              'sphinx>=3.5.3, !=5.1.0',
               'sphinxcontrib-napoleon>=0.7',
               'sphinx-markdown-tables>=0.0.15',
               'sphinx-rtd-theme>=0.5.1',
diff --git a/tests/test_lint.py b/tests/test_lint.py
index f1935add..0d738813 100644
--- a/tests/test_lint.py
+++ b/tests/test_lint.py
@@ -419,7 +419,14 @@ def test_check_condition_df():
     with pytest.raises(AssertionError):
         lint.check_condition_df(condition_df, model)
 
-    # fix:
+    # fix by adding output parameter
+    observable_df = pd.DataFrame({
+        OBSERVABLE_ID: ["obs1"],
+        OBSERVABLE_FORMULA: ["p1"],
+    })
+    lint.check_condition_df(condition_df, model, observable_df)
+
+    # fix by adding parameter
     ss_model.addParameter('p1', 1.0)
     lint.check_condition_df(condition_df, model)
 
@@ -452,7 +459,7 @@ def test_check_ids():
 
 
 def test_check_parameter_df():
-    """Check parameters.normalize_parameter_df."""
+    """Check lint.check_parameter_df."""
 
     parameter_df = pd.DataFrame({
         PARAMETER_ID: ['par0', 'par1', 'par2'],