1279 minor additions to how infections work

cpp/models/abm/analyze_result.h

@@ -187,17 +187,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<ViralShedParameters>(num_runs),
                     [age_group, virus_variant](auto&& model) -> auto& {
-                        return model.parameters.template get<InfectivityDistributions>()[{virus_variant, age_group}];
+                        return model.parameters.template get<ViralShedDistribution>()[{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;

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
@@ -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);
 
     /**
@@ -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.
@@ -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.
@@ -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.
 };
 

cpp/models/abm/model_functions.cpp

@@ -32,21 +32,21 @@ namespace mio
 namespace abm
 {
 
-ScalarType daily_transmissions_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
-                                           const VirusVariant virus, const AgeGroup age_receiver,
-                                           const LocalInfectionParameters& params)
+ScalarType total_exposure_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
+                                      const VirusVariant virus, const AgeGroup age_receiver,
+                                      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) {
-        prob +=
+        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}];
 }
@@ -73,24 +73,26 @@ 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 (auto 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);
-                ScalarType local_indiv_trans_prob_v =
-                    (daily_transmissions_by_contacts(local_contact_exposure, cell_index, virus, 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_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
             }
         }
     }
@@ -110,14 +112,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();

cpp/models/abm/model_functions.h

@@ -32,28 +32,28 @@ namespace abm
 {
 
 /**
- * @brief Compute the number of daily transmissions for contact transmission of a virus in a cell.
+ * @brief Compute the total virus exposure for contact transmission in a cell (unit: per day).
  * @param[in] rates The local exposure rates.
  * @param[in] cell_index Cell index of the Cell.
  * @param[in] virus VirusVariant of interest.
  * @param[in] age_receiver AgeGroup of the receiving Person.
  * @param[in] params The local infection parameters.
- * @return Average amount of Infection%s with the virus from the AgeGroup of the transmitter per day.
+ * @return Total amount of virus exposure by contacts of the receiver.
  */
-ScalarType daily_transmissions_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
-                                           const VirusVariant virus, const AgeGroup age_receiver,
-                                           const LocalInfectionParameters& params);
+ScalarType total_exposure_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
+                                      const VirusVariant virus, const AgeGroup age_receiver,
+                                      const LocalInfectionParameters& params);
 
 /**
- * @brief Compute the number of daily transmissions for aerosol transmission of a virus in a cell.
+ * @brief Compute the total virus exposure for aerosol transmission in a cell (unit: per day).
  * @param[in] rates The local exposure rates.
  * @param[in] cell_index Cell index of the Cell.
  * @param[in] virus VirusVariant of interest.
  * @param[in] global_params The parameter set of the Model.
- * @return Average amount of Infection%s with the virus per day.
+ * @return Total amount of virus exposure by air.
  */
-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);
 
 /**
  * @brief Add the contribution of a person to the local exposure rates.