888 Make Apply_variant function more generic and fix issue with defau…

…lt value (#904) - Fix Bug, where the parameter for the transmission probability was constant zero - apply_variant function is now more generic and not specificied to the delta variant Co-authored-by: reneSchm <[email protected]>
SciCompMod · Jan 25, 2024 · 6391c30 · 6391c30
1 parent 1fa9e20
commit 6391c30
Show file tree

Hide file tree

Showing 8 changed files with 231 additions and 181 deletions.
diff --git a/cpp/benchmarks/flow_simulation_ode_secirvvs.h b/cpp/benchmarks/flow_simulation_ode_secirvvs.h
@@ -30,6 +30,9 @@ using FlowModel = osecirvvs::Model;
 // For comparison benchmarks, an old model version that does not provide computation of flows has been reimplemented here.
 // For more details see the original implementation in:
 // https://github.com/SciCompMod/memilio/blob/13555a6b23177d2d4633c393903461a27ce5762b/cpp/models/ode_secirvvs/model.h
+// Updates from Issue/PR 888:
+// - Apply_variant function has been adjusted to be more generic
+// - Fixed a bug where the transmission probability was always set to zero.
 
 class FlowlessModel : public CompartmentalModel<osecirvvs::InfectionState,
                                                 Populations<AgeGroup, osecirvvs::InfectionState>, osecirvvs::Parameters>
@@ -421,22 +424,23 @@ class Simulation : public Base
     {
     }
 
-    void apply_b161(double t)
+    void apply_variant(const double t, const CustomIndexArray<UncertainValue, AgeGroup> base_infectiousness)
     {
-
-        auto start_day   = this->get_model().parameters.template get<osecirvvs::StartDay>();
-        auto b161_growth = (start_day - get_day_in_year(Date(2021, 6, 6))) * 0.1666667;
-        // 2 equal to the share of the delta variant on June 6
-        double share_new_variant = std::min(1.0, pow(2, t * 0.1666667 + b161_growth) * 0.01);
-        size_t num_groups        = (size_t)this->get_model().parameters.get_num_groups();
-        for (size_t i = 0; i < num_groups; ++i) {
-            double new_transmission =
-                (1 - share_new_variant) *
-                    this->get_model().parameters.template get<osecirvvs::BaseInfectiousnessB117>()[(AgeGroup)i] +
-                share_new_variant *
-                    this->get_model().parameters.template get<osecirvvs::BaseInfectiousnessB161>()[(AgeGroup)i];
-            this->get_model().parameters.template get<osecirvvs::TransmissionProbabilityOnContact>()[(AgeGroup)i] =
-                new_transmission;
+        auto start_day             = this->get_model().parameters.template get<osecirvvs::StartDay>();
+        auto start_day_new_variant = this->get_model().parameters.template get<osecirvvs::StartDayNewVariant>();
+
+        if (start_day + t >= start_day_new_variant - 1e-10) {
+            const double days_variant      = t - (start_day_new_variant - start_day);
+            const double share_new_variant = std::min(1.0, 0.01 * pow(2, (1. / 7) * days_variant));
+            const auto num_groups          = this->get_model().parameters.get_num_groups();
+            for (auto i = AgeGroup(0); i < num_groups; ++i) {
+                double new_transmission =
+                    (1 - share_new_variant) * base_infectiousness[i] +
+                    share_new_variant * base_infectiousness[i] *
+                        this->get_model().parameters.template get<osecirvvs::InfectiousnessNewVariant>()[i];
+                this->get_model().parameters.template get<osecirvvs::TransmissionProbabilityOnContact>()[i] =
+                    new_transmission;
+            }
         }
     }
 
@@ -500,9 +504,13 @@ class Simulation : public Base
     */
     Eigen::Ref<Eigen::VectorXd> advance(double tmax)
     {
-        auto& t_end_dyn_npis   = this->get_model().parameters.get_end_dynamic_npis();
-        auto& dyn_npis         = this->get_model().parameters.template get<osecirvvs::DynamicNPIsInfectedSymptoms>();
-        auto& contact_patterns = this->get_model().parameters.template get<osecirvvs::ContactPatterns>();
+        auto& t_end_dyn_npis    = this->get_model().parameters.get_end_dynamic_npis();
+        auto& dyn_npis          = this->get_model().parameters.template get<osecirvvs::DynamicNPIsInfectedSymptoms>();
+        auto& contact_patterns  = this->get_model().parameters.template get<osecirvvs::ContactPatterns>();
+        const size_t num_groups = (size_t)this->get_model().parameters.get_num_groups();
+
+        auto base_infectiousness =
+            this->get_model().parameters.template get<osecirvvs::TransmissionProbabilityOnContact>();
 
         double delay_lockdown;
         auto t        = Base::get_result().get_last_time();
@@ -516,12 +524,12 @@ class Simulation : public Base
 
             if (t == 0) {
                 //this->apply_vaccination(t); // done in init now?
-                this->apply_b161(t);
+                this->apply_variant(t, base_infectiousness);
             }
             Base::advance(t + dt_eff);
             if (t + 0.5 + dt_eff - std::floor(t + 0.5) >= 1) {
                 this->apply_vaccination(t + 0.5 + dt_eff);
-                this->apply_b161(t);
+                this->apply_variant(t, base_infectiousness);
             }
 
             if (t > 0) {
@@ -560,6 +568,9 @@ class Simulation : public Base
                 m_t_last_npi_check = t;
             }
         }
+
+        this->get_model().parameters.template get<osecirvvs::TransmissionProbabilityOnContact>() = base_infectiousness;
+
         return this->get_result().get_last_value();
     }
 

diff --git a/cpp/models/ode_secirvvs/analyze_result.h b/cpp/models/ode_secirvvs/analyze_result.h
@@ -165,11 +165,7 @@ std::vector<Model> ensemble_params_percentile(const std::vector<std::vector<Mode
             //virus variants
             param_percentil(
                 node, [i](auto&& model) -> auto& {
-                    return model.parameters.template get<BaseInfectiousnessB161>()[i];
-                });
-            param_percentil(
-                node, [i](auto&& model) -> auto& {
-                    return model.parameters.template get<BaseInfectiousnessB117>()[i];
+                    return model.parameters.template get<InfectiousnessNewVariant>()[i];
                 });
         }
         // group independent params

diff --git a/cpp/models/ode_secirvvs/model.h b/cpp/models/ode_secirvvs/model.h
@@ -543,21 +543,34 @@ class Simulation : public BaseT
     {
     }
 
-    void apply_b161(double t)
-    {
+    /**
+    * @brief Applies the effect of a new variant of a disease to the transmission probability of the model.
+    * 
+    * This function adjusts the transmission probability of the disease for each age group based on the share of the new variant.
+    * The share of the new variant is calculated based on the time `t` and the start day of the new variant.
+    * The transmission probability is then updated for each age group in the model.
+    * 
+    * Based on Equation (35) and (36) in doi.org/10.1371/journal.pcbi.1010054
+    * 
+    * @param [in] t The current time.
+    * @param [in] base_infectiousness The base infectiousness of the old variant for each age group.
+    */
 
-        auto start_day   = this->get_model().parameters.template get<StartDay>();
-        auto b161_growth = (start_day - get_day_in_year(Date(2021, 6, 6))) * 0.1666667;
-        // 2 equal to the share of the delta variant on June 6
-        double share_new_variant = std::min(1.0, pow(2, t * 0.1666667 + b161_growth) * 0.01);
-        size_t num_groups        = (size_t)this->get_model().parameters.get_num_groups();
-        for (size_t i = 0; i < num_groups; ++i) {
-            double new_transmission =
-                (1 - share_new_variant) *
-                    this->get_model().parameters.template get<BaseInfectiousnessB117>()[(AgeGroup)i] +
-                share_new_variant * this->get_model().parameters.template get<BaseInfectiousnessB161>()[(AgeGroup)i];
-            this->get_model().parameters.template get<TransmissionProbabilityOnContact>()[(AgeGroup)i] =
-                new_transmission;
+    void apply_variant(const double t, const CustomIndexArray<UncertainValue, AgeGroup> base_infectiousness)
+    {
+        auto start_day             = this->get_model().parameters.template get<StartDay>();
+        auto start_day_new_variant = this->get_model().parameters.template get<StartDayNewVariant>();
+
+        if (start_day + t >= start_day_new_variant - 1e-10) {
+            const double days_variant      = t - (start_day_new_variant - start_day);
+            const double share_new_variant = std::min(1.0, 0.01 * pow(2, (1. / 7) * days_variant));
+            const auto num_groups          = this->get_model().parameters.get_num_groups();
+            for (auto i = AgeGroup(0); i < num_groups; ++i) {
+                double new_transmission = (1 - share_new_variant) * base_infectiousness[i] +
+                                          share_new_variant * base_infectiousness[i] *
+                                              this->get_model().parameters.template get<InfectiousnessNewVariant>()[i];
+                this->get_model().parameters.template get<TransmissionProbabilityOnContact>()[i] = new_transmission;
+            }
         }
     }
 
@@ -622,6 +635,11 @@ class Simulation : public BaseT
         auto& t_end_dyn_npis   = this->get_model().parameters.get_end_dynamic_npis();
         auto& dyn_npis         = this->get_model().parameters.template get<DynamicNPIsInfectedSymptoms>();
         auto& contact_patterns = this->get_model().parameters.template get<ContactPatterns>();
+        // const size_t num_groups = (size_t)this->get_model().parameters.get_num_groups();
+
+        // in the apply_variant function, we adjust the TransmissionProbabilityOnContact parameter. We need to store
+        // the base value to use it in the apply_variant function and also to reset the parameter after the simulation.
+        auto base_infectiousness = this->get_model().parameters.template get<TransmissionProbabilityOnContact>();
 
         double delay_lockdown;
         auto t        = BaseT::get_result().get_last_time();
@@ -635,12 +653,12 @@ class Simulation : public BaseT
 
             if (t == 0) {
                 //this->apply_vaccination(t); // done in init now?
-                this->apply_b161(t);
+                this->apply_variant(t, base_infectiousness);
             }
             BaseT::advance(t + dt_eff);
             if (t + 0.5 + dt_eff - std::floor(t + 0.5) >= 1) {
                 this->apply_vaccination(t + 0.5 + dt_eff);
-                this->apply_b161(t);
+                this->apply_variant(t, base_infectiousness);
             }
 
             if (t > 0) {
@@ -679,6 +697,10 @@ class Simulation : public BaseT
                 m_t_last_npi_check = t;
             }
         }
+        // reset TransmissionProbabilityOnContact. This is important for the graph simulation where the advance
+        // function is called multiple times for the same model.
+        this->get_model().parameters.template get<TransmissionProbabilityOnContact>() = base_infectiousness;
+
         return this->get_result().get_last_value();
     }
 

diff --git a/cpp/models/ode_secirvvs/parameter_space.cpp b/cpp/models/ode_secirvvs/parameter_space.cpp
@@ -149,10 +149,7 @@ Graph<Model, MigrationParameters> draw_sample(Graph<Model, MigrationParameters>&
 
     //infectiousness of virus variants is not sampled independently but depend on base infectiousness
     for (auto i = AgeGroup(0); i < shared_params_model.parameters.get_num_groups(); ++i) {
-        shared_params_model.parameters.template get<BaseInfectiousnessB117>()[i] =
-            shared_params_model.parameters.template get<TransmissionProbabilityOnContact>()[i];
-        shared_params_model.parameters.template get<BaseInfectiousnessB161>()[i] =
-            shared_params_model.parameters.template get<TransmissionProbabilityOnContact>()[i] * delta_fac;
+        shared_params_model.parameters.template get<InfectiousnessNewVariant>()[i] = delta_fac;
     }
 
     for (auto& params_node : graph.nodes()) {