1547 (LCT) implement possibility of death from hospitalization without intensive care

cpp/examples/glct_secir.cpp

@@ -40,7 +40,7 @@ int main()
     // We need to double the choices of the number of subcompartments for some compartments,
     // as we define different strains for different transition possibilities. For both strains, the same number of
     // subcompartments is chosen. The transition probabilities are defined in the StartingProbabilities.
-    constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 6, NumInfectedSymptoms = 2, NumInfectedSevere = 2,
+    constexpr size_t NumExposed = 2, NumInfectedNoSymptoms = 6, NumInfectedSymptoms = 2, NumInfectedSevere = 3,
                      NumInfectedCritical = 10;
     using Model          = mio::glsecir::Model<ScalarType, NumExposed, NumInfectedNoSymptoms, NumInfectedSymptoms,
                                                NumInfectedSevere, NumInfectedCritical>;
@@ -62,6 +62,7 @@ int main()
     const ScalarType recoveredPerInfectedNoSymptoms = 0.09;
     const ScalarType severePerInfectedSymptoms      = 0.2;
     const ScalarType criticalPerSevere              = 0.25;
+    const ScalarType deathsPerSevere                = 0.1;
     const ScalarType deathsPerCritical              = 0.3;
 
     // Define the initial values with the distribution of the population into subcompartments.
@@ -77,7 +78,7 @@ int main()
          10 * (1 - recoveredPerInfectedNoSymptoms), 20 * recoveredPerInfectedNoSymptoms,
          10 * recoveredPerInfectedNoSymptoms, 10 * recoveredPerInfectedNoSymptoms},
         {50 * severePerInfectedSymptoms, 50 * (1 - severePerInfectedSymptoms)},
-        {50 * criticalPerSevere, 50 * (1 - criticalPerSevere)},
+        {50 * criticalPerSevere, 50 * deathsPerSevere, 50 * (1 - criticalPerSevere - deathsPerSevere)},
         {10 * deathsPerCritical, 10 * deathsPerCritical, 5 * deathsPerCritical, 3 * deathsPerCritical,
          2 * deathsPerCritical, 10 * (1 - deathsPerCritical), 10 * (1 - deathsPerCritical), 5 * (1 - deathsPerCritical),
          3 * (1 - deathsPerCritical), 2 * (1 - deathsPerCritical)},
@@ -136,9 +137,8 @@ int main()
     Eigen::VectorX<ScalarType> StartingProbabilitiesInfectedNoSymptoms =
         Eigen::VectorX<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>());
     StartingProbabilitiesInfectedNoSymptoms[0] = 1 - recoveredPerInfectedNoSymptoms;
-    StartingProbabilitiesInfectedNoSymptoms[(
-        Eigen::Index)(LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>() / 2.)] =
-        recoveredPerInfectedNoSymptoms;
+    StartingProbabilitiesInfectedNoSymptoms[(Eigen::Index)(
+        LctState::get_num_subcompartments<InfectionState::InfectedNoSymptoms>() / 2.)] = recoveredPerInfectedNoSymptoms;
     model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedNoSymptoms<ScalarType>>() =
         StartingProbabilitiesInfectedNoSymptoms;
     // Define equal TransitionMatrices for the strains.
@@ -155,10 +155,9 @@ int main()
     // InfectedSymptoms.
     Eigen::VectorX<ScalarType> StartingProbabilitiesInfectedSymptoms =
         Eigen::VectorX<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>());
-    StartingProbabilitiesInfectedSymptoms[0] = severePerInfectedSymptoms;
-    StartingProbabilitiesInfectedSymptoms[(
-        Eigen::Index)(LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>() / 2.)] =
-        1 - severePerInfectedSymptoms;
+    StartingProbabilitiesInfectedSymptoms[0]                                         = severePerInfectedSymptoms;
+    StartingProbabilitiesInfectedSymptoms[(Eigen::Index)(
+        LctState::get_num_subcompartments<InfectionState::InfectedSymptoms>() / 2.)] = 1 - severePerInfectedSymptoms;
     model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedSymptoms<ScalarType>>() =
         StartingProbabilitiesInfectedSymptoms;
     model.parameters.get<mio::glsecir::TransitionMatrixInfectedSymptomsToInfectedSevere<ScalarType>>() =
