Skip to content

Commit 9d72115

Browse files
HenrZujubicker
HenrZu
and
jubicker
authored
1310 Spatially resolved feedback simulation (#1311)
- As used in https://doi.org/10.1016/j.chaos.2025.116782, we extend the already implemented local feedback mechanism (see Issue 978) - Now, we can consider regional and global information for risk perception - Added documentation for readthedocs Co-authored-by: jubicker <[email protected]>
1 parent 2defe8b commit 9d72115

10 files changed

Lines changed: 615 additions & 19 deletions

File tree

‎README.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,8 @@ and, in particular, for
1919
- 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`
2020
- 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`
2121
- 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`
22-
- 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`
22+
- 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`
23+
2324

2425
**Getting started**
2526

‎cpp/examples/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,10 @@ add_executable(ode_secir_feedback ode_secir_feedback.cpp)
8080
target_link_libraries(ode_secir_feedback PRIVATE memilio ode_secir)
8181
target_compile_options(ode_secir_feedback PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
8282

83+
add_executable(ode_secir_feedback_graph ode_secir_feedback_graph.cpp)
84+
target_link_libraries(ode_secir_feedback_graph PRIVATE memilio ode_secir)
85+
target_compile_options(ode_secir_feedback_graph PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
86+
8387
add_executable(ode_secirvvs_example ode_secirvvs.cpp)
8488
target_link_libraries(ode_secirvvs_example PRIVATE memilio ode_secirvvs)
8589
target_compile_options(ode_secirvvs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/ode_secir_feedback_graph.cpp‎

Lines changed: 177 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,177 @@
1+
/*
2+
* Copyright (C) 2020-2025 MEmilio
3+
*
4+
* Authors: Henrik Zunker
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
7+
*
8+
* Licensed under the Apache License, Version 2.0 (the "License");
9+
* you may not use this file except in compliance with the License.
10+
* You may obtain a copy of the License at
11+
*
12+
* http://www.apache.org/licenses/LICENSE-2.0
13+
*
14+
* Unless required by applicable law or agreed to in writing, software
15+
* distributed under the License is distributed on an "AS IS" BASIS,
16+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17+
* See the License for the specific language governing permissions and
18+
* limitations under the License.
19+
*/
20+
#include "memilio/data/analyze_result.h"
21+
#include "ode_secir/model.h"
22+
#include "memilio/compartments/feedback_simulation.h"
23+
#include "memilio/mobility/metapopulation_mobility_instant.h"
24+
#include "memilio/mobility/graph_simulation.h"
25+
#include "memilio/utils/logging.h"
26+
#include <iostream>
27+
28+
// alias for the type of the simulation with feedback
29+
using FeedbackSim = mio::FeedbackSimulation<double, mio::Simulation<double, mio::osecir::Model<double>>,
30+
mio::osecir::ContactPatterns<double>>;
31+
32+
// helper function to initialize the model with population and parameters
33+
void initialize_model(mio::osecir::Model<double>& model, double cont_freq)
34+
{
35+
model.parameters.set<mio::osecir::StartDay>(60);
36+
model.parameters.set<mio::osecir::Seasonality<double>>(0.2);
37+
38+
// Mean stay times per compartment
39+
model.parameters.get<mio::osecir::TimeExposed<double>>() = 3.2;
40+
model.parameters.get<mio::osecir::TimeInfectedNoSymptoms<double>>() = 2.0;
41+
model.parameters.get<mio::osecir::TimeInfectedSymptoms<double>>() = 5.8;
42+
model.parameters.get<mio::osecir::TimeInfectedSevere<double>>() = 9.5;
43+
model.parameters.get<mio::osecir::TimeInfectedCritical<double>>() = 7.1;
44+
45+
// Set transmission and isolation parameters
46+
model.parameters.get<mio::osecir::TransmissionProbabilityOnContact<double>>() = 0.05;
47+
model.parameters.get<mio::osecir::RelativeTransmissionNoSymptoms<double>>() = 0.7;
48+
model.parameters.get<mio::osecir::RecoveredPerInfectedNoSymptoms<double>>() = 0.09;
49+
model.parameters.get<mio::osecir::RiskOfInfectionFromSymptomatic<double>>() = 0.25;
50+
model.parameters.get<mio::osecir::MaxRiskOfInfectionFromSymptomatic<double>>() = 0.45;
51+
model.parameters.get<mio::osecir::TestAndTraceCapacity<double>>() = 35;
52+
model.parameters.get<mio::osecir::SeverePerInfectedSymptoms<double>>() = 0.2;
53+
model.parameters.get<mio::osecir::CriticalPerSevere<double>>() = 0.25;
54+
model.parameters.get<mio::osecir::DeathsPerCritical<double>>() = 0.3;
55+
56+
// contact matrix
57+
mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::osecir::ContactPatterns<double>>();
58+
contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, cont_freq));
59+
}
60+
61+
// helper function to initialize the feedback mechanism parameters for a simulation
62+
void initialize_feedback(FeedbackSim& feedback_simulation)
63+
{
64+
// nominal ICU capacity
65+
feedback_simulation.get_parameters().template get<mio::NominalICUCapacity<double>>() = 10;
66+
67+
// ICU occupancy in the past for memory kernel
68+
auto& icu_occupancy = feedback_simulation.get_parameters().template get<mio::ICUOccupancyHistory<double>>();
69+
Eigen::VectorXd icu_day = Eigen::VectorXd::Constant(1, 1);
70+
const auto cutoff = static_cast<int>(feedback_simulation.get_parameters().template get<mio::GammaCutOff>());
71+
for (int t = -cutoff; t <= 0; ++t) {
72+
icu_occupancy.add_time_point(t, icu_day);
73+
}
74+
75+
// bounds for contact reduction measures
76+
feedback_simulation.get_parameters().template get<mio::ContactReductionMin<double>>() = {0.1};
77+
feedback_simulation.get_parameters().template get<mio::ContactReductionMax<double>>() = {0.8};
78+
79+
// Set blending factors. The global blending factor is implicitly defined as 1 - local - regional.
80+
feedback_simulation.get_parameters().template get<mio::BlendingFactorLocal<double>>() = 0.5;
81+
feedback_simulation.get_parameters().template get<mio::BlendingFactorRegional<double>>() = 0.3;
82+
}
83+
84+
// helper function to create the graph with nodes and edges
85+
mio::Graph<mio::SimulationNode<FeedbackSim>, mio::MobilityEdge<double>>
86+
create_graph(int num_nodes, int total_population, double cont_freq)
87+
{
88+
// Create a graph for the metapopulation simulation
89+
mio::Graph<mio::SimulationNode<FeedbackSim>, mio::MobilityEdge<double>> g;
90+
91+
// Create models and add nodes to the graph
92+
for (int i = 0; i < num_nodes; ++i) {
93+
mio::osecir::Model<double> model(1);
94+
initialize_model(model, cont_freq);
95+
96+
// Set initial populations (infection starts in the first node)
97+
if (i == 0) {
98+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Exposed}] = total_population * 0.1;
99+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedNoSymptoms}] =
100+
total_population * 0.1;
101+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSymptoms}] =
102+
total_population * 0.05;
103+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedSevere}] =
104+
total_population * 0.02;
105+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical}] =
106+
total_population * 0.01;
107+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Recovered}] = 0;
108+
model.populations.set_difference_from_total({mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible},
109+
total_population);
110+
}
111+
else {
112+
model.populations[{mio::AgeGroup(0), mio::osecir::InfectionState::Susceptible}] = total_population;
113+
}
114+
model.apply_constraints();
115+
116+
// Determine the index for the ICU state (InfectedCritical) for the feedback mechanism
117+
auto icu_index = std::vector<size_t>{
118+
model.populations.get_flat_index({mio::AgeGroup(0), mio::osecir::InfectionState::InfectedCritical})};
119+
120+
// Create feedback simulation
121+
auto feedback_sim = FeedbackSim(mio::Simulation<double, mio::osecir::Model<double>>(model), icu_index);
122+
initialize_feedback(feedback_sim);
123+
124+
// Node-ID-Logic: 1001-1005, 2001-2005, ...
125+
const int region_id = i / 5;
126+
const int local_id = i % 5;
127+
const int node_id = (region_id + 1) * 1000 + (local_id + 1);
128+
g.add_node(node_id, std::move(feedback_sim));
129+
}
130+
131+
// Define complete graph, i.e. each node is connected to every other node
132+
std::vector<std::vector<size_t>> mobile_compartments(2);
133+
for (size_t i = 0; i < g.nodes().size(); ++i) {
134+
for (size_t j = 0; j < g.nodes().size(); ++j) {
135+
if (i != j) {
136+
g.add_edge(i, j, Eigen::VectorXd::Constant((size_t)mio::osecir::InfectionState::Count, 0.01));
137+
}
138+
}
139+
}
140+
141+
return g;
142+
}
143+
144+
int main()
145+
{
146+
// This example demonstrates the implementation of a feedback mechanism for a ODE SECIR model in a graph.
147+
// It shows how the perceived risk dynamically impacts contact reduction measures in different regions (nodes).
148+
mio::set_log_level(mio::LogLevel::err);
149+
150+
const auto t0 = 0.;
151+
const auto tmax = 10.;
152+
const auto dt = 0.5;
153+
const int total_population = 1000;
154+
const double cont_freq = 2.7;
155+
const int num_nodes = 10;
156+
157+
// Create the graph
158+
auto g = create_graph(num_nodes, total_population, cont_freq);
159+
160+
// Create and run the simulation
161+
using Graph = decltype(g);
162+
auto sim = mio::FeedbackGraphSimulation<double, Graph>(std::move(g), t0, dt);
163+
sim.advance(tmax);
164+
165+
// The output shows the compartments sizes for a node without any initial infections.
166+
auto& node = sim.get_graph().nodes()[1];
167+
auto& results_node = node.property.get_simulation().get_result();
168+
// interpolate results
169+
auto interpolated_results_node = mio::interpolate_simulation_result(results_node);
170+
171+
// print result with print_table
172+
std::cout << "Node ID: " << node.id << "\n";
173+
std::vector<std::string> cols = {"S", "E", "C", "C_confirmed", "I", "I_confirmed", "H", "U", "R", "D"};
174+
interpolated_results_node.print_table(cols);
175+
176+
return 0;
177+
}

‎cpp/memilio/compartments/feedback_simulation.h‎

Lines changed: 100 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
#ifndef MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H
2121
#define MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H
2222

23+
#include <cassert>
2324
#include "memilio/compartments/simulation.h"
2425
#include "memilio/utils/time_series.h"
2526
#include "memilio/utils/parameter_set.h"
@@ -160,11 +161,44 @@ struct NominalICUCapacity {
160161
}
161162
};
162163

164+
/**
165+
* @brief Blending factor for local ICU occupancy.
166+
* If we only have a local simulation, this factor is set to 1.0 per default.
167+
*/
168+
template <typename FP>
169+
struct BlendingFactorLocal {
170+
using Type = FP;
171+
static Type get_default(AgeGroup)
172+
{
173+
return FP(1.0);
174+
}
175+
static std::string name()
176+
{
177+
return "BlendingFactorLocal";
178+
}
179+
};
180+
181+
/**
182+
* @brief Blending factor for regional ICU occupancy.
183+
*/
184+
template <typename FP>
185+
struct BlendingFactorRegional {
186+
using Type = FP;
187+
static Type get_default(AgeGroup)
188+
{
189+
return FP(0.0);
190+
}
191+
static std::string name()
192+
{
193+
return "BlendingFactorRegional";
194+
}
195+
};
196+
163197
template <typename FP>
164198
using FeedbackSimulationParameters =
165199
ParameterSet<ICUOccupancyHistory<FP>, GammaShapeParameter<FP>, GammaScaleParameter<FP>, GammaCutOff,
166200
ContactReductionMax<FP>, ContactReductionMin<FP>, SoftPlusCurvatureParameter<FP>,
167-
NominalICUCapacity<FP>>;
201+
NominalICUCapacity<FP>, BlendingFactorLocal<FP>, BlendingFactorRegional<FP>>;
168202

169203
/**
170204
* @brief A generic feedback simulation extending existing simulations with a feedback mechanism.
@@ -181,6 +215,8 @@ template <typename FP, typename Sim, typename ContactPatterns>
181215
class FeedbackSimulation
182216
{
183217
public:
218+
using Model = typename Sim::Model;
219+
184220
/**
185221
* @brief Constructs the FeedbackSimulation by taking ownership of an existing simulation instance.
186222
*
@@ -191,7 +227,10 @@ class FeedbackSimulation
191227
: m_simulation(std::move(sim))
192228
, m_icu_indices(icu_indices)
193229
, m_feedback_parameters(m_simulation.get_model().parameters.get_num_groups())
194-
, m_perceived_risk(static_cast<size_t>(m_simulation.get_model().parameters.get_num_groups()))
230+
, m_perceived_risk(static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
231+
, m_regional_icu_occupancy(
232+
static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
233+
, m_global_icu_occupancy(static_cast<Eigen::Index>(m_simulation.get_model().parameters.get_num_groups().get()))
195234
{
196235
}
197236

@@ -238,6 +277,11 @@ class FeedbackSimulation
238277
return m_simulation.get_model();
239278
}
240279

280+
const auto& get_model() const
281+
{
282+
return m_simulation.get_model();
283+
}
284+
241285
/**
242286
* @brief Returns the perceived risk time series.
243287
*/
@@ -268,17 +312,40 @@ class FeedbackSimulation
268312
*/
269313
FP calc_risk_perceived()
270314
{
271-
const auto& icu_occ = m_feedback_parameters.template get<ICUOccupancyHistory<FP>>();
272-
size_t num_time_points = icu_occ.get_num_time_points();
315+
const auto& local_icu_occ = m_feedback_parameters.template get<ICUOccupancyHistory<FP>>();
316+
size_t num_time_points = local_icu_occ.get_num_time_points();
273317
size_t n = std::min(static_cast<size_t>(num_time_points), m_feedback_parameters.template get<GammaCutOff>());
274318
FP perceived_risk = 0.0;
275319
const auto& a = m_feedback_parameters.template get<GammaShapeParameter<FP>>();
276320
const auto& b = m_feedback_parameters.template get<GammaScaleParameter<FP>>();
321+
322+
const auto& blending_local = m_feedback_parameters.template get<BlendingFactorLocal<FP>>();
323+
const auto& blending_regional = m_feedback_parameters.template get<BlendingFactorRegional<FP>>();
324+
const auto& blending_global = 1 - blending_local - blending_regional;
325+
326+
// assert that that each blending factor is positive
327+
assert(blending_local >= 0 && blending_regional >= 0 && blending_global >= 0 &&
328+
"Blending factors must be non-negative.");
329+
277330
for (size_t i = num_time_points - n; i < num_time_points; ++i) {
278331
size_t day = i - (num_time_points - n);
279332
FP gamma = std::pow(b, a) * std::pow(day, a - 1) * std::exp(-b * day) / std::tgamma(a);
280-
FP perc_risk = icu_occ.get_value(i).sum() / m_feedback_parameters.template get<NominalICUCapacity<FP>>();
281-
perc_risk = std::min(perc_risk, 1.0);
333+
FP perc_risk = 0.0;
334+
335+
if (blending_local > 0) {
336+
perc_risk += blending_local * local_icu_occ.get_value(i).sum() /
337+
m_feedback_parameters.template get<NominalICUCapacity<FP>>();
338+
}
339+
if (blending_regional > 0 && m_regional_icu_occupancy.get_num_time_points() > 0) {
340+
perc_risk += blending_regional * m_regional_icu_occupancy.get_value(i).sum() /
341+
m_feedback_parameters.template get<NominalICUCapacity<FP>>();
342+
}
343+
if (blending_global > 0 && m_global_icu_occupancy.get_num_time_points() > 0) {
344+
perc_risk += blending_global * m_global_icu_occupancy.get_value(i).sum() /
345+
m_feedback_parameters.template get<NominalICUCapacity<FP>>();
346+
}
347+
348+
perc_risk = std::min(perc_risk, 1.0);
282349
perceived_risk += perc_risk * gamma;
283350
}
284351
return perceived_risk;
@@ -307,6 +374,27 @@ class FeedbackSimulation
307374
m_feedback_parameters.template get<ICUOccupancyHistory<FP>>().add_time_point(t, icu_occ);
308375
}
309376

377+
/**
378+
* @brief Sets the regional ICU occupancy time series.
379+
* This is used in the graph simulation to initialize and update the regional ICU occupancy.
380+
*
381+
* @param icu_regional The regional ICU occupancy time series.
382+
*/
383+
void set_regional_icu_occupancy(const mio::TimeSeries<FP>& icu_regional)
384+
{
385+
m_regional_icu_occupancy = icu_regional;
386+
}
387+
388+
/**
389+
* @brief Sets the global ICU occupancy time series.
390+
* This is used in the graph simulation to initialize and update the global ICU occupancy.
391+
* @param icu_global The global ICU occupancy time series.
392+
*/
393+
void set_global_icu_occupancy(const mio::TimeSeries<FP>& icu_global)
394+
{
395+
m_global_icu_occupancy = icu_global;
396+
}
397+
310398
/**
311399
* @brief Transforms the perceived risk into a contact reduction factor and applies it to the contact patterns.
312400
*
@@ -349,9 +437,9 @@ class FeedbackSimulation
349437
for (size_t loc = 0; loc < num_locations; ++loc) {
350438

351439
// compute the effective reduction factor using a softplus function.
352-
ScalarType reduc_fac = (contactReductionMax[loc] - contactReductionMin[loc]) * epsilon *
353-
std::log(std::exp(perceived_risk / epsilon) + 1.0) +
354-
contactReductionMin[loc];
440+
FP reduc_fac = (contactReductionMax[loc] - contactReductionMin[loc]) * epsilon *
441+
std::log(std::exp(perceived_risk / epsilon) + 1.0) +
442+
contactReductionMin[loc];
355443

356444
// clamp the reduction factor to the maximum allowed value.
357445
reduc_fac = std::min(reduc_fac, contactReductionMax[loc]);
@@ -370,7 +458,9 @@ class FeedbackSimulation
370458
Sim m_simulation; ///< The simulation instance.
371459
std::vector<size_t> m_icu_indices; ///< The ICU compartment indices from the model.
372460
FeedbackSimulationParameters<FP> m_feedback_parameters; ///< The feedback parameters.
373-
mio::TimeSeries<ScalarType> m_perceived_risk; ///< The perceived risk time series.
461+
mio::TimeSeries<FP> m_perceived_risk; ///< The perceived risk time series.
462+
mio::TimeSeries<FP> m_regional_icu_occupancy; ///< The regional ICU occupancy time series.
463+
mio::TimeSeries<FP> m_global_icu_occupancy; ///< The global ICU occupancy time series.
374464
};
375465

376466
} // namespace mio

‎cpp/memilio/compartments/simulation.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -109,7 +109,7 @@ using is_compartment_model_simulation =
109109
* @tparam Model The particular Model derived from CompartmentModel to simulate.
110110
* @tparam Sim A Simulation that can simulate the model.
111111
*/
112-
template <typename FP, class Model, class Sim = Simulation<FP, Model>>
112+
template <typename FP, class Model, class Sim = mio::Simulation<FP, Model>>
113113
TimeSeries<FP> simulate(FP t0, FP tmax, FP dt, Model const& model,
114114
std::shared_ptr<OdeIntegratorCore<FP>> integrator = nullptr)
115115
{

0 commit comments

Comments
 (0)