Skip to content

Commit

Permalink
323 SECIRTS model with waning and temporary immunity (#963)
Browse files Browse the repository at this point in the history
- Add new ODE model including waning immunity and temporary immunity states build upon the secirvvs model
- New initialization method
- Allow to read booster vaccinations as optional args from vaccination data in cpp

Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com>
Co-authored-by: Martin J. Kühn <62713180+mknaranja@users.noreply.github.com>
  • Loading branch information
3 people authored Dec 10, 2024
1 parent 3a06d82 commit f9b7610
Show file tree
Hide file tree
Showing 41 changed files with 24,661 additions and 17,865 deletions.
1 change: 1 addition & 0 deletions cpp/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,7 @@ add_subdirectory(memilio)
if(MEMILIO_BUILD_MODELS)
add_subdirectory(models/abm)
add_subdirectory(models/ode_secir)
add_subdirectory(models/ode_secirts)
add_subdirectory(models/ode_secirvvs)
add_subdirectory(models/lct_secir)
add_subdirectory(models/glct_secir)
Expand Down
25 changes: 14 additions & 11 deletions cpp/benchmarks/flow_simulation_ode_secirvvs.h
Original file line number Diff line number Diff line change
Expand Up @@ -486,15 +486,17 @@ class Simulation : public Base
double first_vacc;
double full_vacc;
if (t_idx == SimulationDay(0)) {
first_vacc = params.template get<osecirvvs::DailyFirstVaccination<ScalarType>>()[{(AgeGroup)i, t_idx}];
full_vacc = params.template get<osecirvvs::DailyFullVaccination<ScalarType>>()[{(AgeGroup)i, t_idx}];
first_vacc =
params.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>()[{(AgeGroup)i, t_idx}];
full_vacc = params.template get<osecirvvs::DailyFullVaccinations<ScalarType>>()[{(AgeGroup)i, t_idx}];
}
else {
first_vacc = params.template get<osecirvvs::DailyFirstVaccination<ScalarType>>()[{(AgeGroup)i, t_idx}] -
params.template get<osecirvvs::DailyFirstVaccination<ScalarType>>()[{
(AgeGroup)i, t_idx - SimulationDay(1)}];
full_vacc = params.template get<osecirvvs::DailyFullVaccination<ScalarType>>()[{(AgeGroup)i, t_idx}] -
params.template get<osecirvvs::DailyFullVaccination<ScalarType>>()[{
first_vacc =
params.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>()[{(AgeGroup)i, t_idx}] -
params.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>()[{(AgeGroup)i,
t_idx - SimulationDay(1)}];
full_vacc = params.template get<osecirvvs::DailyFullVaccinations<ScalarType>>()[{(AgeGroup)i, t_idx}] -
params.template get<osecirvvs::DailyFullVaccinations<ScalarType>>()[{
(AgeGroup)i, t_idx - SimulationDay(1)}];
}

Expand Down Expand Up @@ -665,10 +667,11 @@ void setup_model(Model& model)

model.parameters.template get<osecirvvs::ICUCapacity<ScalarType>>() = 100;
model.parameters.template get<osecirvvs::TestAndTraceCapacity<ScalarType>>() = 0.0143;
model.parameters.template get<osecirvvs::DailyFirstVaccination<ScalarType>>().resize(SimulationDay(size_t(1000)));
model.parameters.template get<osecirvvs::DailyFirstVaccination<ScalarType>>().array().setConstant(5);
model.parameters.template get<osecirvvs::DailyFullVaccination<ScalarType>>().resize(SimulationDay(size_t(1000)));
model.parameters.template get<osecirvvs::DailyFullVaccination<ScalarType>>().array().setConstant(3);
model.parameters.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>().resize(
SimulationDay(size_t(1000)));
model.parameters.template get<osecirvvs::DailyPartialVaccinations<ScalarType>>().array().setConstant(5);
model.parameters.template get<osecirvvs::DailyFullVaccinations<ScalarType>>().resize(SimulationDay(size_t(1000)));
model.parameters.template get<osecirvvs::DailyFullVaccinations<ScalarType>>().array().setConstant(3);

auto& contacts = model.parameters.template get<osecirvvs::ContactPatterns<ScalarType>>();
auto& contact_matrix = contacts.get_cont_freq_mat();
Expand Down
8 changes: 4 additions & 4 deletions cpp/benchmarks/graph_simulation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -69,10 +69,10 @@ mio::osecirvvs::Model<ScalarType> create_model(size_t num_agegroups, const Scala
const size_t vacc_full = 5;
model.parameters.get<mio::osecirvvs::ICUCapacity<ScalarType>>() = 100;
model.parameters.get<mio::osecirvvs::TestAndTraceCapacity<ScalarType>>() = 0.0143;
model.parameters.get<mio::osecirvvs::DailyFirstVaccination<ScalarType>>().resize(mio::SimulationDay(tmax));
model.parameters.get<mio::osecirvvs::DailyFirstVaccination<ScalarType>>().array().setConstant(vacc_first);
model.parameters.get<mio::osecirvvs::DailyFullVaccination<ScalarType>>().resize(mio::SimulationDay(tmax));
model.parameters.get<mio::osecirvvs::DailyFullVaccination<ScalarType>>().array().setConstant(vacc_full);
model.parameters.get<mio::osecirvvs::DailyPartialVaccinations<ScalarType>>().resize(mio::SimulationDay(tmax));
model.parameters.get<mio::osecirvvs::DailyPartialVaccinations<ScalarType>>().array().setConstant(vacc_first);
model.parameters.get<mio::osecirvvs::DailyFullVaccinations<ScalarType>>().resize(mio::SimulationDay(tmax));
model.parameters.get<mio::osecirvvs::DailyFullVaccinations<ScalarType>>().array().setConstant(vacc_full);

auto& contacts = model.parameters.get<mio::osecirvvs::ContactPatterns<ScalarType>>();
auto& contact_matrix = contacts.get_cont_freq_mat();
Expand Down
4 changes: 4 additions & 0 deletions cpp/examples/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,10 @@ 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})

add_executable(ode_secirts_example ode_secirts.cpp)
target_link_libraries(ode_secirts_example PRIVATE memilio ode_secirts)
target_compile_options(ode_secirts_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

add_executable(ode_secir_ageres_example ode_secir_ageres.cpp)
target_link_libraries(ode_secir_ageres_example PRIVATE memilio ode_secir)
target_compile_options(ode_secir_ageres_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
Expand Down
154 changes: 154 additions & 0 deletions cpp/examples/ode_secirts.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
/*
* Copyright (C) 2020-2024 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 "ode_secirts/analyze_result.h"
#include "ode_secirts/model.h"
#include "ode_secirts/parameters.h"
#include "memilio/compartments/simulation.h"
#include "memilio/utils/logging.h"

int main()
{
// This example demonstrates how to simulate a SECIRTS model.
// The SECIRTS model is an extension of the SECIRVVS model that includes waning and temporary immunity.
// After the simulation, the aggregated size of the temporary immunity states are printed.
mio::set_log_level(mio::LogLevel::debug);

double t0 = 0;
double tmax = 100;
double dt = 0.1;

mio::log_info("Simulating SECIRTS; t={} ... {} with dt = {}.", t0, tmax, dt);

mio::osecirts::Model<double> model(3);
auto nb_groups = model.parameters.get_num_groups();

for (mio::AgeGroup i = 0; i < nb_groups; i++) {
// population
model.populations[{i, mio::osecirts::InfectionState::ExposedNaive}] = 20;
model.populations[{i, mio::osecirts::InfectionState::ExposedImprovedImmunity}] = 20;
model.populations[{i, mio::osecirts::InfectionState::ExposedPartialImmunity}] = 20;
model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsNaive}] = 30;
model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsNaiveConfirmed}] = 0;
model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsPartialImmunity}] = 30;
model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsPartialImmunityConfirmed}] = 0;
model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsImprovedImmunity}] = 30;
model.populations[{i, mio::osecirts::InfectionState::InfectedNoSymptomsImprovedImmunityConfirmed}] = 0;
model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsNaive}] = 40;
model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsNaiveConfirmed}] = 0;
model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsPartialImmunity}] = 40;
model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsPartialImmunityConfirmed}] = 0;
model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsImprovedImmunity}] = 40;
model.populations[{i, mio::osecirts::InfectionState::InfectedSymptomsImprovedImmunityConfirmed}] = 0;
model.populations[{i, mio::osecirts::InfectionState::InfectedSevereNaive}] = 30;
model.populations[{i, mio::osecirts::InfectionState::InfectedSevereImprovedImmunity}] = 30;
model.populations[{i, mio::osecirts::InfectionState::InfectedSeverePartialImmunity}] = 30;
model.populations[{i, mio::osecirts::InfectionState::InfectedCriticalNaive}] = 20;
model.populations[{i, mio::osecirts::InfectionState::InfectedCriticalPartialImmunity}] = 20;
model.populations[{i, mio::osecirts::InfectionState::InfectedCriticalImprovedImmunity}] = 20;
model.populations[{i, mio::osecirts::InfectionState::SusceptibleNaive}] = 1000;
model.populations[{i, mio::osecirts::InfectionState::SusceptiblePartialImmunity}] = 1200;
model.populations[{i, mio::osecirts::InfectionState::SusceptibleImprovedImmunity}] = 1000;
model.populations[{i, mio::osecirts::InfectionState::TemporaryImmunePartialImmunity}] = 60;
model.populations[{i, mio::osecirts::InfectionState::TemporaryImmuneImprovedImmunity}] = 70;
model.populations[{i, mio::osecirts::InfectionState::DeadNaive}] = 0;
model.populations[{i, mio::osecirts::InfectionState::DeadPartialImmunity}] = 0;
model.populations[{i, mio::osecirts::InfectionState::DeadImprovedImmunity}] = 0;

// parameters
//times
model.parameters.get<mio::osecirts::TimeExposed<double>>()[i] = 3.33;
model.parameters.get<mio::osecirts::TimeInfectedNoSymptoms<double>>()[i] = 1.87;
model.parameters.get<mio::osecirts::TimeInfectedSymptoms<double>>()[i] = 7;
model.parameters.get<mio::osecirts::TimeInfectedSevere<double>>()[i] = 6;
model.parameters.get<mio::osecirts::TimeInfectedCritical<double>>()[i] = 7;
model.parameters.get<mio::osecirts::TimeTemporaryImmunityPI<double>>()[i] = 60;
model.parameters.get<mio::osecirts::TimeTemporaryImmunityII<double>>()[i] = 60;
model.parameters.get<mio::osecirts::TimeWaningPartialImmunity<double>>()[i] = 180;
model.parameters.get<mio::osecirts::TimeWaningImprovedImmunity<double>>()[i] = 180;

//probabilities
model.parameters.get<mio::osecirts::TransmissionProbabilityOnContact<double>>()[i] = 0.15;
model.parameters.get<mio::osecirts::RelativeTransmissionNoSymptoms<double>>()[i] = 0.5;
model.parameters.get<mio::osecirts::RiskOfInfectionFromSymptomatic<double>>()[i] = 0.0;
model.parameters.get<mio::osecirts::MaxRiskOfInfectionFromSymptomatic<double>>()[i] = 0.4;
model.parameters.get<mio::osecirts::RecoveredPerInfectedNoSymptoms<double>>()[i] = 0.2;
model.parameters.get<mio::osecirts::SeverePerInfectedSymptoms<double>>()[i] = 0.1;
model.parameters.get<mio::osecirts::CriticalPerSevere<double>>()[i] = 0.1;
model.parameters.get<mio::osecirts::DeathsPerCritical<double>>()[i] = 0.1;

model.parameters.get<mio::osecirts::ReducExposedPartialImmunity<double>>()[i] = 0.8;
model.parameters.get<mio::osecirts::ReducExposedImprovedImmunity<double>>()[i] = 0.331;
model.parameters.get<mio::osecirts::ReducInfectedSymptomsPartialImmunity<double>>()[i] = 0.65;
model.parameters.get<mio::osecirts::ReducInfectedSymptomsImprovedImmunity<double>>()[i] = 0.243;
model.parameters.get<mio::osecirts::ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i] = 0.1;
model.parameters.get<mio::osecirts::ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i] = 0.091;
model.parameters.get<mio::osecirts::ReducTimeInfectedMild<double>>()[i] = 0.9;
}

model.parameters.get<mio::osecirts::ICUCapacity<double>>() = 100;
model.parameters.get<mio::osecirts::TestAndTraceCapacity<double>>() = 0.0143;
const size_t daily_vaccinations = 10;
const size_t num_days = 300;
model.parameters.get<mio::osecirts::DailyPartialVaccinations<double>>().resize(mio::SimulationDay(num_days));
model.parameters.get<mio::osecirts::DailyFullVaccinations<double>>().resize(mio::SimulationDay(num_days));
model.parameters.get<mio::osecirts::DailyBoosterVaccinations<double>>().resize(mio::SimulationDay(num_days));
for (size_t i = 0; i < num_days; ++i) {
for (mio::AgeGroup j = 0; j < nb_groups; ++j) {
auto num_vaccinations = static_cast<double>(i * daily_vaccinations);
model.parameters.get<mio::osecirts::DailyPartialVaccinations<double>>()[{j, mio::SimulationDay(i)}] =
num_vaccinations;
model.parameters.get<mio::osecirts::DailyFullVaccinations<double>>()[{j, mio::SimulationDay(i)}] =
num_vaccinations;
model.parameters.get<mio::osecirts::DailyBoosterVaccinations<double>>()[{j, mio::SimulationDay(i)}] =
num_vaccinations;
}
}

mio::ContactMatrixGroup& contact_matrix = model.parameters.get<mio::osecirts::ContactPatterns<double>>();
const double cont_freq = 10;
const double fact = 1.0 / (double)(size_t)nb_groups;
contact_matrix[0] =
mio::ContactMatrix(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, fact * cont_freq));
contact_matrix.add_damping(Eigen::MatrixXd::Constant((size_t)nb_groups, (size_t)nb_groups, 0.7),
mio::SimulationTime(30.));

model.parameters.get<mio::osecirts::Seasonality<double>>() = 0.2;

model.apply_constraints();

mio::TimeSeries<double> result = simulate(t0, tmax, dt, model);

bool print_to_terminal = true;

if (print_to_terminal) {
auto result_interpolated = mio::interpolate_simulation_result(result);
for (auto t_indx = 0; t_indx < result_interpolated.get_num_time_points(); t_indx++) {
double timm_pi = 0.0;
double timm_ii = 0.0;
for (mio::AgeGroup i = 0; i < nb_groups; i++) {
timm_pi += result_interpolated.get_value(t_indx)[model.populations.get_flat_index(
{i, mio::osecirts::InfectionState::TemporaryImmunePartialImmunity})];
timm_ii += result_interpolated.get_value(t_indx)[model.populations.get_flat_index(
{i, mio::osecirts::InfectionState::TemporaryImmuneImprovedImmunity})];
}
printf("t=%i, timm_pi=%f, timm_ii=%f\n", int(result_interpolated.get_time(t_indx)), timm_pi, timm_ii);
}
}
}
9 changes: 5 additions & 4 deletions cpp/examples/ode_secirvvs.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -67,15 +67,16 @@ int main()
model.parameters.get<mio::osecirvvs::ICUCapacity<double>>() = 100;
model.parameters.get<mio::osecirvvs::TestAndTraceCapacity<double>>() = 0.0143;
const size_t daily_vaccinations = 10;
model.parameters.get<mio::osecirvvs::DailyFirstVaccination<double>>().resize(mio::SimulationDay((size_t)tmax + 1));
model.parameters.get<mio::osecirvvs::DailyFullVaccination<double>>().resize(mio::SimulationDay((size_t)tmax + 1));
model.parameters.get<mio::osecirvvs::DailyPartialVaccinations<double>>().resize(
mio::SimulationDay((size_t)tmax + 1));
model.parameters.get<mio::osecirvvs::DailyFullVaccinations<double>>().resize(mio::SimulationDay((size_t)tmax + 1));
for (size_t i = 0; i < tmax + 1; ++i) {
auto num_vaccinations = static_cast<double>(i * daily_vaccinations);
model.parameters
.get<mio::osecirvvs::DailyFirstVaccination<double>>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
.get<mio::osecirvvs::DailyPartialVaccinations<double>>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
num_vaccinations;
model.parameters
.get<mio::osecirvvs::DailyFullVaccination<double>>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
.get<mio::osecirvvs::DailyFullVaccinations<double>>()[{(mio::AgeGroup)0, mio::SimulationDay(i)}] =
num_vaccinations;
}
model.parameters.get<mio::osecirvvs::DynamicNPIsImplementationDelay<double>>() = 7;
Expand Down
Loading

0 comments on commit f9b7610

Please sign in to comment.