1310 Spatially resolved feedback simulation

README.md

@@ -19,7 +19,8 @@ and, in particular, for
 - Hybrid agent-metapopulation-based models: Bicker J, Schmieding R, Meyer-Hermann M, Kühn MJ. (2025). *Hybrid metapopulation agent-based epidemiological models for efficient insight on the individual scale: A contribution to green computing*. *Infectious Disease Modelling* 10(2): 571-590. `https://doi.org/10.1016/j.idm.2024.12.015`
 - Graph Neural Networks: Schmidt A, Zunker H, Heinlein A, Kühn MJ. (2024). *Towards Graph Neural Network Surrogates Leveraging Mechanistic Expert Knowledge for Pandemic Response*. arXiv. `https://arxiv.org/abs/2411.06500`
 - ODE-based models with Linear Chain Trick: Plötzke L, Wendler A, Schmieding R, Kühn MJ. (2024). *Revisiting the Linear Chain Trick in epidemiological models: Implications of underlying assumptions for numerical solutions*. Submitted for publication. `https://doi.org/10.48550/arXiv.2412.09140`
-- Behavior-based ODE models: Zunker H, Dönges P, Lenz P, Contreras S, Kühn MJ. (2025). *Risk-mediated dynamic regulation of effective contacts de-synchronizes outbreaks in metapopulation epidemic models*. arXiv. `https://arxiv.org/abs/2502.14428`
+- Behavior-based ODE models: Zunker H, Dönges P, Lenz P, Contreras S, Kühn MJ. (2025). *Risk-mediated dynamic regulation of effective contacts de-synchronizes outbreaks in metapopulation epidemic models*. Chaos, Solitons & Fractals. `https://doi.org/10.1016/j.chaos.2025.116782`
+
 
 **Getting started**
 

cpp/examples/CMakeLists.txt

