diff --git a/python/sdist/amici/conserved_quantities_demartino.py b/python/sdist/amici/conserved_quantities_demartino.py
index 43c84fdf5a..cf651241d1 100644
--- a/python/sdist/amici/conserved_quantities_demartino.py
+++ b/python/sdist/amici/conserved_quantities_demartino.py
@@ -139,7 +139,7 @@ def log(*args, **kwargs):
     # print all conserved quantities
     if verbose:
         for i, (coefficients, engaged_species_idxs) in enumerate(
-            zip(species_coefficients, species_indices)
+            zip(species_coefficients, species_indices, strict=True)
         ):
             if not engaged_species_idxs:
                 continue
@@ -148,7 +148,7 @@ def log(*args, **kwargs):
                 "species:"
             )
             for species_idx, coefficient in zip(
-                engaged_species_idxs, coefficients
+                engaged_species_idxs, coefficients, strict=True
             ):
                 name = (
                     species_names[species_idx]
@@ -957,12 +957,12 @@ def _reduce(
                 k2 = order[j]
                 column: list[float] = [0] * num_species
                 for species_idx, coefficient in zip(
-                    cls_species_idxs[k1], cls_coefficients[k1]
+                    cls_species_idxs[k1], cls_coefficients[k1], strict=True
                 ):
                     column[species_idx] = coefficient
                 ok1 = True
                 for species_idx, coefficient in zip(
-                    cls_species_idxs[k2], cls_coefficients[k2]
+                    cls_species_idxs[k2], cls_coefficients[k2], strict=True
                 ):
                     column[species_idx] -= coefficient
                     if column[species_idx] < -_MIN:
diff --git a/python/sdist/amici/cxxcodeprinter.py b/python/sdist/amici/cxxcodeprinter.py
index d9de0ca8f8..69060eb00a 100644
--- a/python/sdist/amici/cxxcodeprinter.py
+++ b/python/sdist/amici/cxxcodeprinter.py
@@ -179,7 +179,9 @@ def format_regular_line(symbol, math, index):
                 # we need toposort to handle the dependencies of extracted
                 #  subexpressions
                 expr_dict = dict(
-                    itertools.chain(zip(symbols, reduced_exprs), replacements)
+                    itertools.chain(
+                        zip(symbols, reduced_exprs, strict=True), replacements
+                    )
                 )
                 sorted_symbols = toposort(
                     {
@@ -192,7 +194,7 @@ def format_regular_line(symbol, math, index):
                     }
                 )
                 symbol_to_idx = {
-                    sym: idx for idx, sym in zip(indices, symbols)
+                    sym: idx for idx, sym in zip(indices, symbols, strict=True)
                 }
 
                 def format_line(symbol: sp.Symbol):
@@ -217,7 +219,9 @@ def format_line(symbol: sp.Symbol):
 
         return [
             format_regular_line(sym, math, index)
-            for index, sym, math in zip(indices, symbols, equations)
+            for index, sym, math in zip(
+                indices, symbols, equations, strict=True
+            )
             if math not in [0, 0.0]
         ]
 
diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py
index 147ae3d381..6b1392a3d1 100644
--- a/python/sdist/amici/de_export.py
+++ b/python/sdist/amici/de_export.py
@@ -717,7 +717,9 @@ def _get_function_body(
             for ipar in range(self.model.num_par()):
                 expressions = []
                 for index, formula in zip(
-                    self.model._x0_fixedParameters_idx, equations[:, ipar]
+                    self.model._x0_fixedParameters_idx,
+                    equations[:, ipar],
+                    strict=True,
                 ):
                     if not formula.is_zero:
                         expressions.extend(
@@ -735,7 +737,7 @@ def _get_function_body(
 
         elif function == "x0_fixedParameters":
             for index, formula in zip(
-                self.model._x0_fixedParameters_idx, equations
+                self.model._x0_fixedParameters_idx, equations, strict=True
             ):
                 lines.append(
                     f"    if(std::find(reinitialization_state_idxs.cbegin(), "
diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py
index f0b9c6cbb0..0caac635a7 100644
--- a/python/sdist/amici/de_model.py
+++ b/python/sdist/amici/de_model.py
@@ -396,7 +396,9 @@ def states(self) -> list[State]:
     def _process_sbml_rate_of(self) -> None:
         """Substitute any SBML-rateOf constructs in the model equations"""
         rate_of_func = sp.core.function.UndefinedFunction("rateOf")
-        species_sym_to_xdot = dict(zip(self.sym("x"), self.sym("xdot")))
+        species_sym_to_xdot = dict(
+            zip(self.sym("x"), self.sym("xdot"), strict=True)
+        )
         species_sym_to_idx = {x: i for i, x in enumerate(self.sym("x"))}
 
         def get_rate(symbol: sp.Symbol):
@@ -427,7 +429,7 @@ def get_rate(symbol: sp.Symbol):
         if made_substitutions:
             # substitute in topological order
             subs = toposort_symbols(
-                dict(zip(self.sym("xdot"), self.eq("xdot")))
+                dict(zip(self.sym("xdot"), self.eq("xdot"), strict=True))
             )
             self._eqs["xdot"] = smart_subs_dict(self.eq("xdot"), subs)
 
@@ -469,6 +471,7 @@ def get_rate(symbol: sp.Symbol):
                     zip(
                         self.sym("w")[self.num_cons_law() :, :],
                         self.eq("w")[self.num_cons_law() :, :],
+                        strict=True,
                     )
                 )
             )
@@ -477,6 +480,7 @@ def get_rate(symbol: sp.Symbol):
                     zip(
                         self.sym("w")[: self.num_cons_law(), :],
                         self.eq("w")[: self.num_cons_law(), :],
+                        strict=True,
                     )
                 )
                 | w_sorted
@@ -1027,7 +1031,9 @@ def static_indices(self, name: str) -> list[int]:
             #  of non-zeros entries of the sparse matrix
             self._static_indices[name] = [
                 i
-                for i, (expr, row_idx) in enumerate(zip(sparseeq, rowvals))
+                for i, (expr, row_idx) in enumerate(
+                    zip(sparseeq, rowvals, strict=True)
+                )
                 # derivative of a static expression is static
                 if row_idx in static_indices_w
                 # constant expressions
@@ -1396,6 +1402,8 @@ def _compute_equation(self, name: str) -> None:
                                 if not s.has_conservation_law()
                             ),
                             self.sym("dx"),
+                            # dx contains extra elements for algebraic states
+                            strict=False,
                         )
                     ]
                     + [eq.get_val() for eq in self._algebraic_equations]
@@ -1556,7 +1564,9 @@ def _compute_equation(self, name: str) -> None:
             # backsubstitution of optimized right-hand side terms into RHS
             # calling subs() is costly. Due to looping over events though, the
             # following lines are only evaluated if a model has events
-            w_sorted = toposort_symbols(dict(zip(self.sym("w"), self.eq("w"))))
+            w_sorted = toposort_symbols(
+                dict(zip(self.sym("w"), self.eq("w"), strict=True))
+            )
             tmp_xdot = smart_subs_dict(self.eq("xdot"), w_sorted)
             self._eqs[name] = self.eq("drootdt")
             if self.num_states_solver():
@@ -1586,7 +1596,7 @@ def _compute_equation(self, name: str) -> None:
                 for event_obs in self._event_observables
             ]
             for (iz, ie), event_obs in zip(
-                enumerate(z2event), self._event_observables
+                enumerate(z2event), self._event_observables, strict=True
             ):
                 event_observables[ie - 1][iz] = event_obs.get_val()
 