@@ -170,25 +169,29 @@ int main()
     // InfectedSevere.
     Eigen::VectorX<ScalarType> StartingProbabilitiesInfectedSevere =
         Eigen::VectorX<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedSevere>());
-    StartingProbabilitiesInfectedSevere[0] = criticalPerSevere;
-    StartingProbabilitiesInfectedSevere[(
-        Eigen::Index)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.)] =
-        1 - criticalPerSevere;
+    StartingProbabilitiesInfectedSevere[0]                                         = criticalPerSevere;
+    StartingProbabilitiesInfectedSevere[(Eigen::Index)(
+        LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 3.)] = deathsPerSevere;
+    StartingProbabilitiesInfectedSevere[2 * (Eigen::Index)(
+                                                LctState::get_num_subcompartments<InfectionState::InfectedSevere>() /
+                                                3.)] = 1 - criticalPerSevere - deathsPerSevere;
     model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedSevere<ScalarType>>() =
         StartingProbabilitiesInfectedSevere;
     model.parameters.get<mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical<ScalarType>>() =
         mio::glsecir::TransitionMatrixInfectedSevereToInfectedCritical<ScalarType>().get_default(
-            (size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.), timeInfectedSevere);
+            (size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 3.), timeInfectedSevere);
+    model.parameters.get<mio::glsecir::TransitionMatrixInfectedSevereToDead<ScalarType>>() =
+        mio::glsecir::TransitionMatrixInfectedSevereToDead<ScalarType>().get_default(
+            (size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 3.), timeInfectedSevere);
     model.parameters.get<mio::glsecir::TransitionMatrixInfectedSevereToRecovered<ScalarType>>() =
         mio::glsecir::TransitionMatrixInfectedSevereToRecovered<ScalarType>().get_default(
-            (size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 2.), timeInfectedSevere);
+            (size_t)(LctState::get_num_subcompartments<InfectionState::InfectedSevere>() / 3.), timeInfectedSevere);
     // InfectedCritical.
     Eigen::VectorX<ScalarType> StartingProbabilitiesInfectedCritical =
         Eigen::VectorX<ScalarType>::Zero(LctState::get_num_subcompartments<InfectionState::InfectedCritical>());
-    StartingProbabilitiesInfectedCritical[0] = deathsPerCritical;
-    StartingProbabilitiesInfectedCritical[(
-        Eigen::Index)(LctState::get_num_subcompartments<InfectionState::InfectedCritical>() / 2.)] =
-        1 - deathsPerCritical;
+    StartingProbabilitiesInfectedCritical[0]                                         = deathsPerCritical;
+    StartingProbabilitiesInfectedCritical[(Eigen::Index)(
+        LctState::get_num_subcompartments<InfectionState::InfectedCritical>() / 2.)] = 1 - deathsPerCritical;
     model.parameters.get<mio::glsecir::StartingProbabilitiesInfectedCritical<ScalarType>>() =
         StartingProbabilitiesInfectedCritical;
     model.parameters.get<mio::glsecir::TransitionMatrixInfectedCriticalToDead<ScalarType>>() =

cpp/examples/lct_secir.cpp

@@ -71,6 +71,7 @@ int main()
     model.parameters.get<mio::lsecir::RecoveredPerInfectedNoSymptoms<ScalarType>>()[0] = 0.09;
     model.parameters.get<mio::lsecir::SeverePerInfectedSymptoms<ScalarType>>()[0]      = 0.2;
     model.parameters.get<mio::lsecir::CriticalPerSevere<ScalarType>>()[0]              = 0.25;