@@ -80,6 +80,10 @@ add_executable(ode_secir_feedback ode_secir_feedback.cpp)
 target_link_libraries(ode_secir_feedback PRIVATE memilio ode_secir)
 target_compile_options(ode_secir_feedback PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(ode_secir_feedback_graph ode_secir_feedback_graph.cpp)
+target_link_libraries(ode_secir_feedback_graph PRIVATE memilio ode_secir)
+target_compile_options(ode_secir_feedback_graph PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 add_executable(ode_secirvvs_example ode_secirvvs.cpp)
 target_link_libraries(ode_secirvvs_example PRIVATE memilio ode_secirvvs)
 target_compile_options(ode_secirvvs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

cpp/examples/ode_secir_feedback_graph.cpp

@@ -0,0 +1,177 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Henrik Zunker
+*
+* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+#include "memilio/data/analyze_result.h"
+#include "ode_secir/model.h"
+#include "memilio/compartments/feedback_simulation.h"
+#include "memilio/mobility/metapopulation_mobility_instant.h"
+#include "memilio/mobility/graph_simulation.h"
+#include "memilio/utils/logging.h"
+#include <iostream>
+
+// alias for the type of the simulation with feedback
+using FeedbackSim = mio::FeedbackSimulation<double, mio::Simulation<double, mio::osecir::Model<double>>,
+                                            mio::osecir::ContactPatterns<double>>;
+
+// helper function to initialize the model with population and parameters
+void initialize_model(mio::osecir::Model<double>& model, double cont_freq)
+{
+    model.parameters.set<mio::osecir::StartDay>(60);
+    model.parameters.set<mio::osecir::Seasonality<double>>(0.2);
+
+    // Mean stay times per compartment
+    model.parameters.get<mio::osecir::TimeExposed<double>>()            = 3.2;
+    model.parameters.get<mio::osecir::TimeInfectedNoSymptoms<double>>() = 2.0;
+    model.parameters.get<mio::osecir::TimeInfectedSymptoms<double>>()   = 5.8;
+    model.parameters.get<mio::osecir::TimeInfectedSevere<double>>()     = 9.5;
+    model.parameters.get<mio::osecir::TimeInfectedCritical<double>>()   = 7.1;
+
+    // Set transmission and isolation parameters
+    model.parameters.get<mio::osecir::TransmissionProbabilityOnContact<double>>()  = 0.05;
+    model.parameters.get<mio::osecir::RelativeTransmissionNoSymptoms<double>>()    = 0.7;
+    model.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<double>>()    = 0.09;
+    model.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic<double>>()    = 0.25;
+    model.parameters.get<mio::osecir::MaxRiskOfInfectionFromSymptomatic<double>>() = 0.45;
+    model.parameters.get<mio::osecir::TestAndTraceCapacity<double>>()              = 35;
+    model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<double>>()         = 0.2;
+    model.parameters.get<mio::osecir::CriticalPerSevere<double>>()                 = 0.25;
+    model.parameters.get<mio::osecir::DeathsPerCritical<double>>()                 = 0.3;
+
+    // contact matrix
+    mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::osecir::ContactPatterns<double>>();
+    contact_matrix[0]                       = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
+}
+
+// helper function to initialize the feedback mechanism parameters for a simulation
+void initialize_feedback(FeedbackSim& feedback_simulation)
+{
+    // nominal ICU capacity
+    feedback_simulation.get_parameters().template get<mio::NominalICUCapacity<double>>() = 10;
+
+    // ICU occupancy in the past for memory kernel
+    auto& icu_occupancy     = feedback_simulation.get_parameters().template get<mio::ICUOccupancyHistory<double>>();
+    Eigen::VectorXd icu_day = Eigen::VectorXd::Constant(1, 1);
+    const auto cutoff       = static_cast<int>(feedback_simulation.get_parameters().template get<mio::GammaCutOff>());
+    for (int t = -cutoff; t <= 0; ++t) {
+        icu_occupancy.add_time_point(t, icu_day);
+    }
+
+    // bounds for contact reduction measures
+    feedback_simulation.get_parameters().template get<mio::ContactReductionMin<double>>() = {0.1};
+    feedback_simulation.get_parameters().template get<mio::ContactReductionMax<double>>() = {0.8};
+
+    // Set blending factors. The global blending factor is implicitly defined as 1 - local - regional.
+    feedback_simulation.get_parameters().template get<mio::BlendingFactorLocal<double>>()    = 0.5;
+    feedback_simulation.get_parameters().template get<mio::BlendingFactorRegional<double>>() = 0.3;
+}
+
+// helper function to create the graph with nodes and edges
+mio::Graph<mio::SimulationNode<FeedbackSim>, mio::MobilityEdge<double>>
+create_graph(int num_nodes, int total_population, double cont_freq)
+{
+    // Create a graph for the metapopulation simulation
+    mio::Graph<mio::SimulationNode<FeedbackSim>, mio::MobilityEdge<double>> g;
+
+    // Create models and add nodes to the graph
+    for (int i = 0; i < num_nodes; ++i) {
+        mio::osecir::Model<double> model(1);
+        initialize_model(model, cont_freq);
+
+        // Set initial populations (infection starts in the first node)
+        if (i == 0) {
+            model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = total_population * 0.1;
+            model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] =
+                total_population * 0.1;
+            model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] =
+                total_population * 0.05;
+            model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSevere}] =
+                total_population * 0.02;
+            model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] =
+                total_population * 0.01;
+            model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = 0;
+            model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
+                                                        total_population);
+        }
+        else {
+            model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}] = total_population;
+        }
+        model.apply_constraints();
+
+        // Determine the index for the ICU state (InfectedCritical) for the feedback mechanism
+        auto icu_index = std::vector<size_t>{
+            model.populations.get_flat_index({mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical})};
+
+        // Create feedback simulation
+        auto feedback_sim = FeedbackSim(mio::Simulation<double, mio::osecir::Model<double>>(model), icu_index);
+        initialize_feedback(feedback_sim);
+
+        // Node-ID-Logic: 1001-1005, 2001-2005, ...
+        const int region_id = i / 5;
+        const int local_id  = i % 5;
+        const int node_id   = (region_id + 1) * 1000 + (local_id + 1);
+        g.add_node(node_id, std::move(feedback_sim));
+    }
+
+    // Define complete graph, i.e. each node is connected to every other node
+    std::vector<std::vector<size_t>> mobile_compartments(2);
+    for (size_t i = 0; i < g.nodes().size(); ++i) {
+        for (size_t j = 0; j < g.nodes().size(); ++j) {
+            if (i != j) {
+                g.add_edge(i, j, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.01));
+            }
+        }
+    }
+
+    return g;
+}
+
+int main()
+{
+    // This example demonstrates the implementation of a feedback mechanism for a ODE SECIR model in a graph.
+    // It shows how the perceived risk dynamically impacts contact reduction measures in different regions (nodes).
+    mio::set_log_level(mio::LogLevel::err);
+
+    const auto t0              = 0.;
+    const auto tmax            = 10.;
+    const auto dt              = 0.5;
+    const int total_population = 1000;
+    const double cont_freq     = 2.7;
+    const int num_nodes        = 10;
+
+    // Create the graph
+    auto g = create_graph(num_nodes, total_population, cont_freq);
+
+    // Create and run the simulation
+    using Graph = decltype(g);
+    auto sim    = mio::FeedbackGraphSimulation<double, Graph>(std::move(g), t0, dt);
+    sim.advance(tmax);
+
+    // The output shows the compartments sizes for a node without any initial infections.
+    auto& node         = sim.get_graph().nodes()[1];
+    auto& results_node = node.property.get_simulation().get_result();
+    // interpolate results
+    auto interpolated_results_node = mio::interpolate_simulation_result(results_node);
+
+    // print result with print_table
+    std::cout << "Node ID: " << node.id << "\n";
+    std::vector<std::string> cols = {"S", "E", "C", "C_confirmed", "I", "I_confirmed", "H", "U", "R", "D"};
+    interpolated_results_node.print_table(cols);
+
+    return 0;
+}

