diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 92346ae059f4e7b327becee82968ffe1468eb190..735641d4a14f90eb3271680c66651f7658f6f50f 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -125,6 +125,9 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
   NNZDerivatives(4, 0),
   temporary_terms_derivatives(4)
 {
+  // Ensure that elements accessed by writeParamsDerivativesFileHelper() exist
+  for (const auto &ord : {pair{0, 1}, pair{1, 1}, pair{0, 2}, pair{1, 2}, pair{2, 1}, pair{3, 1}})
+    params_derivatives.emplace(ord, decltype(params_derivatives)::mapped_type{});
 }
 
 ModelTree::ModelTree(const ModelTree &m) :
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index ffef91435551f6db608226f54e84769c5c459faf..2126c45d3a866be2e9261ad055f6f039780c6659 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -928,7 +928,7 @@ ModelTree::writeParamsDerivativesFileHelper() const
     writeTemporaryTerms<output_type>(tts, temp_term_union, params_derivs_temporary_terms_idxs,
                                      tt_output, tef_terms);
 
-  for (const auto &[indices, d1] : params_derivatives.find({ 0, 1 })->second)
+  for (const auto &[indices, d1] : params_derivatives.at({ 0, 1 }))
     {
       auto [eq, param] { vectorToTuple<2>(indices) };
 
@@ -940,7 +940,7 @@ ModelTree::writeParamsDerivativesFileHelper() const
       rp_output << ";" << endl;
     }
 
-  for (const auto &[indices, d2] : params_derivatives.find({ 1, 1 })->second)
+  for (const auto &[indices, d2] : params_derivatives.at({ 1, 1 }))
     {
       auto [eq, var, param] { vectorToTuple<3>(indices) };
 
@@ -954,7 +954,7 @@ ModelTree::writeParamsDerivativesFileHelper() const
     }
 
   for (int i {1};
-       const auto &[indices, d2] : params_derivatives.find({ 0, 2 })->second)
+       const auto &[indices, d2] : params_derivatives.at({ 0, 2 }))
     {
       auto [eq, param1, param2] { vectorToTuple<3>(indices) };
 
@@ -992,7 +992,7 @@ ModelTree::writeParamsDerivativesFileHelper() const
     }
 
   for (int i {1};
-       const auto &[indices, d2] : params_derivatives.find({ 1, 2 })->second)
+       const auto &[indices, d2] : params_derivatives.at({ 1, 2 }))
     {
       auto [eq, var, param1, param2] { vectorToTuple<4>(indices) };
 
@@ -1035,7 +1035,7 @@ ModelTree::writeParamsDerivativesFileHelper() const
     }
 
   for (int i {1};
-       const auto &[indices, d2] : params_derivatives.find({ 2, 1 })->second)
+       const auto &[indices, d2] : params_derivatives.at({ 2, 1 }))
     {
       auto [eq, var1, var2, param] { vectorToTuple<4>(indices) };
 
@@ -1080,7 +1080,7 @@ ModelTree::writeParamsDerivativesFileHelper() const
   if constexpr(output_type == ExprNodeOutputType::matlabDynamicModel
                || output_type == ExprNodeOutputType::juliaDynamicModel)
     for (int i {1};
-         const auto &[indices, d2] : params_derivatives.find({ 3, 1 })->second)
+         const auto &[indices, d2] : params_derivatives.at({ 3, 1 }))
       {
         auto [eq, var1, var2, var3, param] { vectorToTuple<5>(indices) };