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

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.

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";

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".

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. 20_SIR_full_model.png

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

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

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. 23_SIR_full_model_prm.png

24_SIR_full_reinfect_graph.png

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. 25_SIR_full_reinfect_graph.png

Try setting both mortality and reinfection rates

Try setting the above two parameters at the same time. 28_SIR_full_model_prm.png

29_SIR_full_reinfect_graph.png

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.

Related article

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)

Recommended Posts

Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (reflecting mortality rate and reinfection rate)
Implement the mathematical model "SIR model" of infectious diseases with OpenModelica
Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (see the effect of the vaccine)
Introduction of mathematical prediction model for infectious diseases (SIR model)
Mathematical prediction model for infectious diseases (SIR model): Case study (1)
[SIR model analysis] Peak out of infection numbers in Japan and around the world ♬
Implement the mathematical model "SIR model" of infectious diseases with OpenModelica
Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (reflecting mortality rate and reinfection rate)
Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (see the effect of the vaccine)
Introduction of mathematical prediction model for infectious diseases (SIR model)
Implement a model with state and behavior (3) --Example of implementation by decorator
Mathematical prediction model for infectious diseases (SIR model): Case study (1)
Mathematical model of infectious disease epidemics
Implement a model with state and behavior
Predict infectious disease epidemics with SIR model
Solving the Lorenz 96 model with Julia and Python
Example of reading and writing CSV with Python
[SIR model analysis] Peak out of infection numbers in Japan and around the world ♬
Scraping the schedule of Hinatazaka46 and reflecting it in Google Calendar
Verify the compression rate and time of PIXZ used in practice
Mathematical model of infectious disease epidemics
Implement part of the process in C ++
Sum of variables in a mathematical model
[SIR model analysis] Peak out of infection numbers in Japan and around the world ♬
[SIR model analysis] Determination of γ * (R-1) and peak out of infection number ♬ World edition
[SIR model analysis] Peak out of infections in various parts of Japan ♬
Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (reflecting mortality rate and reinfection rate)
Open an Excel file in Python and color the map of Japan
Display the status of COVID 19 infection in Japan with Splunk (GitHub version)
Implement part of the process in C ++
Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (see the effect of the vaccine)
[SIR model analysis] Transform the formula to determine γ and the effective reproduction number R ♬
I used Python to find out about the role choices of the 51 "Yachts" in the world.