cpp/memilio/compartments/feedback_simulation.h

@@ -20,6 +20,7 @@
 #ifndef FEEDBACK_SIMULATION_H
 #define FEEDBACK_SIMULATION_H
 
+#include <cassert>
 #include "memilio/compartments/simulation.h"
 #include "memilio/utils/time_series.h"
 #include "memilio/utils/parameter_set.h"
@@ -160,11 +161,44 @@ struct NominalICUCapacity {
     }
 };
 
+/**
+ * @brief Blending factor for local ICU occupancy. 
+ * If we only have a local simulation, this factor is set to 1.0 per default.
+ */
+template <typename FP>
+struct BlendingFactorLocal {
+    using Type = FP;
+    static Type get_default(AgeGroup)
+    {
+        return FP(1.0);
+    }
+    static std::string name()
+    {
+        return "BlendingFactorLocal";
+    }
+};
+
+/**
+ * @brief Blending factor for regional ICU occupancy.
+ */
+template <typename FP>
+struct BlendingFactorRegional {
+    using Type = FP;
+    static Type get_default(AgeGroup)
+    {
+        return FP(0.0);
+    }
+    static std::string name()
+    {
+        return "BlendingFactorRegional";
+    }
+};
+
 template <typename FP>
 using FeedbackSimulationParameters =
     ParameterSet<ICUOccupancyHistory<FP>, GammaShapeParameter<FP>, GammaScaleParameter<FP>, GammaCutOff,
                  ContactReductionMax<FP>, ContactReductionMin<FP>, SoftPlusCurvatureParameter<FP>,
-                 NominalICUCapacity<FP>>;
+                 NominalICUCapacity<FP>, BlendingFactorLocal<FP>, BlendingFactorRegional<FP>>;
 
 /**
  * @brief A generic feedback simulation extending existing simulations with a feedback mechanism.
@@ -181,6 +215,8 @@ template <typename FP, typename Sim, typename ContactPatterns>
 class FeedbackSimulation
 {
 public:
+    using Model = typename Sim::Model;
+
     /**
      * @brief Constructs the FeedbackSimulation by taking ownership of an existing simulation instance.
      *
@@ -191,7 +227,10 @@ class FeedbackSimulation
         : m_simulation(std::move(sim))
         , m_icu_indices(icu_indices)
         , m_feedback_parameters(m_simulation.get_model().parameters.get_num_groups())
-        , m_perceived_risk(static_cast<size_t>(m_simulation.get_model().parameters.get_num_groups()))
+        , m_perceived_risk(static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
+        , m_regional_icu_occupancy(
+              static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
+        , m_global_icu_occupancy(static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
     {
     }
 
@@ -238,6 +277,11 @@ class FeedbackSimulation
         return m_simulation.get_model();
     }
 
+    const auto& get_model() const
+    {
+        return m_simulation.get_model();
+    }
+
     /**
      * @brief Returns the perceived risk time series.
      */