@@ -1727,7 +1737,7 @@ def _compute_equation(self, name: str) -> None:
             syms_x = self.sym("x")
             syms_yz = self.sym(name.removeprefix("sigma"))
             xs_in_sigma = {}
-            for sym_yz, eq_yz in zip(syms_yz, self._eqs[name]):
+            for sym_yz, eq_yz in zip(syms_yz, self._eqs[name], strict=True):
                 yz_free_syms = eq_yz.free_symbols
                 if tmp := {x for x in syms_x if x in yz_free_syms}:
                     xs_in_sigma[sym_yz] = tmp
@@ -2229,6 +2239,7 @@ def _collect_heaviside_roots(
                 zip(
                     [expr.get_id() for expr in self._expressions],
                     [expr.get_val() for expr in self._expressions],
+                    strict=True,
                 )
             )
         )
diff --git a/python/sdist/amici/pandas.py b/python/sdist/amici/pandas.py
index b61cfa0c79..e83e524e08 100644
--- a/python/sdist/amici/pandas.py
+++ b/python/sdist/amici/pandas.py
@@ -155,7 +155,7 @@ def getSimulationObservablesAsDataFrame(
 
     # aggregate records
     dicts = []
-    for edata, rdata in zip(edata_list, rdata_list):
+    for edata, rdata in zip(edata_list, rdata_list, strict=True):
         for i_time, timepoint in enumerate(rdata["t"]):
             datadict = {
                 "time": timepoint,
@@ -212,7 +212,7 @@ def getSimulationStatesAsDataFrame(
 
     # aggregate records
     dicts = []
-    for edata, rdata in zip(edata_list, rdata_list):
+    for edata, rdata in zip(edata_list, rdata_list, strict=True):
         for i_time, timepoint in enumerate(rdata["t"]):
             datadict = {
                 "time": timepoint,
@@ -268,7 +268,7 @@ def get_expressions_as_dataframe(
 
     # aggregate records
     dicts = []
-    for edata, rdata in zip(edata_list, rdata_list):
+    for edata, rdata in zip(edata_list, rdata_list, strict=True):
         for i_time, timepoint in enumerate(rdata["t"]):
             datadict = {
                 "time": timepoint,
diff --git a/python/sdist/amici/petab/conditions.py b/python/sdist/amici/petab/conditions.py
index 61b8679765..831420f68e 100644
--- a/python/sdist/amici/petab/conditions.py
+++ b/python/sdist/amici/petab/conditions.py
@@ -69,7 +69,9 @@ def fill_in_parameters(
             RuntimeWarning,
         )
 
-    for edata, mapping_for_condition in zip(edatas, parameter_mapping):
+    for edata, mapping_for_condition in zip(
+        edatas, parameter_mapping, strict=True
+    ):
         fill_in_parameters_for_condition(
             edata,
             problem_parameters,
diff --git a/python/sdist/amici/petab/parameter_mapping.py b/python/sdist/amici/petab/parameter_mapping.py
index 369ad13a96..54930073da 100644
--- a/python/sdist/amici/petab/parameter_mapping.py
+++ b/python/sdist/amici/petab/parameter_mapping.py
@@ -385,7 +385,7 @@ def create_parameter_mapping(
 
     parameter_mapping = ParameterMapping()
     for (_, condition), prelim_mapping_for_condition in zip(
-        simulation_conditions.iterrows(), prelim_parameter_mapping
+        simulation_conditions.iterrows(), prelim_parameter_mapping, strict=True
     ):
         mapping_for_condition = create_parameter_mapping_for_condition(
             prelim_mapping_for_condition, condition, petab_problem, amici_model
diff --git a/python/sdist/amici/petab/pysb_import.py b/python/sdist/amici/petab/pysb_import.py
index b770603935..f5b1c84dbf 100644
--- a/python/sdist/amici/petab/pysb_import.py
+++ b/python/sdist/amici/petab/pysb_import.py
@@ -55,6 +55,7 @@ def _add_observation_model(
         petab_problem.observable_df.index,
         petab_problem.observable_df[OBSERVABLE_FORMULA],
         petab_problem.observable_df[NOISE_FORMULA],
+        strict=True,
     ):
         obs_symbol = sp.sympify(observable_formula, locals=local_syms)
         if observable_id in pysb_model.expressions.keys():
diff --git a/python/sdist/amici/petab/simulations.py b/python/sdist/amici/petab/simulations.py
index 77f14ccb35..e80e40a94b 100644
--- a/python/sdist/amici/petab/simulations.py
+++ b/python/sdist/amici/petab/simulations.py
@@ -172,6 +172,7 @@ def simulate_petab(
             zip(
                 petab_problem.x_ids,
                 petab_problem.x_nominal_scaled,
+                strict=True,
             )
         )
         # depending on `fill_fixed_parameters` for parameter mapping, the
@@ -311,7 +312,7 @@ def aggregate_sllh(
             )
 
     for condition_parameter_mapping, edata, rdata in zip(
-        parameter_mapping, edatas, rdatas
+        parameter_mapping, edatas, rdatas, strict=True
     ):
         for sllh_parameter_index, condition_parameter_sllh in enumerate(
             rdata.sllh
@@ -433,7 +434,9 @@ def rdatas_to_measurement_df(
     observable_ids = model.getObservableIds()
     rows = []
     # iterate over conditions
-    for (_, condition), rdata in zip(simulation_conditions.iterrows(), rdatas):
+    for (_, condition), rdata in zip(
+        simulation_conditions.iterrows(), rdatas, strict=True
+    ):
         # current simulation matrix
         y = rdata.y
         # time array used in rdata
diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py
index 651ba0e23e..3067d26bc5 100644
--- a/python/sdist/amici/plotting.py
+++ b/python/sdist/amici/plotting.py
@@ -65,7 +65,7 @@ def plot_state_trajectories(
     else:
         labels = np.asarray(rdata.ptr.state_ids)[list(state_indices)]
 
-    for ix, label in zip(state_indices, labels):
+    for ix, label in zip(state_indices, labels, strict=True):
         ax.plot(rdata["t"], rdata["x"][:, ix], marker=marker, label=label)
 
     ax.set_xlabel("$t$")
@@ -131,7 +131,7 @@ def plot_observable_trajectories(
     else:
         labels = np.asarray(rdata.ptr.observable_ids)[list(observable_indices)]
 
-    for iy, label in zip(observable_indices, labels):
+    for iy, label in zip(observable_indices, labels, strict=True):
         (l,) = ax.plot(
             rdata["t"], rdata["y"][:, iy], marker=marker, label=label
         )
diff --git a/python/sdist/amici/pysb_import.py b/python/sdist/amici/pysb_import.py
index 46476a753d..1a21fef1ca 100644
--- a/python/sdist/amici/pysb_import.py
+++ b/python/sdist/amici/pysb_import.py
@@ -563,7 +563,11 @@ def _add_expression(
         cost_fun_expr = sp.sympify(
             cost_fun_str,
             locals=dict(
-                zip(_get_str_symbol_identifiers(name), (y, my, sigma))
+                zip(
+                    _get_str_symbol_identifiers(name),
+                    (y, my, sigma),
+                    strict=True,
+                )
             ),
         )
         ode_model.add_component(
diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py
index 40c4881a09..843b954a8f 100644
--- a/python/sdist/amici/sbml_import.py
+++ b/python/sdist/amici/sbml_import.py
@@ -556,9 +556,18 @@ def _build_ode_model(
         dxdt = smart_multiply(
             self.stoichiometric_matrix, MutableDenseMatrix(fluxes)
         )
+        # dxdt has algebraic states at the end
+        assert dxdt.shape[0] - len(self.symbols[SymbolId.SPECIES]) == len(
+            self.symbols.get(SymbolId.ALGEBRAIC_STATE, [])
+        ), (
+            self.symbols[SymbolId.SPECIES],
+            dxdt,
+            self.symbols[SymbolId.SPECIES],
+        )
+
         # correct time derivatives for compartment changes
         for ix, ((species_id, species), formula) in enumerate(
-            zip(self.symbols[SymbolId.SPECIES].items(), dxdt)
+            zip(self.symbols[SymbolId.SPECIES].items(), dxdt, strict=False)
         ):
             # rate rules and amount species don't need to be updated
             if "dt" in species:
@@ -604,7 +613,7 @@ def _build_ode_model(
 
         # add fluxes as expressions, this needs to happen after base
         # expressions from symbols have been parsed
-        for flux_id, flux in zip(fluxes, self.flux_vector):
+        for flux_id, flux in zip(fluxes, self.flux_vector, strict=True):
             # replace splines inside fluxes
             flux = flux.subs(spline_subs)
             ode_model.add_component(
@@ -981,7 +990,9 @@ def _process_species_initial(self):
             self.symbols[SymbolId.SPECIES], "init"
         )
         for species, rateof_dummies in zip(
-            self.symbols[SymbolId.SPECIES].values(), all_rateof_dummies
+            self.symbols[SymbolId.SPECIES].values(),
+            all_rateof_dummies,
+            strict=True,
         ):
             species["init"] = _dummy_to_rateof(
                 smart_subs_dict(species["init"], sorted_species, "init"),
@@ -1945,6 +1956,7 @@ def _process_log_likelihood(
         for (obs_id, obs), (sigma_id, sigma) in zip(
             self.symbols[obs_symbol].items(),
             self.symbols[sigma_symbol].items(),
+            strict=True,
         ):
             symbol = symbol_with_assumptions(f"J{obs_id}")
             dist = noise_distributions.get(str(obs_id), "normal")
@@ -1955,6 +1967,7 @@ def _process_log_likelihood(
                     zip(
                         _get_str_symbol_identifiers(obs_id),
                         (obs_id, obs["measurement_symbol"], sigma_id),
+                        strict=True,
                     )
                 ),
             )
@@ -2173,9 +2186,9 @@ def _get_conservation_laws_demartino(
             len(cls_coefficients), len(ode_model._differential_states)
         )
         for i_cl, (cl, coefficients) in enumerate(
-            zip(cls_state_idxs, cls_coefficients)
+            zip(cls_state_idxs, cls_coefficients, strict=True)
         ):
-            for i, c in zip(cl, coefficients):
+            for i, c in zip(cl, coefficients, strict=True):
                 A[i_cl, i] = sp.Rational(c)
         rref, pivots = A.rref()
 
@@ -2319,7 +2332,10 @@ def _add_conservation_for_non_constant_species(
                     "coefficients": {
                         state_id: coeff * compartment
                         for state_id, coeff, compartment in zip(
-                            state_ids, coefficients, compartment_sizes
+                            state_ids,
+                            coefficients,
+                            compartment_sizes,
+                            strict=True,
                         )
                     },
                 }
diff --git a/python/sdist/amici/splines.py b/python/sdist/amici/splines.py
index cd500a4570..5d423a19e5 100644
--- a/python/sdist/amici/splines.py
+++ b/python/sdist/amici/splines.py
@@ -594,7 +594,9 @@ def check_if_valid(self, importer: sbml_import.SbmlImporter) -> None:
             importer.symbols[SymbolId.FIXED_PARAMETER][fp]["value"]
             for fp in fixed_parameters
         ]
-        subs = dict(zip(fixed_parameters, fixed_parameters_values))
+        subs = dict(
+            zip(fixed_parameters, fixed_parameters_values, strict=True)
+        )
         nodes_values = [sp.simplify(x.subs(subs)) for x in self.nodes]
         for x in nodes_values:
             assert x.is_Number
@@ -1093,7 +1095,7 @@ def add_to_sbml_model(
                     # It makes no sense to give a single nominal value:
                     # grid values must all be different
                     raise TypeError("x_nominal must be a Sequence!")
-                for _x, _val in zip(self.nodes, x_nominal):
+                for _x, _val in zip(self.nodes, x_nominal, strict=True):
                     if _x.is_Symbol and not model.getParameter(_x.name):
                         add_parameter(
                             model, _x.name, value=_val, units=x_units
@@ -1116,7 +1118,7 @@ def add_to_sbml_model(
                 else:
                     y_constant = len(self.values_at_nodes) * [y_constant]
                 for _y, _val, _const in zip(
-                    self.values_at_nodes, y_nominal, y_constant
+                    self.values_at_nodes, y_nominal, y_constant, strict=True
                 ):
                     if _y.is_Symbol and not model.getParameter(_y.name):
                         add_parameter(
diff --git a/python/sdist/amici/sympy_utils.py b/python/sdist/amici/sympy_utils.py
index 2c8599ba7e..9794fadef0 100644
--- a/python/sdist/amici/sympy_utils.py
+++ b/python/sdist/amici/sympy_utils.py
@@ -184,7 +184,11 @@ def _parallel_applyfunc(obj: sp.Matrix, func: Callable) -> sp.Matrix:
             elif isinstance(obj, sp.SparseMatrix):
                 dok = obj.todok()
                 mapped = p.map(func, dok.values())
-                dok = {k: v for k, v in zip(dok.keys(), mapped) if v != 0}
+                dok = {
+                    k: v
+                    for k, v in zip(dok.keys(), mapped, strict=True)
+                    if v != 0
+                }
                 return obj._new(obj.rows, obj.cols, dok)
             else:
                 raise ValueError(f"Unsupported matrix type {type(obj)}")
diff --git a/python/tests/petab/test_petab_problem.py b/python/tests/petab/test_petab_problem.py
index 5a8a299bb9..736d895f8f 100644
--- a/python/tests/petab/test_petab_problem.py
+++ b/python/tests/petab/test_petab_problem.py
@@ -81,7 +81,7 @@ def test_amici_petab_problem_pregenerate_equals_on_demand():
     app_store_false.set_parameters(parameter_update, scaled_parameters=True)
 
     for edata_store_true, edata_store_false in zip(
-        app_store_true.get_edatas(), app_store_false.get_edatas()
+        app_store_true.get_edatas(), app_store_false.get_edatas(), strict=True
     ):
         assert edata_store_true is not edata_store_false
         assert edata_store_true == edata_store_false
diff --git a/python/tests/splines_utils.py b/python/tests/splines_utils.py
index 6565ac80a4..99c0da880c 100644
--- a/python/tests/splines_utils.py
+++ b/python/tests/splines_utils.py
@@ -216,7 +216,7 @@ def create_petab_problem(
     zz_true = np.array(
         [
             integrate_spline(spline, params_true, tt_obs, iv)
-            for (spline, iv) in zip(splines, initial_values)
+            for (spline, iv) in zip(splines, initial_values, strict=True)
         ],
         dtype=float,
     )
@@ -480,7 +480,7 @@ def compute_ground_truth(
     x_true_sym = sp.Matrix(
         [
             integrate_spline(spline, None, times, iv)
-            for (spline, iv) in zip(splines, initial_values)
+            for (spline, iv) in zip(splines, initial_values, strict=True)
         ]
     ).transpose()
     groundtruth = {
@@ -604,7 +604,7 @@ def param_by_name(id):
         x_true_sym = sp.Matrix(
             [
                 integrate_spline(spline, None, tt, iv)
-                for (spline, iv) in zip(splines, initial_values)
+                for (spline, iv) in zip(splines, initial_values, strict=True)
             ]
         ).transpose()
         x_true = np.asarray(x_true_sym.subs(params_true), dtype=float)
@@ -896,7 +896,7 @@ def example_spline_1(
     yy = list(sp.symbols(f"y{idx}_0:{len(yy_true)}"))
 
     if fixed_values is None:
-        params = dict(zip(yy, yy_true))
+        params = dict(zip(yy, yy_true, strict=True))
     elif fixed_values == "all":
         params = {}
         for i in range(len(yy_true)):
@@ -939,7 +939,7 @@ def example_spline_2(idx: int = 0):
     xx = UniformGrid(0, 25, number_of_nodes=len(yy_true))
     yy = list(sp.symbols(f"y{idx}_0:{len(yy_true) - 1}"))
     yy.append(yy[0])
-    params = dict(zip(yy, yy_true))
+    params = dict(zip(yy, yy_true, strict=True))
     spline = CubicHermiteSpline(
         f"y{idx}",
         nodes=xx,
@@ -960,7 +960,7 @@ def example_spline_3(idx: int = 0):
     yy_true = [0.0, 2.0, 5.0, 6.0, 5.0, 4.0, 2.0, 3.0, 4.0, 6.0]
     xx = UniformGrid(0, 25, number_of_nodes=len(yy_true))
     yy = list(sp.symbols(f"y{idx}_0:{len(yy_true)}"))
-    params = dict(zip(yy, yy_true))
+    params = dict(zip(yy, yy_true, strict=True))
     spline = CubicHermiteSpline(
         f"y{idx}",
         nodes=xx,
diff --git a/python/tests/test_petab_objective.py b/python/tests/test_petab_objective.py
index 5d29ad88ff..07770ac413 100755
--- a/python/tests/test_petab_objective.py
+++ b/python/tests/test_petab_objective.py
@@ -40,10 +40,7 @@ def test_simulate_petab_sensitivities(lotka_volterra):
     amici_solver.setMaxSteps(int(1e5))
 
     problem_parameters = dict(
-        zip(
-            petab_problem.x_ids,
-            petab_problem.x_nominal,
-        )
+        zip(petab_problem.x_ids, petab_problem.x_nominal, strict=True)
     )
 
     results = {}
diff --git a/python/tests/test_splines.py b/python/tests/test_splines.py
index a7fe01e84a..66104b824b 100644
--- a/python/tests/test_splines.py
+++ b/python/tests/test_splines.py
@@ -62,7 +62,7 @@ def test_multiple_splines(**kwargs):
 
     tols = []
     for t0, t1, t2, t3, t4, t5 in zip(
-        tols0, tols1, tols2, tols3, tols4, tols5
+        tols0, tols1, tols2, tols3, tols4, tols5, strict=True
     ):
         keys = set().union(
             t0.keys(), t1.keys(), t2.keys(), t3.keys(), t4.keys(), t5.keys()
diff --git a/python/tests/test_splines_python.py b/python/tests/test_splines_python.py
index 4c4de5ccfc..539fb1dd4d 100644
--- a/python/tests/test_splines_python.py
+++ b/python/tests/test_splines_python.py
@@ -215,8 +215,11 @@ def test_SplineNonUniformPeriodicExtrapolation():
 @skip_on_valgrind
 def check_gradient(spline, t, params, params_values, expected, rel_tol=1e-9):
     value = spline.evaluate(t)
-    subs = {pname: pvalue for (pname, pvalue) in zip(params, params_values)}
-    for p, exp in zip(params, expected):
+    subs = {
+        pname: pvalue
+        for (pname, pvalue) in zip(params, params_values, strict=True)
+    }
+    for p, exp in zip(params, expected, strict=True):
         assert math.isclose(
             float(value.diff(p).subs(subs)), exp, rel_tol=rel_tol
         )
diff --git a/python/tests/test_splines_short.py b/python/tests/test_splines_short.py
index 59e54a3279..43b0ccd294 100644
--- a/python/tests/test_splines_short.py
+++ b/python/tests/test_splines_short.py
@@ -44,7 +44,7 @@ def test_two_splines(**kwargs):
         tols1 = (tols1, tols1, tols1)
 
     tols = []
-    for t0, t1 in zip(tols0, tols1):
+    for t0, t1 in zip(tols0, tols1, strict=True):
         keys = set().union(t0.keys(), t1.keys())
         t = {
             key: max(
diff --git a/python/tests/util.py b/python/tests/util.py
index 58e17137a6..b450501a81 100644
--- a/python/tests/util.py
+++ b/python/tests/util.py
@@ -109,7 +109,7 @@ def create_event_assignment(target, assignment):
 
         if isinstance(event_def["target"], list):
             for event_target, event_assignment in zip(
-                event_def["target"], event_def["assignment"]
+                event_def["target"], event_def["assignment"], strict=True
             ):
                 create_event_assignment(event_target, event_assignment)