+    model.parameters.get<mio::lsecir::DeathsPerSevere<ScalarType>>()[0]                = 0.1;
     model.parameters.get<mio::lsecir::DeathsPerCritical<ScalarType>>()[0]              = 0.3;
 
     if (use_initializer_flows) {

cpp/examples/lct_secir_2_diseases.cpp

@@ -89,6 +89,8 @@ int main()
     model.parameters.get<mio::lsecir2d::SeverePerInfectedSymptoms_b<ScalarType>>()[0]      = 0.2;
     model.parameters.get<mio::lsecir2d::CriticalPerSevere_a<ScalarType>>()[0]              = 0.25;
     model.parameters.get<mio::lsecir2d::CriticalPerSevere_b<ScalarType>>()[0]              = 0.25;
+    model.parameters.get<mio::lsecir2d::DeathsPerSevere_a<ScalarType>>()[0]                = 0.1;
+    model.parameters.get<mio::lsecir2d::DeathsPerSevere_b<ScalarType>>()[0]                = 0.1;
     model.parameters.get<mio::lsecir2d::DeathsPerCritical_a<ScalarType>>()[0]              = 0.3;
     model.parameters.get<mio::lsecir2d::DeathsPerCritical_b<ScalarType>>()[0]              = 0.3;
 

cpp/models/glct_secir/model.h

@@ -268,23 +268,32 @@ class Model
             params.template get<TransitionMatrixInfectedSevereToInfectedCritical<FP>>().transpose() *
             y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>(),
                       dimensionInfectedSevereToInfectedCritical);
+        // Flow from InfectedSevere To Dead.
+        size_t dimensionInfectedSevereToDead = params.template get<TransitionMatrixInfectedSevereToDead<FP>>().rows();
+        dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
+                         dimensionInfectedSevereToInfectedCritical,
+                     dimensionInfectedSevereToDead) +=
+            params.template get<TransitionMatrixInfectedSevereToDead<FP>>().transpose() *
+            y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
+                          dimensionInfectedSevereToInfectedCritical,
+                      dimensionInfectedSevereToDead);
         // Flow from InfectedSevere To Recovered.
         size_t dimensionInfectedSevereToRecovered =
             params.template get<TransitionMatrixInfectedSevereToRecovered<FP>>().rows();
         dydt.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
-                         dimensionInfectedSevereToInfectedCritical,
+                         dimensionInfectedSevereToInfectedCritical + dimensionInfectedSevereToDead,
                      dimensionInfectedSevereToRecovered) +=
             params.template get<TransitionMatrixInfectedSevereToRecovered<FP>>().transpose() *
             y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
-                          dimensionInfectedSevereToInfectedCritical,
+                          dimensionInfectedSevereToInfectedCritical + dimensionInfectedSevereToDead,
                       dimensionInfectedSevereToRecovered);
         // Add flow directly to Recovered compartment.
         dydt[LctState::template get_first_index<InfectionState::Recovered>()] +=
             -(params.template get<TransitionMatrixInfectedSevereToRecovered<FP>>() *
               Eigen::VectorX<FP>::Ones(dimensionInfectedSevereToRecovered))
                  .transpose() *
             y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
-                          dimensionInfectedSevereToInfectedCritical,
+                          dimensionInfectedSevereToInfectedCritical + dimensionInfectedSevereToDead,
                       dimensionInfectedSevereToRecovered);
 
         // --- InfectedCritical. ---
@@ -325,7 +334,14 @@ class Model
                       dimensionInfectedCriticalToRecovered);
 
         // --- Dead. ---
-        dydt[LctState::template get_first_index<InfectionState::Dead>()] =
+        dydt[LctState::template get_first_index<InfectionState::Dead>()] +=
+            -(params.template get<TransitionMatrixInfectedSevereToDead<FP>>() *
+              Eigen::VectorX<FP>::Ones(dimensionInfectedSevereToDead))
+                 .transpose() *
+            y.segment(LctState::template get_first_index<InfectionState::InfectedSevere>() +
+                          dimensionInfectedSevereToInfectedCritical,
+                      dimensionInfectedSevereToDead);
+        dydt[LctState::template get_first_index<InfectionState::Dead>()] +=
             -(params.template get<TransitionMatrixInfectedCriticalToDead<FP>>() *
               Eigen::VectorX<FP>::Ones(dimensionInfectedCriticalToDead))
                  .transpose() *

cpp/models/glct_secir/parameters.h

@@ -263,6 +263,31 @@ struct TransitionMatrixInfectedSevereToInfectedCritical {
     }
 };
 
