Skip to content

1279 minor additions to how infections work #1320

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 10 commits into
base: main
Choose a base branch
from
Draft

1279 minor additions to how infections work #1320

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
ee81819
refactor
xsaschako Jul 17, 2025
2fb9b00
Improve infection draw, rename infectivity and related functions
DavidKerkmann Jul 21, 2025
24f4e4d
Quick naming fixes
DavidKerkmann Jul 21, 2025
98b415b
merge main
DavidKerkmann Jul 21, 2025
77ae766
Fix vaccination time dependency
DavidKerkmann Jul 21, 2025
763be59
Changed multiply function in time.h
DavidKerkmann Jul 22, 2025
f3abc80
Small update to function description
DavidKerkmann Jul 22, 2025
8f66c0b
Merge branch 'main' into 1279-minor-additions-to-how-infections-work
DavidKerkmann Nov 3, 2025
42cbb92
merge main
DavidKerkmann Nov 3, 2025
6bc0825
Improvements according to review --- IGNORE ---
DavidKerkmann Nov 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 4 additions & 8 deletions cpp/models/abm/analyze_result.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,10 +169,6 @@ std::vector<Model> ensemble_params_percentile(const std::vector<std::vector<Mode
return model.parameters.template get<DeathsPerInfectedCritical>()[{virus_variant, age_group}];
});

param_percentile(node, [age_group, virus_variant](auto&& model) -> auto& {
return model.parameters.template get<DetectInfection>()[{virus_variant, age_group}];
});

