Skip to content

Commit f864b19

Browse files
nijawareneSchmmknaranja
authored
1032 implementation of stochastic two variant seir model (#1047)
- Implemented stochastic SEIR model with two virus variants. - Add a model specific Simulation using the Euler-Maruyama method. Co-authored-by: reneSchm <49305466+reneSchm@users.noreply.github.com> Co-authored-by: Martin J. Kühn <62713180+mknaranja@users.noreply.github.com>
1 parent cb5affd commit f864b19

File tree

12 files changed

+1082
-0
lines changed

12 files changed

+1082
-0
lines changed

‎cpp/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -153,6 +153,7 @@ if(MEMILIO_BUILD_MODELS)
153153
add_subdirectory(models/ode_sir)
154154
add_subdirectory(models/sde_sir)
155155
add_subdirectory(models/sde_sirs)
156+
add_subdirectory(models/sde_seirvv)
156157
endif()
157158

158159
if(MEMILIO_BUILD_EXAMPLES)

‎cpp/examples/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,10 @@ add_executable(sde_sirs_example sde_sirs.cpp)
3434
target_link_libraries(sde_sirs_example PRIVATE memilio sde_sirs)
3535
target_compile_options(sde_sirs_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
3636

37+
add_executable(sde_seirvv_example sde_seirvv.cpp)
38+
target_link_libraries(sde_seirvv_example PRIVATE memilio sde_seirvv)
39+
target_compile_options(sde_seirvv_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
40+
3741
add_executable(ode_seair_example ode_seair.cpp)
3842
target_link_libraries(ode_seair_example PRIVATE memilio ode_seair AD::AD)
3943
target_compile_options(ode_seair_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/examples/sde_seirvv.cpp‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
/*
2+
* Copyright (C) 2020-2024 MEmilio
3+
*
4+
* Authors: Nils Wassmuth, Rene Schmieding, Martin J. Kuehn
5+
*
6+
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
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/utils/logging.h"
21+
#include "memilio/utils/uncertain_value.h"
22+
#include "sde_seirvv/model.h"
23+
#include "sde_seirvv/simulation.h"
24+
25+
#include <vector>
26+
27+
int main()
28+
{
29+
mio::set_log_level(mio::LogLevel::debug);
30+
31+
ScalarType t0 = 0.;
32+
ScalarType tmid = 100.;
33+
ScalarType tmax = 400.;
34+
ScalarType dt = 0.1;
35+
36+
mio::log_info("Simulating SEIRVV; t={} ... {} with dt = {}.", t0, tmax, dt);
37+
38+
mio::sseirvv::Model model;
39+
40+
ScalarType total_population = 180000;
41+
42+
model.populations[{mio::sseirvv::InfectionState::ExposedV1}] = 0;
43+
model.populations[{mio::sseirvv::InfectionState::ExposedV2}] = 0;
44+
model.populations[{mio::sseirvv::InfectionState::InfectedV1}] = 7200;
45+
model.populations[{mio::sseirvv::InfectionState::InfectedV2}] = 0;
46+
model.populations[{mio::sseirvv::InfectionState::RecoveredV1}] = 0;
47+
model.populations[{mio::sseirvv::InfectionState::RecoveredV2}] = 0;
48+
model.populations[{mio::sseirvv::InfectionState::ExposedV1V2}] = 0;
49+
model.populations[{mio::sseirvv::InfectionState::InfectedV1V2}] = 0;
50+
model.populations[{mio::sseirvv::InfectionState::RecoveredV1V2}] = 0;
51+
model.populations[{mio::sseirvv::InfectionState::Susceptible}] =
52+
total_population -
53+
model.populations[{mio::sseirvv::InfectionState::ExposedV1}] -
54+
model.populations[{mio::sseirvv::InfectionState::ExposedV2}] -
55+
model.populations[{mio::sseirvv::InfectionState::InfectedV1}] -
56+
model.populations[{mio::sseirvv::InfectionState::InfectedV2}] -
57+
model.populations[{mio::sseirvv::InfectionState::RecoveredV1}] -
58+
model.populations[{mio::sseirvv::InfectionState::RecoveredV2}] -
59+
model.populations[{mio::sseirvv::InfectionState::ExposedV1V2}] -
60+
model.populations[{mio::sseirvv::InfectionState::InfectedV1V2}] -
61+
model.populations[{mio::sseirvv::InfectionState::RecoveredV1V2}];
62+
63+
// It is assumed that both variants have the same transmission probability
64+
// on contact and the same time exposed. The time infected is scaled by
65+
// 1.35 for the second variant.
66+
model.parameters.get<mio::sseirvv::ContactPatterns>().get_baseline()(0, 0) = 1;
67+
model.parameters.set<mio::sseirvv::TransmissionProbabilityOnContactV1>(0.076);
68+
model.parameters.set<mio::sseirvv::TransmissionProbabilityOnContactV2>(0.076);
69+
model.parameters.set<mio::sseirvv::TimeExposedV1>(5.33);
70+
model.parameters.set<mio::sseirvv::TimeExposedV2>(5.33);
71+
model.parameters.set<mio::sseirvv::TimeInfectedV1>(17.2);
72+
model.parameters.set<mio::sseirvv::TimeInfectedV2>(17.2 * 1.35);
73+
74+
model.check_constraints();
75+
76+
// Simulate the model up until tmid, with only the first variant.
77+
auto sseirv = mio::sseirvv::simulate(t0, tmid, dt, model);
78+
// Set the model population to the simulation result, so it is used as initial value for the second simulation.
79+
model.populations.array() = sseirv.get_last_value().cast<mio::UncertainValue<ScalarType>>();
80+
// The second variant enters with 100 individuals. This increases the model population to total_population + 100.
81+
model.populations[{mio::sseirvv::InfectionState::InfectedV2}] = 100;
82+
// Simulate the model from tmid to tmax, now with both variants.
83+
auto sseirv2 = mio::sseirvv::simulate(tmid, tmax, dt, model);
84+
85+
sseirv.print_table({"Susceptible", "ExposedV1", "InfectedV1", "RecoveredV1", "ExposedV2", "InfectedV2", "RecoveredV2", "ExposedV1V2", "InfectedV1V2", "RecoveredV1V2"});
86+
sseirv2.print_table({"Susceptible", "ExposedV1", "InfectedV1", "RecoveredV1", "ExposedV2", "InfectedV2", "RecoveredV2", "ExposedV1V2", "InfectedV1V2", "RecoveredV1V2"});
87+
}

‎cpp/models/sde_seirvv/CMakeLists.txt‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
add_library(sde_seirvv
2+
infection_state.h
3+
model.h
4+
model.cpp
5+
parameters.h
6+
simulation.h
7+
)
8+
target_link_libraries(sde_seirvv PUBLIC memilio)
9+
target_include_directories(sde_seirvv PUBLIC
10+
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
11+
$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
12+
)
13+
target_compile_options(sde_seirvv PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})

‎cpp/models/sde_seirvv/README.md‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
# SDE-based SEIR-type model with two variants
2+
3+
This module models and simulates the epidemic using an SDE-based SEIR-type model approach with two variants. After a recovered infection from the first variant (considered the wild type) you can still get infected with the second variant. Infection with the second variant immunizes you against the first variant. Unlike the agent based model that uses particular agents, this model simulates the spread of a communicable disease in a population with subpopulations being in different compartments. The model also introduces stochasticity compared to an ODE-based model. The compartments are as follows: `Susceptible`, `ExposedV1`, `ExposedV2`, `ExposedV1V2`, `InfectedV1`, `InfectedV2`, `InfectedV1V2`, `RecoveredV1`, `RecoveredV2` and `RecoveredV1V2`. The compartment names are appended by the relevant variants. The addendum `V1` means that the individual is (or has been) infected with the first variant, the addendum `V2` means that the individual is (or has been) infected with the second variant with no prior infections and the addendum `V1V2` means that the individual is (or has been) infected with the second variant after a successful recovery from the first variant. Only individuals in the infected compartments are infectious.
4+
5+
Below is an overview of the model architecture.
6+
7+
8+
9+
## Simulation
10+
11+
The simulation runs in discrete time steps using an Euler-Maruyama integration scheme. The Simulation class handles the parameters and the numerical integrator. It also stores the result.
12+
![SEIRVV_model](https://github.com/user-attachments/assets/55258e5d-05f5-4b16-93b0-f089f8f70782)
13+
| Mathematical variable | C++ variable name | Description |
14+
|---------------------------- | --------------- | -------------------------------------------------------------------------------------------------- |
15+
| $\phi$ | `ContactPatterns` | Daily contact rate / Number of daily contacts. |
16+
| $\rho_1$ | `TransmissionProbabilityOnContactV1` | Transmission risk for people located in the Susceptible compartment (susceptible to infection with variant 1). |
17+
| $\rho_2$ | `TransmissionProbabilityOnContactV2` | Transmission risk for people located in the Susceptible compartment or in the RecoveredV1 compartment (susceptivle to infection with variant 2). |
18+
| $N$ | `populations.get_total()` | Total population. |
19+
| $T_{E,1}$ | `TimeExposed1` | Average time in days an individual stays in the ExposedV1 compartment. |
20+
| $T_{E,2}$ | `TimeExposed2` | Average time in days an individual stays in the ExposedV2 or in the ExposedV1V2 compartment. |
21+
| $T_{I,1}$ | `TimeInfected1` | Average time in days an individual stays in the InfectedV1 compartment. |
22+
| $T_{I,2}$ | `TimeInfected2` | Average time in days an individual stays in the InfectedV2 or in the InfectedV1V2 compartment. |
23+
24+
An example can be found in [examples/sde_seirvv.cpp](../../examples/sde_seirvv.cpp)

‎cpp/models/sde_seirvv/infection_state.h‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
/*
2+
* Copyright (C) 2020-2024 MEmilio
3+
*
4+
* Authors: Nils Wassmuth, Rene Schmieding, Martin J. Kuehn
5+
*
6+
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
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+
#ifndef MIO_SDE_SEIRVV_INFECTIONSTATE_H
22+
#define MIO_SDE_SEIRVV_INFECTIONSTATE_H
23+
24+
namespace mio
25+
{
26+
namespace sseirvv
27+
{
28+
29+
/**
30+
* @brief The InfectionState enum describes the possible
31+
* categories for the infectious state of persons
32+
*/
33+
enum class InfectionState
34+
{
35+
Susceptible,
36+
ExposedV1,
37+
InfectedV1,
38+
RecoveredV1,
39+
ExposedV2,
40+
InfectedV2,
41+
RecoveredV2,
42+
ExposedV1V2,
43+
InfectedV1V2,
44+
RecoveredV1V2,
45+
Count
46+
};
47+
48+
} // namespace sseirvv
49+
} // namespace mio
50+
51+
#endif // MIO_SDE_SEIRVV_INFECTIONSTATE_H

‎cpp/models/sde_seirvv/model.cpp‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
/*
2+
* Copyright (C) 2020-2024 MEmilio
3+
*
4+
* Authors: Nils Wassmuth, Rene Schmieding, Martin J. Kuehn
5+
*
6+
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
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 "sde_seirvv/model.h"
22+
23+
namespace mio
24+
{
25+
namespace sseirvv
26+
{
27+
28+
} // namespace sseirvv
29+
} // namespace mio

0 commit comments

Comments
 (0)