1279 minor additions to how infections work

cpp/models/abm/analyze_result.h

@@ -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}];
                 });
@@ -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.viral_shed_alpha < dist2.viral_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;

cpp/models/abm/infection.h

@@ -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
@@ -58,6 +67,13 @@ struct ViralLoad {
     }
 };
 
+/**
+ * @brief Distributions of the relative time that people have been in their initial infection state at the beginning of the simulation.
+ * Values have to be within [0, 1].
+ * This makes it possible to draw from a user-defined distribution instead of drawing from a uniform distribution.
+ */
+using InitialInfectionStateDistribution = CustomIndexArray<AbstractParameterDistribution, VirusVariant, AgeGroup>;
+
 class Infection
 {
 public:
@@ -77,24 +93,43 @@ class Infection
               TimePoint start_date, InfectionState start_state = InfectionState::Exposed,
               ProtectionEvent latest_protection = {ProtectionType::NoProtection, TimePoint(0)}, bool detected = false);
 
+    /**
+     * @brief Create an Infection for a single Person with a time spent in the given initial state that is drawn from the given distribution.
+     * @param[inout] rng PersonalRandomNumberGenerator of the Person.
+     * @param[in] virus Virus type of the Infection.
+     * @param[in] age AgeGroup to determine the ViralLoad course.
+     * @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] init_state_dist Distribution to draw the relative time spent in the initial state from. Values have to be within [0, 1].
+     * @param[in] latest_protection The pair value of last ProtectionType (previous Infection/Vaccination) and TimePoint of that protection.
+     * @param[in] detected If the Infection is detected.
+     */
+    Infection(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age, const Parameters& params,
+              TimePoint init_date, InfectionState init_state, const InitialInfectionStateDistribution& init_state_dist,
+              ProtectionEvent latest_protection = {ProtectionType::NoProtection, TimePoint(0)}, bool detected = false);
+
     /**
      * @brief Gets the ViralLoad of the Infection at a given TimePoint.
      * @param[in] t TimePoint of querry.
      */
     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.
@@ -133,7 +168,7 @@ class Infection
             .add("viral_load", m_viral_load)
             .add("log_norm_alpha", m_log_norm_alpha)
             .add("log_norm_beta", m_log_norm_beta)
-            .add("individual_virus_shed_factor", m_individual_virus_shed_factor)
+            .add("individual_viral_shed_factor", m_individual_viral_shed_factor)
             .add("detected", m_detected);
     }
 
@@ -164,8 +199,7 @@ 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.
      * @param[inout] rng PersonalRandomNumberGenerator of the Person.
      * @param[in] age AgeGroup of the Person.
@@ -192,12 +226,97 @@ 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] latest_protection Latest protection against Infection.
+     */
+    void initialize_viral_load(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age,
+                               const Parameters& params, 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.
 };
 

cpp/models/abm/model_functions.cpp

@@ -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}];
 }
@@ -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::Exposed,
+                                                   person.get_latest_protection(t + dt / 2),
+                                                   false)); // Starting time in second order approximation
             }
         }
     }
@@ -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();