diff --git a/README.md b/README.md index 722e8423cf..5333ab438a 100644 --- a/README.md +++ b/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** diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt index fea58a84cf..980164d324 100644 --- a/cpp/examples/CMakeLists.txt +++ b/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}) diff --git a/cpp/examples/ode_secir_feedback_graph.cpp b/cpp/examples/ode_secir_feedback_graph.cpp new file mode 100644 index 0000000000..1be0c2234e --- /dev/null +++ b/cpp/examples/ode_secir_feedback_graph.cpp @@ -0,0 +1,177 @@ +/* +* Copyright (C) 2020-2025 MEmilio +* +* Authors: Henrik Zunker +* +* Contact: Martin J. Kuehn +* +* 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 + +// alias for the type of the simulation with feedback +using FeedbackSim = mio::FeedbackSimulation>, + mio::osecir::ContactPatterns>; + +// helper function to initialize the model with population and parameters +void initialize_model(mio::osecir::Model& model, double cont_freq) +{ + model.parameters.set(60); + model.parameters.set>(0.2); + + // time-related parameters + model.parameters.get>() = 3.2; + model.parameters.get>() = 2.0; + model.parameters.get>() = 5.8; + model.parameters.get>() = 9.5; + model.parameters.get>() = 7.1; + + // Set transmission and isolation parameters + model.parameters.get>() = 0.05; + model.parameters.get>() = 0.7; + model.parameters.get>() = 0.09; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.45; + model.parameters.get>() = 35; + model.parameters.get>() = 0.2; + model.parameters.get>() = 0.25; + model.parameters.get>() = 0.3; + + // contact matrix + mio::ContactMatrixGroup& contact_matrix = model.parameters.get>(); + 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>() = 10; + + // ICU occupancy in the past for memory kernel + auto& icu_occupancy = feedback_simulation.get_parameters().template get>(); + Eigen::VectorXd icu_day = Eigen::VectorXd::Constant(1, 1); + const auto cutoff = static_cast(feedback_simulation.get_parameters().template get()); + 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>() = {0.1}; + feedback_simulation.get_parameters().template get>() = {0.8}; + + // Set blending factors. The global blending factor is implicitly defined as 1 - local - regional. + feedback_simulation.get_parameters().template get>() = 0.5; + feedback_simulation.get_parameters().template get>() = 0.3; +} + +// helper function to create the graph with nodes and edges +mio::Graph, mio::MobilityEdge> +create_graph(int num_nodes, int total_population, double cont_freq) +{ + // Create a graph for the metapopulation simulation + mio::Graph, mio::MobilityEdge> g; + + // Create models and add nodes to the graph + for (int i = 0; i < num_nodes; ++i) { + mio::osecir::Model 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{ + model.populations.get_flat_index({mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical})}; + + // Create feedback simulation + auto feedback_sim = FeedbackSim(mio::Simulation>(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> 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(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 cols = {"S", "E", "C", "C_confirmed", "I", "I_confirmed", "H", "U", "R", "D"}; + interpolated_results_node.print_table(cols); + + return 0; +} diff --git a/cpp/memilio/compartments/feedback_simulation.h b/cpp/memilio/compartments/feedback_simulation.h index 4f0d035175..87767dce36 100644 --- a/cpp/memilio/compartments/feedback_simulation.h +++ b/cpp/memilio/compartments/feedback_simulation.h @@ -20,6 +20,7 @@ #ifndef FEEDBACK_SIMULATION_H #define FEEDBACK_SIMULATION_H +#include #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 +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 +struct BlendingFactorRegional { + using Type = FP; + static Type get_default(AgeGroup) + { + return FP(0.0); + } + static std::string name() + { + return "BlendingFactorRegional"; + } +}; + template using FeedbackSimulationParameters = ParameterSet, GammaShapeParameter, GammaScaleParameter, GammaCutOff, ContactReductionMax, ContactReductionMin, SoftPlusCurvatureParameter, - NominalICUCapacity>; + NominalICUCapacity, BlendingFactorLocal, BlendingFactorRegional>; /** * @brief A generic feedback simulation extending existing simulations with a feedback mechanism. @@ -181,6 +215,8 @@ template 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(m_simulation.get_model().parameters.get_num_groups())) + , m_perceived_risk(static_cast(m_simulation.get_model().parameters.get_num_groups().get())) + , m_regional_icu_occupancy( + static_cast(m_simulation.get_model().parameters.get_num_groups().get())) + , m_global_icu_occupancy(static_cast(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>(); - size_t num_time_points = icu_occ.get_num_time_points(); + const auto& local_icu_occ = m_feedback_parameters.template get>(); + size_t num_time_points = local_icu_occ.get_num_time_points(); size_t n = std::min(static_cast(num_time_points), m_feedback_parameters.template get()); FP perceived_risk = 0.0; const auto& a = m_feedback_parameters.template get>(); const auto& b = m_feedback_parameters.template get>(); + + const auto& blending_local = m_feedback_parameters.template get>(); + const auto& blending_regional = m_feedback_parameters.template get>(); + 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>(); - 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>(); + } + 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>(); + } + 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>(); + } + + 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>().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& 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& 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 m_icu_indices; ///< The ICU compartment indices from the model. FeedbackSimulationParameters m_feedback_parameters; ///< The feedback parameters. - mio::TimeSeries m_perceived_risk; ///< The perceived risk time series. + mio::TimeSeries m_perceived_risk; ///< The perceived risk time series. + mio::TimeSeries m_regional_icu_occupancy; ///< The regional ICU occupancy time series. + mio::TimeSeries m_global_icu_occupancy; ///< The global ICU occupancy time series. }; } // namespace mio diff --git a/cpp/memilio/compartments/simulation.h b/cpp/memilio/compartments/simulation.h index fc9cee4435..a00dce23e1 100644 --- a/cpp/memilio/compartments/simulation.h +++ b/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 > +template > TimeSeries simulate(FP t0, FP tmax, FP dt, Model const& model, std::shared_ptr> integrator = nullptr) { diff --git a/cpp/memilio/mobility/graph_simulation.h b/cpp/memilio/mobility/graph_simulation.h index b911b43fbd..cc42c4ecf7 100644 --- a/cpp/memilio/mobility/graph_simulation.h +++ b/cpp/memilio/mobility/graph_simulation.h @@ -1,7 +1,7 @@ /* * Copyright (C) 2020-2025 MEmilio * -* Authors: Daniel Abele +* Authors: Daniel Abele, Henrik Zunker * * Contact: Martin J. Kuehn * @@ -22,6 +22,8 @@ #include "memilio/mobility/graph.h" #include "memilio/utils/random_number_generator.h" +#include "memilio/compartments/feedback_simulation.h" +#include "memilio/geography/regions.h" namespace mio { @@ -266,5 +268,206 @@ auto make_graph_sim_stochastic(FP t0, FP dt, Graph&& g, NodeF&& node_func, EdgeF t0, dt, std::forward(g), std::forward(node_func), std::forward(edge_func)); } +// FeedbackGraphSimulation is only allowed to be used with local FeedbackSimulation. +// Therefore, we use type traits to check if the type is a specialization of FeedbackSimulation +template +struct is_feedback_simulation : std::false_type { +}; + +template +struct is_feedback_simulation> : std::true_type { +}; + +template +class FeedbackGraphSimulation +{ +public: + FeedbackGraphSimulation(Graph&& g, FP t0, FP dt) + : m_graph(std::move(g)) + , m_t(t0) + , m_dt(dt) + , m_initialized(false) + , m_global_icu_occupancy(6) + { + using SimT = decltype(m_graph.nodes()[0].property.get_simulation()); + static_assert(is_feedback_simulation>::value, + "Graph node simulation must be a FeedbackSimulation."); + } + + void advance(FP t_max) + { + // Initialize global and regional ICU occupancy if not done yet + if (!m_initialized) { + calculate_global_icu_occupancy(); + calculate_regional_icu_occupancy(); + m_initialized = true; + } + + while (m_t < t_max) { + FP dt_eff = std::min(m_dt, t_max - m_t); + + // assign the regional and global ICU occupancy to each node + distribute_icu_data(); + + // apply feedback and advance simulation for each node + for (auto& node : m_graph.nodes()) { + node.property.get_simulation().advance(m_t + dt_eff, m_dt); + } + + // apply mobility between nodes + for (auto& edge : m_graph.edges()) { + auto& node_from = m_graph.nodes()[edge.start_node_idx]; + auto& node_to = m_graph.nodes()[edge.end_node_idx]; + edge.property.apply_mobility(m_t, m_dt, node_from.property, node_to.property); + } + + // update ICU occupancy for each node + update_global_icu_occupancy(); + update_regional_icu_occupancy(); + + m_t += dt_eff; + } + } + + Graph& get_graph() + { + return m_graph; + } + +private: + /** + * @brief Calculates the global ICU occupancy based on the ICU history of all nodes. + * The result is stored in m_global_icu_occupancy. + */ + void calculate_global_icu_occupancy() + { + auto& first_node_sim = m_graph.nodes()[0].property.get_simulation(); + auto& first_node_icu = first_node_sim.get_parameters().template get>(); + m_global_icu_occupancy = mio::TimeSeries(first_node_icu.get_num_elements()); + m_global_icu_occupancy.add_time_point(0.0, Eigen::VectorXd::Zero(first_node_icu.get_num_elements())); + + FP total_population = 0; + for (Eigen::Index i = 0; i < first_node_icu.get_num_time_points(); ++i) { + Eigen::VectorXd sum = Eigen::VectorXd::Zero(first_node_icu.get_num_elements()); + for (auto& node : m_graph.nodes()) { + auto& sim = node.property.get_simulation(); + auto& icu_history = sim.get_parameters().template get>(); + sum += icu_history.get_value(i) * node.property.get_simulation().get_model().populations.get_total(); + total_population += node.property.get_simulation().get_model().populations.get_total(); + } + m_global_icu_occupancy.add_time_point(first_node_icu.get_time(i), sum / total_population); + } + } + + /** + * @brief Calculates the regional ICU occupancy for each region. + * The calculation is based on the ICU history of the nodes in each region. + * The results are stored in m_regional_icu_occupancy. + */ + void calculate_regional_icu_occupancy() + { + std::unordered_map regional_population; + for (auto& node : m_graph.nodes()) { + auto region_id = mio::regions::get_state_id(node.id).get(); + regional_population[region_id] += node.property.get_simulation().get_model().populations.get_total(); + } + + auto& first_node_sim = m_graph.nodes()[0].property.get_simulation(); + auto& first_node_icu = first_node_sim.get_parameters().template get>(); + Eigen::Index num_time_points = first_node_icu.get_num_time_points(); + Eigen::Index num_elements = first_node_icu.get_num_elements(); + + // For each region: vector of summed ICU values for all time points + std::unordered_map> regional_sums; + for (auto const& [region_id, _] : regional_population) { + regional_sums[region_id] = + std::vector(num_time_points, Eigen::VectorXd::Zero(num_elements)); + } + + // Sum up + for (auto& node : m_graph.nodes()) { + auto region_id = mio::regions::get_state_id(node.id).get(); + auto& sim = node.property.get_simulation(); + auto& icu_history = sim.get_parameters().template get>(); + FP pop = node.property.get_simulation().get_model().populations.get_total(); + for (Eigen::Index i = 0; i < num_time_points; ++i) { + regional_sums[region_id][i] += icu_history.get_value(i) * pop; + } + } + + // Initialize TimeSeries and insert values + m_regional_icu_occupancy.clear(); + for (auto const& [region_id, sum_vec] : regional_sums) { + mio::TimeSeries ts(num_elements); + for (Eigen::Index i = 0; i < num_time_points; ++i) { + ts.add_time_point(first_node_icu.get_time(i), sum_vec[i] / regional_population[region_id]); + } + m_regional_icu_occupancy.emplace(region_id, std::move(ts)); + } + } + + /** + * @brief Updates the global ICU occupancy with the latest values from all nodes. + * A new time point is added to m_global_icu_occupancy. + */ + void update_global_icu_occupancy() + { + FP total_population = 0; + Eigen::VectorXd sum = Eigen::VectorXd::Zero(m_global_icu_occupancy.get_num_elements()); + for (auto& node : m_graph.nodes()) { + auto& sim = node.property.get_simulation(); + auto& icu_history = sim.get_parameters().template get>(); + sum += icu_history.get_last_value() * node.property.get_simulation().get_model().populations.get_total(); + total_population += node.property.get_simulation().get_model().populations.get_total(); + } + m_global_icu_occupancy.add_time_point(m_t, sum / total_population); + } + + /** + * @brief Updates the regional ICU occupancy for each region with the latest values. + * A new time point is added to each TimeSeries in m_regional_icu_occupancy. + */ + void update_regional_icu_occupancy() + { + + std::unordered_map regional_population; + for (auto& [region_id, regional_data] : m_regional_icu_occupancy) { + Eigen::VectorXd sum = Eigen::VectorXd::Zero(regional_data.get_num_elements()); + for (auto& node : m_graph.nodes()) { + if (mio::regions::get_state_id(node.id).get() == region_id) { + auto& sim = node.property.get_simulation(); + auto& icu_history = sim.get_parameters().template get>(); + sum += icu_history.get_last_value() * + node.property.get_simulation().get_model().populations.get_total(); + regional_population[region_id] += + node.property.get_simulation().get_model().populations.get_total(); + } + } + regional_data.add_time_point(m_t, sum / regional_population[region_id]); + } + } + + /** + * @brief Distributes the calculated global and regional ICU occupancy data to each node. + * This makes the data available for feedback mechanisms within the node-local simulations. + */ + void distribute_icu_data() + { + for (auto& node : m_graph.nodes()) { + auto region_id = mio::regions::get_state_id(node.id).get(); + auto& sim = node.property.get_simulation(); + sim.set_regional_icu_occupancy(m_regional_icu_occupancy.at(region_id)); + sim.set_global_icu_occupancy(m_global_icu_occupancy); + } + } + + Graph m_graph; + FP m_t; + FP m_dt; + bool m_initialized; + mio::TimeSeries m_global_icu_occupancy; + std::unordered_map> m_regional_icu_occupancy; +}; + } // namespace mio #endif //MIO_MOBILITY_GRAPH_SIMULATION_H diff --git a/cpp/tests/test_feedback.cpp b/cpp/tests/test_feedback.cpp index a3adca9574..df41d78df3 100644 --- a/cpp/tests/test_feedback.cpp +++ b/cpp/tests/test_feedback.cpp @@ -108,6 +108,48 @@ TEST_F(TestFeedbackSimulation, CalcPerceivedRisk) EXPECT_NEAR(risk, 0.27357588823428847, 1e-10); } +TEST_F(TestFeedbackSimulation, CalcPerceivedRiskWithBlending) +{ + auto& fb_params = feedback_sim->get_parameters(); + fb_params.template get>() = 1; + fb_params.template get>() = 1; + fb_params.template get>() = 10; + + // add a single local ICU occupancy time point with value 2. + Eigen::VectorXd icu_value_local(1); + icu_value_local << 2; + auto& icu_occ_local = fb_params.template get>(); + icu_occ_local.add_time_point(0.0, icu_value_local); + + // create regional and global ICU occupancy data + mio::TimeSeries icu_occ_regional(1), icu_occ_global(1); + Eigen::VectorXd icu_value_regional(1), icu_value_global(1); + icu_value_regional << 4; + icu_value_global << 6; + icu_occ_regional.add_time_point(0.0, icu_value_regional); + icu_occ_global.add_time_point(0.0, icu_value_global); + + feedback_sim->set_regional_icu_occupancy(icu_occ_regional); + feedback_sim->set_global_icu_occupancy(icu_occ_global); + + // Scenario 1: Test with local, regional, and global blending + fb_params.template get>() = 0.2; + fb_params.template get>() = 0.3; + // global blending factor is 1 - 0.2 - 0.3 = 0.5 + + double risk = feedback_sim->calc_risk_perceived(); + // Expected risk: 0.2 * (2/10) + 0.3 * (4/10) + 0.5 * (6/10) = 0.46 + EXPECT_NEAR(risk, 0.46, 1e-10); + + // Scenario 2: Test with only regional blending + fb_params.template get>() = 0.0; + fb_params.template get>() = 1.0; + + risk = feedback_sim->calc_risk_perceived(); + // Expected risk: 1.0 * (4/10) = 0.4 + EXPECT_NEAR(risk, 0.4, 1e-10); +} + TEST_F(TestFeedbackSimulation, ApplyFeedback) { // bounds for contact reduction measures diff --git a/cpp/tests/test_graph_simulation.cpp b/cpp/tests/test_graph_simulation.cpp index 485373b285..46c7074e8c 100644 --- a/cpp/tests/test_graph_simulation.cpp +++ b/cpp/tests/test_graph_simulation.cpp @@ -23,6 +23,7 @@ #include "memilio/mobility/metapopulation_mobility_instant.h" #include "memilio/mobility/metapopulation_mobility_stochastic.h" #include "memilio/compartments/simulation.h" +#include "memilio/compartments/feedback_simulation.h" #include "ode_seir/model.h" #include "gtest/gtest.h" #include "load_test_data.h" @@ -120,8 +121,8 @@ TEST(TestGraphSimulation, stopsAtTmax) const auto tmax = 3.123; const auto dt = 0.076; - auto sim = - mio::make_graph_sim(t0, dt, g, [](auto&&, auto&&, auto&&) {}, [](auto&&, auto&&, auto&&, auto&&, auto&&) {}); + auto sim = mio::make_graph_sim( + t0, dt, g, [](auto&&, auto&&, auto&&) {}, [](auto&&, auto&&, auto&&, auto&&, auto&&) {}); sim.advance(tmax); @@ -333,15 +334,82 @@ TEST(TestGraphSimulation, consistencyFlowMobility) } } +TEST(TestGraphSimulation, feedbackSimulation) +{ + using Model = mio::oseir::Model; + using Simulation = mio::Simulation; + using FeedbackSim = mio::FeedbackSimulation>; + using Node = mio::SimulationNode; + using Edge = mio::MobilityEdge; + using Graph = mio::Graph; + using GraphFeedback = mio::FeedbackGraphSimulation; + + double t0 = 0; + double tmax = 5.0; + double dt = 1.0; + + Model model(1); + model.populations.set_total(1000); + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Susceptible}] = 900; + model.populations[{mio::AgeGroup(0), mio::oseir::InfectionState::Exposed}] = 100; + + std::vector icu_indices = {(size_t)mio::oseir::InfectionState::Infected}; + + Graph g; + + const auto num_nodes = 2; + for (int i = 0; i < num_nodes; ++i) { + + auto feedback_sim = FeedbackSim(Simulation(model, t0), icu_indices); + + // set feedback parameters + feedback_sim.get_parameters().get>() = 10; + auto& icu_occupancy = feedback_sim.get_parameters().get>(); + Eigen::VectorXd icu_day = Eigen::VectorXd::Constant(1, 1); + const auto cutoff = static_cast(feedback_sim.get_parameters().get()); + for (int t = -cutoff; t <= 0; ++t) { + icu_occupancy.add_time_point(t, icu_day); + } + + // bounds for contact reduction measures + feedback_sim.get_parameters().get>() = {0.1}; + feedback_sim.get_parameters().get>() = {0.8}; + + // Set blending factors. The global blending factor is implicitly defined as 1 - local - regional. + feedback_sim.get_parameters().get>() = 0.5; + feedback_sim.get_parameters().get>() = 0.3; + + g.add_node(i, std::move(feedback_sim)); + } + g.add_edge(0, 1, Eigen::VectorXd::Constant(4, 0.01)); + + double total_pop_before = 0; + for (auto& node : g.nodes()) { + total_pop_before += node.property.get_simulation().get_model().populations.get_total(); + } + + GraphFeedback sim(std::move(g), t0, dt); + + sim.advance(tmax); + + double total_pop_after = 0; + for (auto& node : sim.get_graph().nodes()) { + total_pop_after += node.property.get_simulation().get_model().populations.get_total(); + } + + EXPECT_NEAR(total_pop_before, total_pop_after, 1e-10); + EXPECT_NEAR(sim.get_graph().nodes()[0].property.get_simulation().get_result().get_last_time(), tmax, 1e-10); +} + namespace { struct MoveOnly { MoveOnly(); - MoveOnly(const MoveOnly&) = delete; + MoveOnly(const MoveOnly&) = delete; MoveOnly& operator=(const MoveOnly&) = delete; MoveOnly(MoveOnly&&) = default; - MoveOnly& operator=(MoveOnly&&) = default; + MoveOnly& operator=(MoveOnly&&) = default; }; using MoveOnlyGraph = mio::Graph; using MoveOnlyGraphSim = mio::GraphSimulation; diff --git a/docs/source/cpp/graph_metapop.rst b/docs/source/cpp/graph_metapop.rst index b5498fe946..eaa57fa1c3 100644 --- a/docs/source/cpp/graph_metapop.rst +++ b/docs/source/cpp/graph_metapop.rst @@ -43,6 +43,17 @@ For further details, please refer to: In the stochastic mobility approach, transitions of individuals between regions, i.e. graph nodes, are modeled as stochastic jumps. The frequency of individuals transitioning from region i to region j is determined by a rate :math:`\lambda_{ij}` which can be interpreted as the (expected) fraction of the population in region i that commutes to region j within one day. In contrast to the instant and detailed mobility approach, the time points and the number of transitions between two regions are stochastic. +Behavior-based ODE models +------------------------- + +The graph-based simulation can be combined with the feedback mechanism described in the :doc:`ODE models documentation `. +This allows for modeling behavioral changes in response to the epidemiological situation not only on a local level but also considering regional and global information. The ``FeedbackGraphSimulation`` extends the local feedback by incorporating a weighted blend of local, regional, and global ICU occupancy to calculate the perceived risk. This blended risk then influences the contact patterns in each node, allowing for a more nuanced and realistic simulation of public health responses across a metapopulation. + +The extension and inclusion of regional and global information is based on the following paper: + +- 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` + + How to: Set up a graph and run a graph simulation ------------------------------------------------- diff --git a/docs/source/cpp/ode.rst b/docs/source/cpp/ode.rst index ee435716c7..c7a0cf57a9 100644 --- a/docs/source/cpp/ode.rst +++ b/docs/source/cpp/ode.rst @@ -110,6 +110,8 @@ You can run a simulation using either fixed or adaptive solution schemes with an default, the simulation uses an adaptive solution scheme of the boost library and an absolute tolerance of 1e-10 and a relative tolerance of 1e-5. For more details on the possible integration schemes, see . +Additionally, a feedback simulation is available, which allows for dynamic adjustments of the simulation based on its own output. This is particularly useful for modeling behavioral changes in response to the epidemiological situation. The local feedback simulation is based on the work of Dönges et al. (doi.org/10.3389/fphy.2022.842180). It uses the history of ICU occupancy to calculate a "perceived risk", which then translates into a reduction of contact rates. The ``FeedbackSimulation`` class wraps an existing simulation and applies this feedback mechanism at regular intervals. This concept can be extended to a metapopulation model using a ``FeedbackGraphSimulation``, which is described in more detail in the "Graph-based metapopulation model" section. + Output ------ @@ -142,5 +144,3 @@ List of models models/osecir models/osecirvvs models/osecirts - - \ No newline at end of file