param_percentile(node, [virus_variant](auto&& model) -> auto& {
return model.parameters.template get<AerosolTransmissionRates>()[{virus_variant}];
});
Expand All @@ -187,17 +183,17 @@ std::vector<Model> ensemble_params_percentile(const std::vector<std::vector<Mode
return dist1.viral_load_peak < dist2.viral_load_peak;
});
param_percentile_dist(
node, std::vector<InfectivityDistributionsParameters>(num_runs),
node, std::vector<ViralShedTuple>(num_runs),
[age_group, virus_variant](auto&& model) -> auto& {
return model.parameters.template get<InfectivityDistributions>()[{virus_variant, age_group}];
return model.parameters.template get<ViralShedParameters>()[{virus_variant, age_group}];
},
[](auto& dist1, auto& dist2) {
return dist1.infectivity_alpha < dist2.infectivity_alpha;
return dist1.virus_shed_alpha < dist2.virus_shed_alpha;
});
param_percentile_dist(
node, std::vector<mio::AbstractParameterDistribution>(num_runs),
[age_group, virus_variant](auto&& model) -> auto& {
return model.parameters.template get<VirusShedFactor>()[{virus_variant, age_group}];
return model.parameters.template get<ViralShedFactor>()[{virus_variant, age_group}];
},
[](auto& dist1, auto& dist2) {
return dist1 < dist2;
Expand Down
524 changes: 295 additions & 229 deletions cpp/models/abm/infection.cpp
View file
Open in desktop

Large diffs are not rendered by default.

123 changes: 112 additions & 11 deletions cpp/models/abm/infection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,15 @@ namespace mio
namespace abm
{

/**
* @brief Represents a transition between infection states with duration and probability.
*/
struct StateTransition {
InfectionState from_state;
InfectionState to_state;
TimeSpan duration;
};

/**
* @brief Models the ViralLoad for an Infection, modelled on a log_10 scale.
* Based on https://www.science.org/doi/full/10.1126/science.abi5273
Expand Down Expand Up @@ -74,7 +83,7 @@ class Infection
* @param[in] detected [Default: false] If the Infection is detected.
*/
Infection(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age, const Parameters& params,
TimePoint start_date, InfectionState start_state = InfectionState::Exposed,
TimePoint start_date, InfectionState start_state = InfectionState::Susceptible,
ProtectionEvent latest_protection = {ProtectionType::NoProtection, TimePoint(0)}, bool detected = false);

/**
Expand All @@ -84,17 +93,20 @@ class Infection
ScalarType get_viral_load(TimePoint t) const;

/**
* @brief Get infectivity at a given time.
* @brief Get viral shed at a specific time.
* Computed depending on current ViralLoad and individual invlogit function of each Person
* corresponding to https://www.science.org/doi/full/10.1126/science.abi5273
* The mapping corresponds to Fig. 2 C.
* Formula of invlogit function can be found here:
* https://github.com/VirologyCharite/SARS-CoV-2-VL-paper/tree/main
* in ExtendedMethods.html, Section 3.1.2.1.
* * Also in accordance to Fig. 3d of another publication:
* https://www.nature.com/articles/s41564-022-01105-z/figures/3
* The result is in arbitrary units and has to be scaled to the rate "infections per day".
* @param[in] t TimePoint of the querry.
* @return Infectivity at given TimePoint.
* @return Viral shed at given TimePoint.
*/
ScalarType get_infectivity(TimePoint t) const;
ScalarType get_viral_shed(TimePoint t) const;

/**
* @brief: Get VirusVariant.
Expand Down Expand Up @@ -163,19 +175,22 @@ class Infection
* InfectedSymptoms -> Infected_Severe or InfectedSymptoms -> Recovered,
* InfectedSevere -> InfectedCritical or InfectedSevere -> Recovered or InfectedSevere -> Dead,
* InfectedCritical -> Recovered or InfectedCritical -> Dead,
* with artifical, hardcoded probabilites, until either Recoverd or Dead is reached.
* This is subject to change when parameter distributions for these transitions are implemented.
* until either Recoverd or Dead is reached.
* The duration in each #InfectionState is taken from the respective parameter.
* The first transition has a random, uniformly distributed length from 0 to the respective parameter
* for persons that are initialized somewhere in the middle of the infection, i.e. not a newly infected
* person nor a Recovered/Dead person, to reflect uncertainty in the current state.
* @param[inout] rng PersonalRandomNumberGenerator of the Person.
* @param[in] age AgeGroup of the Person.
* @param[in] params Parameters of the Model.
* @param[in] init_date Date of initializing the Infection.
* @param[in] init_state #InfectionState at time of initializing the Infection.
* @param[in] latest_protection Latest protection against Infection, has an influence on transition probabilities.
* @return The starting date of the init_state.
*/
void draw_infection_course_forward(PersonalRandomNumberGenerator& rng, AgeGroup age, const Parameters& params,
TimePoint init_date, InfectionState init_state,
ProtectionEvent latest_protection);
TimePoint draw_infection_course_forward(PersonalRandomNumberGenerator& rng, AgeGroup age, const Parameters& params,
TimePoint init_date, InfectionState init_state,
ProtectionEvent latest_protection);

/**
* @brief Determine Infection course prior to the given #InfectionState start_state.
Expand All @@ -191,12 +206,98 @@ class Infection
TimePoint draw_infection_course_backward(PersonalRandomNumberGenerator& rng, AgeGroup age, const Parameters& params,
TimePoint init_date, InfectionState init_state);

/**
* @brief Initialize the viral load parameters for the infection.
* @param[inout] rng PersonalRandomNumberGenerator of the Person.
* @param[in] virus Virus type of the Infection.
* @param[in] age AgeGroup of the Person.
* @param[in] params Parameters of the Model.
* @param[in] init_date Date of initializing the Infection.
* @param[in] latest_protection Latest protection against Infection.
*/
void initialize_viral_load(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age,
const Parameters& params, TimePoint init_date, ProtectionEvent latest_protection);

/**
* @brief Initialize the viral shed parameters and individual factor for the infection.
* @param[inout] rng PersonalRandomNumberGenerator of the Person.
* @param[in] virus Virus type of the Infection.
* @param[in] age AgeGroup of the Person.
* @param[in] params Parameters of the Model.
*/
void initialize_viral_shed(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age,
const Parameters& params);

/**
* @brief Get the forward transition from a given infection state.
* @param[inout] rng PersonalRandomNumberGenerator of the Person.
* @param[in] age AgeGroup of the Person.
* @param[in] params Parameters of the Model.
* @param[in] current_state Current infection state.
* @param[in] current_time Current time point.
* @param[in] latest_protection Latest protection against Infection.
* @return StateTransition representing the next transition.
*/
StateTransition get_forward_transition(PersonalRandomNumberGenerator& rng, AgeGroup age, const Parameters& params,
InfectionState current_state, TimePoint current_time,
ProtectionEvent latest_protection) const;

/**
* @brief Get the backward transition from a given infection state.
* @param[inout] rng PersonalRandomNumberGenerator of the Person.
* @param[in] age AgeGroup of the Person.
* @param[in] params Parameters of the Model.
* @param[in] current_state Current infection state.
* @return StateTransition representing the previous transition.
*/
StateTransition get_backward_transition(PersonalRandomNumberGenerator& rng, AgeGroup age, const Parameters& params,
InfectionState current_state) const;

/**
* @brief Get the backward transition from recovered state.
* @param[inout] rng PersonalRandomNumberGenerator of the Person.
* @param[in] age AgeGroup of the Person.
* @param[in] params Parameters of the Model.
* @return StateTransition representing the transition that led to recovery.
*/
StateTransition get_recovered_backward_transition(PersonalRandomNumberGenerator& rng, AgeGroup age,
const Parameters& params) const;

/**
* @brief Get the backward transition from dead state.
* @param[inout] rng PersonalRandomNumberGenerator of the Person.
* @param[in] age AgeGroup of the Person.
* @param[in] params Parameters of the Model.
* @return StateTransition representing the transition that led to death.
*/
StateTransition get_dead_backward_transition(PersonalRandomNumberGenerator& rng, AgeGroup age,
const Parameters& params) const;

/**
* @brief Calculate the overall death probability for the infection.
* @param[in] age AgeGroup of the Person.
* @param[in] params Parameters of the Model.
* @return The probability of death for this infection.
*/
ScalarType calculate_death_probability(AgeGroup age, const Parameters& params) const;

/**
* @brief Get the severity protection factor based on latest protection.
* @param[in] params Parameters of the Model.
* @param[in] latest_protection Latest protection against Infection.
* @param[in] age AgeGroup of the Person.
* @param[in] current_time Current time point.
* @return The protection factor against severe outcomes.
*/
ScalarType get_severity_protection_factor(const Parameters& params, ProtectionEvent latest_protection, AgeGroup age,
TimePoint current_time) const;

std::vector<std::pair<TimePoint, InfectionState>> m_infection_course; ///< Start date of each #InfectionState.
VirusVariant m_virus_variant; ///< Variant of the Infection.
ViralLoad m_viral_load; ///< ViralLoad of the Infection.
ScalarType m_log_norm_alpha,
m_log_norm_beta; ///< Parameters for the infectivity mapping, which is modelled through an invlogit function.
ScalarType m_individual_virus_shed_factor; ///< Individual virus shed factor.
m_log_norm_beta; ///< Parameters for the viral shed mapping, which is modelled through an invlogit function.
ScalarType m_individual_viral_shed_factor; ///< Individual viral shed factor.
bool m_detected; ///< Whether an Infection is detected or not.
};

Expand Down
54 changes: 28 additions & 26 deletions cpp/models/abm/model_functions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,30 +32,30 @@ namespace mio
namespace abm
{

ScalarType daily_transmissions_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
const VirusVariant virus, const AgeGroup age_receiver,
size_t age_receiver_group_size, const LocalInfectionParameters& params)
ScalarType total_exposure_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
const VirusVariant virus, const AgeGroup age_receiver,
size_t age_receiver_group_size, const LocalInfectionParameters& params)
{
assert(age_receiver < rates.size<AgeGroup>());
ScalarType prob = 0;
ScalarType total_exposure = 0;
for (AgeGroup age_transmitter(0); age_transmitter < rates.size<AgeGroup>(); ++age_transmitter) {
if (age_receiver == age_transmitter &&
age_receiver_group_size > 1) // adjust for the person not meeting themself
{
prob += rates[{cell_index, virus, age_transmitter}] *
params.get<ContactRates>()[{age_receiver, age_transmitter}] * age_receiver_group_size /
(age_receiver_group_size - 1);
total_exposure += rates[{cell_index, virus, age_transmitter}] *
params.get<ContactRates>()[{age_receiver, age_transmitter}] * age_receiver_group_size /
(age_receiver_group_size - 1);
}
else {
prob += rates[{cell_index, virus, age_transmitter}] *
params.get<ContactRates>()[{age_receiver, age_transmitter}];
total_exposure += rates[{cell_index, virus, age_transmitter}] *
params.get<ContactRates>()[{age_receiver, age_transmitter}];
}
}
return prob;
return total_exposure;
}

ScalarType daily_transmissions_by_air(const AirExposureRates& rates, const CellIndex cell_index,
const VirusVariant virus, const Parameters& global_params)
ScalarType total_exposure_by_air(const AirExposureRates& rates, const CellIndex cell_index, const VirusVariant virus,
const Parameters& global_params)
{
return rates[{cell_index, virus}] * global_params.get<AerosolTransmissionRates>()[{virus}];
}
Expand Down Expand Up @@ -83,25 +83,27 @@ void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const
ScalarType mask_protection = person.get_mask_protective_factor(global_parameters);
assert(person.get_cells().size() && "Person is in multiple cells. Interact logic is incorrect at the moment.");
for (CellIndex cell_index : person.get_cells()) {
std::pair<VirusVariant, ScalarType> local_indiv_trans_prob[static_cast<uint32_t>(VirusVariant::Count)];
std::pair<VirusVariant, ScalarType> local_indiv_expected_trans[static_cast<uint32_t>(VirusVariant::Count)];
for (uint32_t v = 0; v != static_cast<uint32_t>(VirusVariant::Count); ++v) {
VirusVariant virus = static_cast<VirusVariant>(v);
size_t local_population_by_age_receiver = local_population_by_age[{cell_index, age_receiver}];
ScalarType local_indiv_trans_prob_v =
(daily_transmissions_by_contacts(local_contact_exposure, cell_index, virus, age_receiver,
local_population_by_age_receiver, local_parameters) +
daily_transmissions_by_air(local_air_exposure, cell_index, virus, global_parameters)) *
ScalarType exposed_viral_shed =
(total_exposure_by_contacts(local_contact_exposure, cell_index, virus, age_receiver,
local_population_by_age_receiver, local_parameters) +
total_exposure_by_air(local_air_exposure, cell_index, virus, global_parameters)) *
(1 - mask_protection) * (1 - person.get_protection_factor(t, virus, global_parameters));

local_indiv_trans_prob[v] = std::make_pair(virus, local_indiv_trans_prob_v);
ScalarType infection_rate =
global_parameters.get<InfectionRateFromViralShed>()[{virus}] * exposed_viral_shed;
local_indiv_expected_trans[v] = std::make_pair(virus, infection_rate);
}
VirusVariant virus =
random_transition(personal_rng, VirusVariant::Count, dt,
local_indiv_trans_prob); // use VirusVariant::Count for no virus submission
local_indiv_expected_trans); // use VirusVariant::Count for no virus submission
if (virus != VirusVariant::Count) {
person.add_new_infection(Infection(personal_rng, virus, age_receiver, global_parameters, t + dt / 2,
mio::abm::InfectionState::Exposed, person.get_latest_protection(),
false)); // Starting time in first approximation
mio::abm::InfectionState::Susceptible,
person.get_latest_protection(t + dt / 2),
false)); // Starting time in second order approximation
}
}
}
Expand All @@ -121,14 +123,14 @@ void add_exposure_contribution(AirExposureRates& local_air_exposure, ContactExpo
auto& infection = person.get_infection();
auto virus = infection.get_virus_variant();
auto age = person.get_age();
// average infectivity over the time step to second order accuracy using midpoint rule
const auto infectivity = infection.get_infectivity(t + dt / 2);
// average viral shed over the time step to second order accuracy using midpoint rule
const auto viral_shed = infection.get_viral_shed(t + dt / 2);
const auto quarantine_factor =
person.is_in_quarantine(t, params) ? (1.0 - params.get<QuarantineEffectiveness>()) : 1.0;

for (CellIndex cell : person.get_cells()) {
auto air_contribution = infectivity * quarantine_factor;
auto contact_contribution = infectivity * quarantine_factor;
auto air_contribution = viral_shed * quarantine_factor;
auto contact_contribution = viral_shed * quarantine_factor;

if (location.get_infection_parameters().get<UseLocationCapacityForTransmissions>()) {
air_contribution *= location.get_cells()[cell.get()].compute_space_per_person_relative();
Expand Down
Loading