Skip to content

Commit c18828e

Browse files
johapauHenrZu
johapau
and
HenrZu
authored
677 Computation of ode secir model reproduction number (#858)
Co-authored-by: HenrZu <henrik.zunker@dlr.de>
1 parent 3b8cf9d commit c18828e

File tree

2 files changed

+426
-0
lines changed

2 files changed

+426
-0
lines changed

‎cpp/models/ode_secir/model.h‎

Lines changed: 259 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
#include "ode_secir/parameters.h"
2828
#include "memilio/math/smoother.h"
2929
#include "memilio/math/eigen_util.h"
30+
#include "memilio/math/interpolation.h"
3031

3132
namespace mio
3233
{
@@ -328,6 +329,264 @@ double get_infections_relative(const Simulation<Base>& sim, double /*t*/, const
328329
return inf_rel;
329330
}
330331

332+
/**
333+
*@brief Computes the reproduction number at a given index time of the Model output obtained by the Simulation.
334+
*@param t_idx The index time at which the reproduction number is computed.
335+
*@param sim The simulation holding the SECIR model
336+
*@tparam Base simulation type that uses a SECIR compartment model. see Simulation.
337+
*@returns The computed reproduction number at the provided index time.
338+
*/
339+
340+
template <class Base>
341+
IOResult<ScalarType> get_reproduction_number(size_t t_idx, const Simulation<Base>& sim)
342+
{
343+
344+
if (!(t_idx < static_cast<size_t>(sim.get_result().get_num_time_points()))) {
345+
return mio::failure(mio::StatusCode::OutOfRange, "t_idx is not a valid index for the TimeSeries");
346+
}
347+
348+
auto const& params = sim.get_model().parameters;
349+
const size_t num_groups = (size_t)sim.get_model().parameters.get_num_groups();
350+
//The infected compartments in the SECIR Model are: Exposed, Carrier, Infected, Hospitalized, ICU in respective agegroups
351+
const size_t num_infected_compartments = 5;
352+
const size_t total_infected_compartments = num_infected_compartments * num_groups;
353+
const double pi = std::acos(-1);
354+
355+
//F encodes new Infections and V encodes transition times in the next-generation matrix calculation of R_t
356+
Eigen::MatrixXd F(total_infected_compartments, total_infected_compartments);
357+
Eigen::MatrixXd V(total_infected_compartments, total_infected_compartments);
358+
F = Eigen::MatrixXd::Zero(total_infected_compartments,
359+
total_infected_compartments); //Initialize matrices F and V with zeroes
360+
V = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments);
361+
362+
auto test_and_trace_required = 0.0;
363+
auto icu_occupancy = 0.0;
364+
for (auto i = AgeGroup(0); i < (mio::AgeGroup)num_groups; ++i) {
365+
auto rateINS = 0.5 / (params.template get<IncubationTime>()[i] - params.template get<SerialInterval>()[i]);
366+
test_and_trace_required +=
367+
(1 - params.template get<RecoveredPerInfectedNoSymptoms>()[i]) * rateINS *
368+
sim.get_result().get_value(
369+
t_idx)[sim.get_model().populations.get_flat_index({i, InfectionState::InfectedNoSymptoms})];
370+
icu_occupancy += sim.get_result().get_value(
371+
t_idx)[sim.get_model().populations.get_flat_index({i, InfectionState::InfectedCritical})];
372+
}
373+
374+
double season_val = (1 + params.template get<Seasonality>() *
375+
sin(pi * (std::fmod((sim.get_model().parameters.template get<StartDay>() +
376+
sim.get_result().get_time(t_idx)),
377+
365.0) /
378+
182.5 +
379+
0.5)));
380+
ContactMatrixGroup const& contact_matrix = sim.get_model().parameters.template get<ContactPatterns>();
381+
382+
Eigen::MatrixXd cont_freq_eff(num_groups, num_groups);
383+
Eigen::MatrixXd riskFromInfectedSymptomatic_derivatives(num_groups, num_groups);
384+
Eigen::VectorXd divN(num_groups);
385+
Eigen::VectorXd riskFromInfectedSymptomatic(num_groups);
386+
Eigen::VectorXd rateINS(num_groups);
387+
388+
for (mio::AgeGroup k = 0; k < (mio::AgeGroup)num_groups; k++) {
389+
double temp = sim.get_result().get_value(
390+
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Susceptible})] +
391+
sim.get_result().get_value(
392+
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Exposed})] +
393+
sim.get_result().get_value(
394+
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedNoSymptoms})] +
395+
sim.get_result().get_value(
396+
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedSymptoms})] +
397+
sim.get_result().get_value(
398+
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedSevere})] +
399+
sim.get_result().get_value(
400+
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedCritical})] +
401+
sim.get_result().get_value(
402+
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Recovered})];
403+
if (temp == 0) {
404+
temp = 1;
405+
}
406+
divN[(size_t)k] = 1 / temp;
407+
408+
riskFromInfectedSymptomatic[(size_t)k] = smoother_cosine(
409+
test_and_trace_required, params.template get<TestAndTraceCapacity>(),
410+
(params.template get<TestAndTraceCapacity>()) * 5, params.template get<RiskOfInfectionFromSymptomatic>()[k],
411+
params.template get<MaxRiskOfInfectionFromSymptomatic>()[k]);
412+
413+
rateINS[(size_t)k] =
414+
0.5 / (params.template get<IncubationTime>()[k] - params.template get<SerialInterval>()[(mio::AgeGroup)k]);
415+
416+
for (mio::AgeGroup l = 0; l < (mio::AgeGroup)num_groups; l++) {
417+
if (test_and_trace_required < params.template get<TestAndTraceCapacity>() ||
418+
test_and_trace_required > 5 * params.template get<TestAndTraceCapacity>()) {
419+
riskFromInfectedSymptomatic_derivatives((size_t)k, (size_t)l) = 0;
420+
}
421+
else {
422+
riskFromInfectedSymptomatic_derivatives((size_t)k, (size_t)l) =
423+
-0.5 * pi *
424+
(params.template get<MaxRiskOfInfectionFromSymptomatic>()[k] -
425+
params.template get<RiskOfInfectionFromSymptomatic>()[k]) /
426+
(4 * params.template get<TestAndTraceCapacity>()) *
427+
(1 - params.template get<RecoveredPerInfectedNoSymptoms>()[l]) * rateINS[(size_t)l] *
428+
std::sin(pi / (4 * params.template get<TestAndTraceCapacity>()) *
429+
(test_and_trace_required - params.template get<TestAndTraceCapacity>()));
430+
}
431+
}
432+
433+
for (Eigen::Index l = 0; l < (Eigen::Index)num_groups; l++) {
434+
cont_freq_eff(l, (size_t)k) =
435+
season_val * contact_matrix.get_matrix_at(static_cast<double>(t_idx))(
436+
static_cast<Eigen::Index>((size_t)l), static_cast<Eigen::Index>((size_t)k));
437+
}
438+
}
439+
440+
//Check criterion if matrix V will be invertible by checking if subblock J is invertible
441+
Eigen::MatrixXd J(num_groups, num_groups);
442+
J = Eigen::MatrixXd::Zero(num_groups, num_groups);
443+
for (size_t i = 0; i < num_groups; i++) {
444+
J(i, i) = 1 / (params.template get<TimeInfectedCritical>()[(mio::AgeGroup)i]);
445+
446+
if (!(icu_occupancy < 0.9 * params.template get<ICUCapacity>() ||
447+
icu_occupancy > (double)(params.template get<ICUCapacity>()))) {
448+
for (size_t j = 0; j < num_groups; j++) {
449+
J(i, j) -= sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
450+
{(mio::AgeGroup)i, InfectionState::InfectedSevere})] /
451+
params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i] * 5 *
452+
params.template get<CriticalPerSevere>()[(mio::AgeGroup)i] * pi /
453+
(params.template get<ICUCapacity>()) *
454+
std::sin(pi / (0.1 * params.template get<ICUCapacity>()) *
455+
(icu_occupancy - 0.9 * params.template get<ICUCapacity>()));
456+
}
457+
}
458+
}
459+
460+
//Check, if J is invertible
461+
if (J.determinant() == 0) {
462+
return mio::failure(mio::StatusCode::UnknownError, "Matrix V is not invertible");
463+
}
464+
465+
//Initialize the matrix F
466+
for (size_t i = 0; i < num_groups; i++) {
467+
for (size_t j = 0; j < num_groups; j++) {
468+
469+
double temp = 0;
470+
for (Eigen::Index k = 0; k < (Eigen::Index)num_groups; k++) {
471+
temp += cont_freq_eff(i, k) *
472+
sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
473+
{(mio::AgeGroup)k, InfectionState::InfectedSymptoms})] *
474+
riskFromInfectedSymptomatic_derivatives(k, j) * divN[k];
475+
}
476+
477+
F(i, j + num_groups) =
478+
sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
479+
{(mio::AgeGroup)i, InfectionState::Susceptible})] *
480+
params.template get<TransmissionProbabilityOnContact>()[(mio::AgeGroup)i] *
481+
(cont_freq_eff(i, j) * params.template get<RelativeTransmissionNoSymptoms>()[(mio::AgeGroup)j] *
482+
divN[(size_t)j] +
483+
temp);
484+
}
485+
486+
for (size_t j = 0; j < num_groups; j++) {
487+
F(i, j + 2 * num_groups) = sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
488+
{(mio::AgeGroup)i, InfectionState::Susceptible})] *
489+
params.template get<TransmissionProbabilityOnContact>()[(mio::AgeGroup)i] *
490+
cont_freq_eff(i, j) * riskFromInfectedSymptomatic[(size_t)j] * divN[(size_t)j];
491+
}
492+
}
493+
494+
//Initialize the matrix V
495+
for (Eigen::Index i = 0; i < (Eigen::Index)num_groups; i++) {
496+
497+
double rateE = 1.0 / (2 * params.template get<SerialInterval>()[(mio::AgeGroup)i] -
498+
params.template get<IncubationTime>()[(mio::AgeGroup)i]);
499+
500+
double criticalPerSevereAdjusted = smoother_cosine(
501+
icu_occupancy, 0.90 * params.template get<ICUCapacity>(), params.template get<ICUCapacity>(),
502+
params.template get<CriticalPerSevere>()[(mio::AgeGroup)i], 0);
503+
504+
V(i, i) = rateE;
505+
V(i + num_groups, i) = -rateE;
506+
V(i + num_groups, i + num_groups) = rateINS[i];
507+
V(i + 2 * num_groups, i + num_groups) =
508+
-(1 - params.template get<RecoveredPerInfectedNoSymptoms>()[(mio::AgeGroup)i]) * rateINS[i];
509+
V(i + 2 * num_groups, i + 2 * num_groups) = (1 / params.template get<TimeInfectedSymptoms>()[(mio::AgeGroup)i]);
510+
V(i + 3 * num_groups, i + 2 * num_groups) =
511+
-params.template get<SeverePerInfectedSymptoms>()[(mio::AgeGroup)i] /
512+
params.template get<TimeInfectedSymptoms>()[(mio::AgeGroup)i];
513+
V(i + 3 * num_groups, i + 3 * num_groups) = 1 / (params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i]);
514+
V(i + 4 * num_groups, i + 3 * num_groups) =
515+
-criticalPerSevereAdjusted / (params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i]);
516+
517+
for (size_t j = 0; j < num_groups; j++) {
518+
V(i + 4 * num_groups, j + 4 * num_groups) = J(i, j);
519+
}
520+
}
521+
522+
V = V.inverse();
523+
524+
//Compute F*V
525+
Eigen::MatrixXd NextGenMatrix(num_infected_compartments * num_groups, 5 * num_groups);
526+
NextGenMatrix = F * V;
527+
528+
//Compute the largest eigenvalue in absolute value
529+
Eigen::ComplexEigenSolver<Eigen::MatrixXd> ces;
530+
531+
ces.compute(NextGenMatrix);
532+
const Eigen::VectorXcd eigen_vals = ces.eigenvalues();
533+
534+
Eigen::VectorXd eigen_vals_abs;
535+
eigen_vals_abs.resize(eigen_vals.size());
536+
537+
for (int i = 0; i < eigen_vals.size(); i++) {
538+
eigen_vals_abs[i] = std::abs(eigen_vals[i]);
539+
}
540+
return mio::success(eigen_vals_abs.maxCoeff());
541+
}
542+
/**
543+
*@brief Computes the reproduction number for all time points of the Model output obtained by the Simulation.
544+
*@param sim The Model Simulation.
545+
*@tparam Base simulation type that uses a SECIR compartment model. see Simulation.
546+
*@returns Eigen::Vector containing all reproduction numbers
547+
*/
548+
549+
template <class Base>
550+
Eigen::VectorXd get_reproduction_numbers(const Simulation<Base>& sim)
551+
{
552+
Eigen::VectorXd temp(sim.get_result().get_num_time_points());
553+
for (int i = 0; i < sim.get_result().get_num_time_points(); i++) {
554+
temp[i] = get_reproduction_number((size_t)i, sim).value();
555+
}
556+
return temp;
557+
}
558+
559+
/**
560+
*@brief @brief Computes the reproduction number at a given time point of the Simulation. If the particular time point is not part of the output, a linearly interpolated value is returned.
561+
*@param t_value The time point at which the reproduction number should be computed.
562+
*@param sim The Model Simulation.
563+
*@tparam Base simulation type that uses a SECIR compartment model. see Simulation.
564+
*@returns The computed reproduction number at the provided time point, potentially using linear interpolation.
565+
*/
566+
template <class Base>
567+
IOResult<ScalarType> get_reproduction_number(ScalarType t_value, const Simulation<Base>& sim)
568+
{
569+
if (t_value < sim.get_result().get_time(0) || t_value > sim.get_result().get_last_time()) {
570+
return mio::failure(mio::StatusCode::OutOfRange,
571+
"Cannot interpolate reproduction number outside computed horizon of the TimeSeries");
572+
}
573+
574+
if (t_value == sim.get_result().get_time(0)) {
575+
return mio::success(get_reproduction_number((size_t)0, sim).value());
576+
}
577+
578+
const auto& times = sim.get_result().get_times();
579+
580+
auto time_late = std::distance(times.begin(), std::lower_bound(times.begin(), times.end(), t_value));
581+
582+
ScalarType y1 = get_reproduction_number(static_cast<size_t>(time_late - 1), sim).value();
583+
ScalarType y2 = get_reproduction_number(static_cast<size_t>(time_late), sim).value();
584+
585+
auto result = linear_interpolation(t_value, sim.get_result().get_time(time_late - 1),
586+
sim.get_result().get_time(time_late), y1, y2);
587+
return mio::success(static_cast<ScalarType>(result));
588+
}
589+
331590
/**
332591
* Get migration factors.
333592
* Used by migration graph simulation.

0 commit comments

Comments
 (0)