+/**
+ * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedSevere
+ *      compartment before dying.
+ */
+template <typename FP>
+struct TransitionMatrixInfectedSevereToDead {
+    using Type = Eigen::MatrixX<FP>;
+    /**
+     * @brief Default parameters can be used to get an Erlang distributed stay time in InfectedSevere compartment
+     *   before dying.
+     * @param[in] dimension Number of rows/columns of the transition matrix.
+     * @param[in] time Average time spent in InfectedSevere before dying in day unit.
+     */
+    static Type get_default(size_t dimension, FP time = 1.)
+    {
+        Eigen::MatrixX<FP> def = Eigen::VectorX<FP>::Constant(dimension, -(FP)dimension / time).asDiagonal();
+        def.diagonal(1).setConstant((FP)dimension / time);
+        return def;
+    }
+    static std::string name()
+    {
+        return "TransitionMatrixInfectedSevereToDead";
+    }
+};
+
 /**
  * @brief Transition matrix of the phase-type distribution describing the stay time in the InfectedSevere
  *      compartment before recovery.
@@ -461,10 +486,11 @@ using ParametersBase =
                  TransitionMatrixInfectedNoSymptomsToRecovered<FP>, StartingProbabilitiesInfectedSymptoms<FP>,
                  TransitionMatrixInfectedSymptomsToInfectedSevere<FP>, TransitionMatrixInfectedSymptomsToRecovered<FP>,
                  StartingProbabilitiesInfectedSevere<FP>, TransitionMatrixInfectedSevereToInfectedCritical<FP>,
-                 TransitionMatrixInfectedSevereToRecovered<FP>, StartingProbabilitiesInfectedCritical<FP>,
-                 TransitionMatrixInfectedCriticalToDead<FP>, TransitionMatrixInfectedCriticalToRecovered<FP>,
-                 TransmissionProbabilityOnContact<FP>, ContactPatterns<FP>, RelativeTransmissionNoSymptoms<FP>,
-                 RiskOfInfectionFromSymptomatic<FP>, StartDay<FP>, Seasonality<FP>>;
+                 TransitionMatrixInfectedSevereToDead<FP>, TransitionMatrixInfectedSevereToRecovered<FP>,
+                 StartingProbabilitiesInfectedCritical<FP>, TransitionMatrixInfectedCriticalToDead<FP>,
+                 TransitionMatrixInfectedCriticalToRecovered<FP>, TransmissionProbabilityOnContact<FP>,
+                 ContactPatterns<FP>, RelativeTransmissionNoSymptoms<FP>, RiskOfInfectionFromSymptomatic<FP>,
+                 StartDay<FP>, Seasonality<FP>>;
 
 /// @brief Parameters of an GLCT-SECIR model.
 template <typename FP>
@@ -523,6 +549,8 @@ class Parameters : public ParametersBase<FP>
              this->template get<TransitionMatrixInfectedSymptomsToRecovered<FP>>().rows()) ||
             (this->template get<TransitionMatrixInfectedSevereToInfectedCritical<FP>>().cols() !=
              this->template get<TransitionMatrixInfectedSevereToInfectedCritical<FP>>().rows()) ||
+            (this->template get<TransitionMatrixInfectedSevereToDead<FP>>().cols() !=
+             this->template get<TransitionMatrixInfectedSevereToDead<FP>>().rows()) ||
             (this->template get<TransitionMatrixInfectedSevereToRecovered<FP>>().cols() !=
              this->template get<TransitionMatrixInfectedSevereToRecovered<FP>>().rows()) ||
             (this->template get<TransitionMatrixInfectedCriticalToDead<FP>>().cols() !=
@@ -559,6 +587,7 @@ class Parameters : public ParametersBase<FP>
 
         if (this->template get<StartingProbabilitiesInfectedSevere<FP>>().rows() !=
             this->template get<TransitionMatrixInfectedSevereToInfectedCritical<FP>>().rows() +
+                this->template get<TransitionMatrixInfectedSevereToDead<FP>>().rows() +
                 this->template get<TransitionMatrixInfectedSevereToRecovered<FP>>().rows()) {
             log_error("Constraint check: Dimensions of StartingProbabilitiesInfectedSevere and "
                       "TransitionMatrices of InfectedSevere compartment are not matching.");
@@ -653,6 +682,14 @@ class Parameters : public ParametersBase<FP>
                 "flow InfectedSevereToInfectedCritical.");
             return true;
         }
+        if (((this->template get<TransitionMatrixInfectedSevereToDead<FP>>() *
+              Eigen::VectorX<FP>::Ones(this->template get<TransitionMatrixInfectedSevereToDead<FP>>().rows()))
+                 .array() > 1e-10)
+                .any()) {
+            log_warning("Constraint check: The entries of TransitionMatrixInfectedSevereToDead lead to a negative "
+                        "flow InfectedSevereToDead.");
+            return true;
+        }
         if (((this->template get<TransitionMatrixInfectedSevereToRecovered<FP>>() *
               Eigen::VectorX<FP>::Ones(this->template get<TransitionMatrixInfectedSevereToRecovered<FP>>().rows()))
                  .array() > 1e-10)

cpp/models/ide_secir/parameters_io.h

@@ -56,9 +56,6 @@ namespace isecir
 * using the means of the respective TransitionDistribution. 
 * The flow InfectedNoSymptomsToInfectedSymptoms is calculated with the standard formula from the IDE model
 * using the results for ExposedToInfectedNoSymptoms.
-* Throughout this RKI-based initialization, we assume that individuals can only die from InfectedCritical and not from 
-* InfectedSevere, i.e. this initialization routine assumes that the probability from transitioning from InfectedSevere 
-* to Dead is 0. 
 *
 * The number of deaths used in the model is computed by shifting the reported RKI data according to the mean values of the transitions 
 * InfectedSymptomsToInfectedSevere, InfectedSevereToInfectedCritical and InfectedCriticalToDead.

cpp/models/lct_secir/infection_state.h

@@ -56,10 +56,11 @@ enum class InfectionTransition
     InfectedSymptomsToInfectedSevere     = 4,
     InfectedSymptomsToRecovered          = 5,
     InfectedSevereToInfectedCritical     = 6,
-    InfectedSevereToRecovered            = 7,
-    InfectedCriticalToDead               = 8,
-    InfectedCriticalToRecovered          = 9,
-    Count                                = 10
+    InfectedSevereToDead                 = 7,
+    InfectedSevereToRecovered            = 8,
+    InfectedCriticalToDead               = 9,
+    InfectedCriticalToRecovered          = 10,
+    Count                                = 11
 };
 
 } // namespace lsecir

cpp/models/lct_secir/model.h

@@ -260,7 +260,9 @@ class Model : public CompartmentalModel<FP, InfectionState, LctPopulations<FP, L
         }
 
         // Calculate derivative of the InfectedCritical compartment.
-        dydt[Ri] += dydt[ICri_first_index] * (1 - params.template get<CriticalPerSevere<FP>>()[Group]);
+        dydt[Ri] += dydt[ICri_first_index] * (1 - params.template get<CriticalPerSevere<FP>>()[Group] -
+                                              params.template get<DeathsPerSevere<FP>>()[Group]);
+        dydt[Di]               = dydt[ICri_first_index] * params.template get<DeathsPerSevere<FP>>()[Group];
         dydt[ICri_first_index] = dydt[ICri_first_index] * params.template get<CriticalPerSevere<FP>>()[Group];
         for (size_t subcomp = 0;
              subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() - 1;
@@ -276,7 +278,7 @@ class Model : public CompartmentalModel<FP, InfectionState, LctPopulations<FP, L
                (1 / params.template get<TimeInfectedCritical<FP>>()[Group]) * y[Ri - 1];
         dydt[Ri - 1] -= flow;
         dydt[Ri] = dydt[Ri] + (1 - params.template get<DeathsPerCritical<FP>>()[Group]) * flow;
-        dydt[Di] = params.template get<DeathsPerCritical<FP>>()[Group] * flow;
+        dydt[Di] = dydt[Di] + params.template get<DeathsPerCritical<FP>>()[Group] * flow;
 
         // Function call for next group if applicable.
         if constexpr (Group + 1 < num_groups) {