Try adding mortality and reinfection rates to the basic SIR model

Mortality and reinfection in class cl_SIRbetaInput created in the previous article, "Implementing a mathematical model of infectious diseases" SIR model "in OpenModelica" Let's create a full version with the rate added.

I found an error in the mortality formula in the first post, so I am correcting it.
- (Corrected on May 13, 2020) Revised again.

Prerequisites for this model ――A certain percentage of people die from infection to recovery ――Even those who have acquired immunity once lose immunity at a certain rate and return to infectious people over time. --Only survivors return from recoverers to infectious individuals

** Please note that all of the following simulations are only verification of the model with temporary parameters and do not reflect or anticipate the actual values. ** **

Create class cl_SIRreinfect with mortality added to class cl_SIRbetaInput.

It is possible to set "dying_rate = mortality rate" as a parameter. The default value is 0.01 (1%).

parameter Real dying_rate = 0.01 "Dying rate";

equation
  der(Rdead) = Iy * gamma_const * dying_rate;
//  der(Rdead) = der(Ry) * dying_rate;
  Rlive = Ry - Rdead;

Add more reinfection rates

For the formula considering reinfection, refer to "R and Pandemic Mathematical Model New Coronavirus (COVID-19) Research as an Example". I did. If ρ is the reinfection rate, the formula is as follows.

\begin{align}

\frac{dS}{dt} &= -\beta SI + ρ R\\
\frac{dI}{dt} &=  \beta SI -\gamma I \\
\frac{dR}{dt} &=  \gamma I - ρ R\\

\end{align}

Allows you to set "reinfect_rate = reinfection rate" as a parameter. The default value is 0.001 (0.1%).

parameter Real reinfect_rate = 0.001 "Reinfection rate";

The additional equation parts are as follows When transitioning from I (infected person) to R (recoverer), the number of recoverers becomes death in proportion to the mortality rate, and the rest become "survival recovery people". Among "survival recoverers", S (potentially infectious) shifts in proportion to the reinfection rate, and R (recovery) decreases accordingly.

equation
  der(Sy) = (-contact_rate * Sy * Iy) + Rlive * reinfect_rate;
  der(Iy) = contact_rate * Sy * Iy - Iy * gamma_const;
  der(Ry) = Iy * gamma_const - Rlive * reinfect_rate;
  Inew = contact_rate * Sy * Iy;

The whole code looks like this.

class cl_SIRreinfect
  extends Modelica.Blocks.Icons.Block;
  parameter Real gamma_const = 0.2 "Recovery rate (1/day)";
  parameter Real reinfect_rate = 0.001 "Reinfection rate";
  parameter Real dying_rate = 0.01 "Dying rate";
  Modelica.Blocks.Interfaces.RealInput Si "Connector of initial S(Susceptible)" annotation(
    Placement(visible = true, transformation(origin = {-120, 60}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, 62}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealInput Ii "Connector of initial I (Infectious)" annotation(
    Placement(visible = true, transformation(origin = {-120, 0}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, 0}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealInput Ri "Connector of initial R (Removed)" annotation(
    Placement(visible = true, transformation(origin = {-120, -60}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, -62}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealOutput Sy "Output connector of S(Susceptible)" annotation(
    Placement(visible = true, transformation(origin = {110, 60}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, 62}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealOutput Iy "Output connector of I (Infectious)" annotation(
    Placement(visible = true, transformation(origin = {110, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealOutput Ry "Output connector of R (Removed)" annotation(
    Placement(visible = true, transformation(origin = {110, -60}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, -60}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealInput contact_rate "Connector of infectious rate" annotation(
    Placement(visible = true, transformation(origin = {0, 120}, extent = {{-20, -20}, {20, 20}}, rotation = -90), iconTransformation(origin = {0, 120}, extent = {{-20, -20}, {20, 20}}, rotation = -90)));
  Modelica.Blocks.Interfaces.RealOutput Inew "Output connector of New I(Infectious)" annotation(
    Placement(visible = true, transformation(origin = {-60, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90), iconTransformation(origin = {-60, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));
  Modelica.Blocks.Interfaces.RealOutput Rlive "Output connector of R (Live)" annotation(
    Placement(visible = true, transformation(origin = {60, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90), iconTransformation(origin = {60, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));
  Modelica.Blocks.Interfaces.RealOutput Rdead "Output connector of R (Dead)" annotation(
    Placement(visible = true, transformation(origin = {0, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90), iconTransformation(origin = {0, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));
initial equation
  Sy = Si;
  Iy = Ii;
  Ry = Ri;
  Rdead = 0;
equation
  der(Sy) = (-contact_rate * Sy * Iy) + Rlive * reinfect_rate;
  der(Iy) = contact_rate * Sy * Iy - Iy * gamma_const;
  der(Ry) = Iy * gamma_const - Rlive * reinfect_rate;
  Inew = contact_rate * Sy * Iy;
  der(Rdead) = Iy * gamma_const * dying_rate;
//  der(Rdead) = der(Ry) * dying_rate;
  Rlive = Ry - Rdead;
  annotation(
    Diagram,
    Icon(coordinateSystem(initialScale = 0.1), graphics = {Text(origin = {-62, 53}, extent = {{-52, 29}, {174, -133}}, textString = "SIR\nReI")}));
end cl_SIRreinfect;

Simulation results including mortality rate and reinfection rate

Simulation with only mortality set

Use the created cl_SIRreinfect to create and simulate the model mdl_SIRmodel_reinfect. Vaccines are not used here.

For clarity, we will simulate with a fairly high mortality rate of 20%.

As a result, the number of deaths has converged to about 18% of the total.

Simulation by setting only the reinfection rate

Next, let's simulate with a model in which only the reinfection rate is set. Try a higher 1% reinfection rate.

In the simulation up to 100 days, the number of people who can be infected continues to increase and it is not possible to see what will happen in the future. When simulating up to 500 days, S, I, and R all oscillate and gradually converge. Looking at I (number of infected people), we can see that there are 2nd and 3rd waves.

Try setting both mortality and reinfection rates

Try setting the above two parameters at the same time.

- (Revised on May 13, 2020 and restored) In the simulation, the number of deaths gradually increased while re-infecting, and about 59% died, and the number of people who could be infected increased while I (infected person) and RLive (survival recovery person) converged to zero. It has converged to just under 40%. It should be a good formula, but it feels different from an intuitive graph. If you have any problems with how to formulate the formula, please point it out.

Summary

By modeling with OpenModelica, I think it was easy to add elements to the basic SIR model little by little and simulate it. Actually, it seems that both the mortality rate and the reinfection rate are more than an order of magnitude lower than this simulation, but if more accurate actual data can be obtained, it may be possible to identify the parameters.

Implementing the mathematical model "SIR model" of infectious diseases with OpenModelica Implementing a mathematical model "SIR model" of infectious diseases with OpenModelica (see the effect of the vaccine) Introduction of mathematical prediction model for infectious diseases (SIR model)

[PYTHON] Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (reflecting mortality rate and reinfection rate)