diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index e1401abfeaf8e3bcaed32df9e26f6ea6bedc3912..7721cdffe3feca922df1567e5e5b7031f84ec834 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -86,7 +86,8 @@ DynamicModel::DynamicModel(const DynamicModel &m) :
   nonzero_hessian_eqs{m.nonzero_hessian_eqs},
   variableMapping{m.variableMapping},
   blocks_jacob_cols_endo{m.blocks_jacob_cols_endo},
-  var_expectation_functions_to_write{m.var_expectation_functions_to_write}
+  var_expectation_functions_to_write{m.var_expectation_functions_to_write},
+  mfs{m.mfs}
 {
   copyHelper(m);
 }
@@ -133,6 +134,7 @@ DynamicModel::operator=(const DynamicModel &m)
   blocks_jacob_cols_endo = m.blocks_jacob_cols_endo;
 
   var_expectation_functions_to_write = m.var_expectation_functions_to_write;
+  mfs = m.mfs;
 
   copyHelper(m);
 
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 7e53abb52b222bf41dfacacbebd55f26eaf57940..0d2187b12ec8948b36bb77309b2a1e4b017441be 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -109,6 +109,9 @@ private:
   //! Used for var_expectation and var_model
   map<string, set<int>> var_expectation_functions_to_write;
 
+  // Value of the “mfs” option of “model” block (or ”model_options” command)
+  int mfs{1};
+
   // Writes dynamic model file (MATLAB/Octave version, legacy representation)
   void writeDynamicMFile(const string &basename) const;
   //! Writes the code of the block-decomposed model in virtual machine bytecode
@@ -655,6 +658,18 @@ public:
 
   // Returns the set of equations (as numbers) which have a pac_expectation operator
   set<int> findPacExpectationEquationNumbers() const;
+
+  int
+  getMFS() const override
+  {
+    return mfs;
+  }
+
+  void
+  setMFS(int mfs_arg)
+  {
+    mfs = mfs_arg;
+  }
 };
 
 template<bool julia>
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index dbc65069ac3572f053f9234c7c61ef7399f43828..b7a510fd10260e7bd997ca629c84f97e8025a852 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -172,8 +172,7 @@ ModelTree::ModelTree(const ModelTree &m) :
   eq2block{m.eq2block},
   blocks_jacobian_sparse_colptr{m.blocks_jacobian_sparse_colptr},
   endo2eq{m.endo2eq},
-  cutoff{m.cutoff},
-  mfs{m.mfs}
+  cutoff{m.cutoff}
 {
   copyHelper(m);
 }
@@ -221,7 +220,6 @@ ModelTree::operator=(const ModelTree &m)
   blocks_jacobian_sparse_colptr = m.blocks_jacobian_sparse_colptr;
   endo2eq = m.endo2eq;
   cutoff = m.cutoff;
-  mfs = m.mfs;
 
   user_set_add_flags = m.user_set_add_flags;
   user_set_subst_flags = m.user_set_subst_flags;
@@ -539,7 +537,7 @@ ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &fi
               try
                 {
                   normalized_eq = equations[eq]->normalizeEquation(symbol_table.getID(SymbolType::endogenous, var), 0);
-                  if ((mfs == 2 && variable_not_in_derivative) || mfs == 3)
+                  if ((getMFS() == 2 && variable_not_in_derivative) || getMFS() == 3)
                     Equation_Simulation_Type = EquationType::evaluateRenormalized;
                 }
               catch (ExprNode::NormalizationFailed &e)
@@ -706,7 +704,7 @@ ModelTree::computeBlockDecomposition(int prologue, int epilogue)
              || variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0
              || equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0
              || equation_lag_lead[eq_idx_block2orig[i+prologue]].second > 0))
-        || mfs == 0)
+        || getMFS() == 0)
       add_edge(vertex(i, G), vertex(i, G), G);
 
   const vector<int> old_eq_idx_block2orig(eq_idx_block2orig), old_endo_idx_block2orig(endo_idx_block2orig);
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index c463986010236c07e2a39d743a368a73805edf07..aa0f6fd8971fa5bcb575d2dd8656a7e7487f16d2 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -619,8 +619,6 @@ protected:
 public:
   //! Absolute value under which a number is considered to be zero
   double cutoff{1e-15};
-  // Setting for minimum feedback set computation (see the reference manual)
-  int mfs{1};
   //! Declare a node as an equation of the model; also give its line number
   void addEquation(expr_t eq, optional<int> lineno);
   //! Declare a node as an equation of the model, also giving its tags
@@ -703,6 +701,10 @@ public:
         return "UNKNOWN                      ";
       }
   }
+
+  /* Returns the minimum feedback set value (see the “mfs” option of the
+     “model” block in the reference manual for the possible values) */
+  virtual int getMFS() const = 0;
 };
 
 template<ExprNodeOutputType output_type>
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 8868f3cf1701e3a1ac0f4fe8715233354dcb99c4..161cc25ea349f3d6f26c89314bf1ee5c5bb49303 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -697,7 +697,7 @@ void
 ParsingDriver::mfs(const string &value)
 {
   int val = stoi(value);
-  mod_file->dynamic_model.mfs = val;
+  mod_file->dynamic_model.setMFS(val);
 }
 
 void
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 35076bb01773f08cadfe92bc3d9ddada0b364ea1..bfaf0bdaf1fd6ff1b5e9da469bd555df3735fe5a 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -90,7 +90,6 @@ StaticModel::StaticModel(const DynamicModel &m) :
     addAuxEquation(aux_eq->toStatic(*this));
 
   cutoff = m.cutoff;
-  mfs = m.mfs;
 
   user_set_add_flags = m.user_set_add_flags;
   user_set_subst_flags = m.user_set_subst_flags;
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 5e98f28f0db2b67e6b901ef19dd62df7d142e049..f940c50b2856babc24492e34302f25ac08e7ad9a 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -190,6 +190,12 @@ public:
 
   // Writes ramsey_multipliers_derivatives (C version)
   void writeRamseyMultipliersDerivativesCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, int ramsey_orig_endo_nbr) const;
+
+  int
+  getMFS() const override
+  {
+    return 0;
+  }
 };
 
 template<bool julia>