@@ -268,17 +312,40 @@ class FeedbackSimulation
      */
     FP calc_risk_perceived()
     {
-        const auto& icu_occ    = m_feedback_parameters.template get<ICUOccupancyHistory<FP>>();
-        size_t num_time_points = icu_occ.get_num_time_points();
+        const auto& local_icu_occ = m_feedback_parameters.template get<ICUOccupancyHistory<FP>>();
+        size_t num_time_points    = local_icu_occ.get_num_time_points();
         size_t n = std::min(static_cast<size_t>(num_time_points), m_feedback_parameters.template get<GammaCutOff>());
         FP perceived_risk = 0.0;
         const auto& a     = m_feedback_parameters.template get<GammaShapeParameter<FP>>();
         const auto& b     = m_feedback_parameters.template get<GammaScaleParameter<FP>>();
+
+        const auto& blending_local    = m_feedback_parameters.template get<BlendingFactorLocal<FP>>();
+        const auto& blending_regional = m_feedback_parameters.template get<BlendingFactorRegional<FP>>();
+        const auto& blending_global   = 1 - blending_local - blending_regional;
+
+        // assert that that each blending factor is positive
+        assert(blending_local >= 0 && blending_regional >= 0 && blending_global >= 0 &&
+               "Blending factors must be non-negative.");
+
         for (size_t i = num_time_points - n; i < num_time_points; ++i) {
             size_t day   = i - (num_time_points - n);
             FP gamma     = std::pow(b, a) * std::pow(day, a - 1) * std::exp(-b * day) / std::tgamma(a);
-            FP perc_risk = icu_occ.get_value(i).sum() / m_feedback_parameters.template get<NominalICUCapacity<FP>>();
-            perc_risk    = std::min(perc_risk, 1.0);
+            FP perc_risk = 0.0;
+
+            if (blending_local > 0) {
+                perc_risk += blending_local * local_icu_occ.get_value(i).sum() /
+                             m_feedback_parameters.template get<NominalICUCapacity<FP>>();
+            }
+            if (blending_regional > 0 && m_regional_icu_occupancy.get_num_time_points() > 0) {
+                perc_risk += blending_regional * m_regional_icu_occupancy.get_value(i).sum() /
+                             m_feedback_parameters.template get<NominalICUCapacity<FP>>();
+            }
+            if (blending_global > 0 && m_global_icu_occupancy.get_num_time_points() > 0) {
+                perc_risk += blending_global * m_global_icu_occupancy.get_value(i).sum() /
+                             m_feedback_parameters.template get<NominalICUCapacity<FP>>();
+            }
+
+            perc_risk = std::min(perc_risk, 1.0);
             perceived_risk += perc_risk * gamma;
         }
         return perceived_risk;
@@ -307,6 +374,27 @@ class FeedbackSimulation
         m_feedback_parameters.template get<ICUOccupancyHistory<FP>>().add_time_point(t, icu_occ);
     }
 
+    /**
+     * @brief Sets the regional ICU occupancy time series. 
+     * This is used in the graph simulation to initialize and update the regional ICU occupancy.
+     * 
+     * @param icu_regional The regional ICU occupancy time series.
+     */
+    void set_regional_icu_occupancy(const mio::TimeSeries<FP>& icu_regional)
+    {
+        m_regional_icu_occupancy = icu_regional;
+    }
+
+    /**
+     * @brief Sets the global ICU occupancy time series.
+     * This is used in the graph simulation to initialize and update the global ICU occupancy.
+     * @param icu_global The global ICU occupancy time series.
+     */
+    void set_global_icu_occupancy(const mio::TimeSeries<FP>& icu_global)
+    {
+        m_global_icu_occupancy = icu_global;
+    }
+
     /**
      * @brief Transforms the perceived risk into a contact reduction factor and applies it to the contact patterns.
      *
@@ -349,9 +437,9 @@ class FeedbackSimulation
         for (size_t loc = 0; loc < num_locations; ++loc) {
 
             // compute the effective reduction factor using a softplus function.
-            ScalarType reduc_fac = (contactReductionMax[loc] - contactReductionMin[loc]) * epsilon *
-                                       std::log(std::exp(perceived_risk / epsilon) + 1.0) +
-                                   contactReductionMin[loc];
+            FP reduc_fac = (contactReductionMax[loc] - contactReductionMin[loc]) * epsilon *
+                               std::log(std::exp(perceived_risk / epsilon) + 1.0) +
+                           contactReductionMin[loc];
 
             // clamp the reduction factor to the maximum allowed value.
             reduc_fac = std::min(reduc_fac, contactReductionMax[loc]);
@@ -370,7 +458,9 @@ class FeedbackSimulation
     Sim m_simulation; ///< The simulation instance.
     std::vector<size_t> m_icu_indices; ///< The ICU compartment indices from the model.
     FeedbackSimulationParameters<FP> m_feedback_parameters; ///< The feedback parameters.
-    mio::TimeSeries<ScalarType> m_perceived_risk; ///< The perceived risk time series.
+    mio::TimeSeries<FP> m_perceived_risk; ///< The perceived risk time series.
+    mio::TimeSeries<FP> m_regional_icu_occupancy; ///< The regional ICU occupancy time series.
+    mio::TimeSeries<FP> m_global_icu_occupancy; ///< The global ICU occupancy time series.
 };
 
 } // namespace mio

cpp/memilio/compartments/simulation.h

@@ -208,7 +208,7 @@ using is_compartment_model_simulation =
  * @tparam Model The particular Model derived from CompartmentModel to simulate.
  * @tparam Sim A Simulation that can simulate the model.
  */
-template <typename FP, class Model, class Sim = Simulation<FP, Model>>
+template <typename FP, class Model, class Sim = mio::Simulation<FP, Model>>
 TimeSeries<FP> simulate(FP t0, FP tmax, FP dt, Model const& model,
                         std::shared_ptr<IntegratorCore<FP>> integrator = nullptr)
 {