|
| 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 | + |
| 21 | +#include "ode_seirv/model.h" |
| 22 | +#include "memilio/compartments/simulation.h" |
| 23 | +#include "memilio/utils/logging.h" |
| 24 | + |
| 25 | +/** |
| 26 | + * @brief set_initial_population sets the initial population of the model |
| 27 | + * |
| 28 | + * We assume that no one is initially infected or exposed at the beginning of the season. |
| 29 | + * Infections are seeded later via the OutsideFoI parameter. |
| 30 | + * The distinction into the 2 layers of susceptibility is done via the parameters. |
| 31 | + * @tparam FP floating point type |
| 32 | + * @param[in, out] model SEIRV model |
| 33 | + * @param[in] total_pop total population size |
| 34 | + */ |
| 35 | +template <class FP> |
| 36 | +void set_initial_population(mio::oseirv::Model<FP>& model, const FP total_pop) |
| 37 | +{ |
| 38 | + auto& params = model.parameters; |
| 39 | + auto& pop = model.populations; |
| 40 | + const size_t num_groups = (size_t)params.get_num_groups(); |
| 41 | + |
| 42 | + for (size_t i = 0; i < num_groups; ++i) { |
| 43 | + pop.template set_difference_from_group_total<mio::AgeGroup>( |
| 44 | + {mio::AgeGroup(i), mio::oseirv::InfectionState::Susceptible}, total_pop / num_groups); |
| 45 | + |
| 46 | + // Total population N_i as currently stored in group totals |
| 47 | + FP Ni = pop.get_group_total(mio::AgeGroup(i)); |
| 48 | + |
| 49 | + FP s_age = params.template get<mio::oseirv::SusceptibilityByAge<FP>>()[mio::AgeGroup(i)]; |
| 50 | + FP s_frac = params.template get<mio::oseirv::SusceptibleFraction<FP>>(); |
| 51 | + FP vc = params.template get<mio::oseirv::VaccineCoverage<FP>>()[mio::AgeGroup(i)]; |
| 52 | + FP ve = params.template get<mio::oseirv::VaccineEffectiveness<FP>>()[mio::AgeGroup(i)]; |
| 53 | + |
| 54 | + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Exposed}] = 0; |
| 55 | + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::ExposedVaccinated}] = 0; |
| 56 | + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Infected}] = 0; |
| 57 | + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::InfectedVaccinated}] = 0; |
| 58 | + |
| 59 | + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Susceptible}] = s_frac * s_age * (FP(1) - vc) * Ni; |
| 60 | + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::SusceptibleVaccinated}] = |
| 61 | + s_frac * s_age * (FP(1) - ve) * vc * Ni; |
| 62 | + |
| 63 | + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::Recovered}] = (FP(1) - s_frac * s_age) * (FP(1) - vc) * Ni; |
| 64 | + pop[{mio::AgeGroup(i), mio::oseirv::InfectionState::RecoveredVaccinated}] = |
| 65 | + (FP(1) - s_frac * s_age * (FP(1) - ve)) * vc * Ni; |
| 66 | + } |
| 67 | +} |
| 68 | + |
| 69 | +// Example usage of the SEIRV ODE model presented in https://doi.org/10.1186/s12879-017-2344-6 |
| 70 | +int main() |
| 71 | +{ |
| 72 | + using FP = double; |
| 73 | + mio::set_log_level(mio::LogLevel::debug); |
| 74 | + |
| 75 | + const FP t0 = 0., tmax = 42.; |
| 76 | + const FP dt = 0.1; |
| 77 | + |
| 78 | + const size_t num_groups = 1; |
| 79 | + const auto total_pop = 1e5; |
| 80 | + mio::oseirv::Model<FP> model((int)num_groups); |
| 81 | + auto& parameters = model.parameters; |
| 82 | + |
| 83 | + FP cont_freq = 10.0; |
| 84 | + FP cont_freq_group = FP(1) / FP(num_groups); |
| 85 | + |
| 86 | + // Healthy contacts |
| 87 | + mio::ContactMatrixGroup<ScalarType>& contacts_healthy = parameters.get<mio::oseirv::ContactPatternsHealthy<FP>>(); |
| 88 | + contacts_healthy[0] = mio::ContactMatrix<FP>( |
| 89 | + Eigen::MatrixXd::Constant((Eigen::Index)num_groups, (Eigen::Index)num_groups, cont_freq_group * cont_freq)); |
| 90 | + |
| 91 | + // Sick contacts (here 20% reduced) |
| 92 | + mio::ContactMatrixGroup<ScalarType>& contacts_sick = parameters.get<mio::oseirv::ContactPatternsSick<FP>>(); |
| 93 | + contacts_sick[0] = mio::ContactMatrix<FP>(Eigen::MatrixXd::Constant( |
| 94 | + (Eigen::Index)num_groups, (Eigen::Index)num_groups, cont_freq_group * cont_freq * 0.8)); |
| 95 | + |
| 96 | + // Parameters |
| 97 | + parameters.get<mio::oseirv::BaselineTransmissibility<FP>>() = 1.2; |
| 98 | + parameters.get<mio::oseirv::TimeExposed<FP>>() = 2.0; |
| 99 | + parameters.get<mio::oseirv::TimeInfected<FP>>() = 2.0; // same duration for E->I and I->R |
| 100 | + parameters.get<mio::oseirv::SeasonalityAmplitude<FP>>() = 1.0; |
| 101 | + parameters.get<mio::oseirv::SeasonalityShiftPerSubtype<FP>>() = 0.0; |
| 102 | + parameters.get<mio::oseirv::SeasonalityShiftPerSeason<FP>>() = 0.0; |
| 103 | + parameters.get<mio::oseirv::OutsideFoI<FP>>() = 1e-6; |
| 104 | + parameters.get<mio::oseirv::ClusteringExponent<FP>>() = 0.9; |
| 105 | + parameters.get<mio::oseirv::SickMixing<FP>>() = 2.0; |
| 106 | + |
| 107 | + for (size_t i = 0; i < num_groups; ++i) { |
| 108 | + parameters.get<mio::oseirv::SusceptibilityByAge<FP>>()[mio::AgeGroup(i)] = 1.0; |
| 109 | + parameters.get<mio::oseirv::VaccineCoverage<FP>>()[mio::AgeGroup(i)] = 0.3; |
| 110 | + parameters.get<mio::oseirv::VaccineEffectiveness<FP>>()[mio::AgeGroup(i)] = 0.5; |
| 111 | + model.populations.set_difference_from_group_total<mio::AgeGroup>( |
| 112 | + {mio::AgeGroup(i), mio::oseirv::InfectionState::Susceptible}, 1e5 / num_groups); |
| 113 | + } |
| 114 | + |
| 115 | + set_initial_population(model, total_pop); |
| 116 | + |
| 117 | + auto seirv = mio::simulate<FP>(t0, tmax, dt, model); |
| 118 | + |
| 119 | + std::vector<std::string> vars = {"S", "SV", "E", "EV", "I", "IV", "R", "RV"}; |
| 120 | + seirv.print_table(vars); |
| 121 | +} |
0 commit comments