The Modelica epidemiological model is wrong

[mbaudin47]

The epidemiological model implemented in epid.mo is wrong. Indeed the simplest SIR equations are (see ¹ page 455):

$$ \begin{aligned} &\dfrac{dS}{dt} = -\beta \dfrac{SI}{N}\\ &\dfrac{dI}{dt} = \beta \dfrac{SI}{N} - \gamma I\\ &\dfrac{dR}{dt} = \gamma I \end{aligned} $$

But the equation for $\frac{dS}{dt}$ is implemented as:

otfmi/otfmi/example/file/epid.mo

Line 18 in 6716944

der(susceptible) = - infection_rate*infected*susceptible;

The significant fact is that the right hand side of the ODE in the model does involve N: the model uses infection_rate * infected * susceptible (i.e. $\beta SI$) instead of infection_rate * infected * susceptible / total_pop (i.e. $\beta SI / N$). This is wrong because it is not consistent with the initial number of infected $I(0)$ which is set to 1:

otfmi/otfmi/example/file/epid.mo

Line 13 in 6716944

infected = 1;

Actually, there are two different consistent choices.

• We set the total population $N$ to 1. In this case, S, I and R represent the fraction of the population in each compartment. Under this hypothesis, the initial number of infected is set to $I(0) = \frac{1}{N}$ and $S(0) = 1 - I(0)$. In this case, the line 13 of the model is wrong and the expression for der(susceptible), der(infected) and der(removed) in the .mo file is correct.
• We use the actual total population. In this case, we must divide by $N$. Under this hypothesis, the expression for der(susceptible), der(infected) and der(removed) in the current .mo file is wrong. This is the choice consistent with ² figure 2 page 18 (the error must have been introduced later, because the Modelica model presented in ² page 16 is correct).

The consequence of the bug is the the beta variable does not have the value it should: the current code uses $\beta / N$ where it should use $\beta$. Hence, the basic reproduction number $R_0 = \beta / \gamma$ is wrong. It is currently close to 0.005456, which is much smaller than the expected value of influenza. The correct value should be close to 4.163, i.e. 763 times much larger.

All in all, I suggest to use the following model, where the parameters are the infectious period (which is 1 / gamma) and the reproduction number (R0 = beta / gamma). These parameters are much easier to interpret and standard in the bibliography. Below, I use the parameters found in ¹ figure 1 page 455. This leads to beta = 0.2143 (infection rate) and gamma = 1.0 / 14 = 0.07143 (healing rate).

model epid

parameter Real total_pop = 1000.0;
parameter Real healing_rate = 0.07143;
parameter Real infection_rate = 0.2143;

Real infected;
Real susceptible;
Real removed;

initial equation
infected = 1;
removed = 0;
total_pop = infected + susceptible + removed;

equation
der(susceptible) = - infection_rate * infected * susceptible / total_pop;
der(infected) = infection_rate * infected * susceptible / total_pop - healing_rate * infected;
der(removed) = healing_rate * infected;

annotation(
    experiment(StartTime = 0.0, StopTime = 200.0, Tolerance = 1e-6, Interval = 0.1));
end epid;

When we simulate the model, we can compare to the results found at https://shiny.bcgsc.ca/posepi1/.

If we want to reproduce the results from ², we may use the following settings, obtained after calibrating from the data found in ³ using non linear least squares, and rounded to 1 significant digits (I believe that the extra digits are not accurate, given the uncertainties in the data):

parameter Real total_pop = 763.0;
parameter Real healing_rate = 0.5;
parameter Real infection_rate = 2.0;

Updating these model will require to update the .fmu for Linux and Windows platforms available in the FMU directory.

Footnotes

1. Bjørnstad, O. N., Shea, K., Krzywinski, M., & Altman, N. (2020). Modeling infectious epidemics. Nature methods, 17(5), 455-457. ↩ ↩²

2. S. Girard, “A probabilistic take on system modeling with Modelica and Python”, Février 2019, https://sylvaingirard.net/pdf/girard17-probabilistic_modelica_python.pdf ↩ ↩² ↩³

3. Anonymous (1978). ‘Influenza in a boarding school’. In: British Medical Journal 587.1. url: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1603269/?page=2. ↩

[SG-phimeca]

The model as programmed in ¹ is taken from ², page 713 (page 15 of pdf refered below).
The system of equation is indeed the same as the one you wrote in the beginning of the issue.

The Modelica model used (for direct simulation -- the model for parameter estimation as a slightly different initialization) in ¹ was

model SIR "epidemic model"
  parameter Real total_population=763;
  Real susceptible;
  Real infected;
  Real removed;
  input Real contact_rate (start=0.5);
  input Real average_infectious_period (start=2.5, fixed=true);
  parameter Real infected_initial=1;
  initial equation
  total_population = susceptible + infected + removed;
  infected = infected_initial;

  equation

  der(susceptible) = -contact_rate * susceptible * infected /
    total_population;
  der(infected) = contact_rate * susceptible * infected / total_population -
    infected / average_infectious_period;
  der(removed) = infected / average_infectious_period;

end SIR;

I think the total population (N) should be a number of people (not 1).

I believe the bug is actually an “inclusion” of the 1/N factor into the Modelica parameter infection_rate. It may have been intentional;
I do not remember the rationale.
I agree that sticking to the classic parameterization would be more straightforward.

Footnotes

1. S. Girard, “A probabilistic take on system modeling with Modelica and Python”, Février 2019, https://sylvaingirard.net/pdf/girard17-probabilistic_modelica_python.pdf ↩ ↩²

2. Kermack, William O. and McKendrick, Anderson G., “A contribution to the mathematical theory of epidemics”, 1927
   Kermack27_-_A_contribution_to_the_mathematical_theory_of_epidemics.pdf ↩

[mbaudin47]

The model you suggest seems to be excellent from my point of view. Moreover, the fact that we can set infected_initial in the new model you suggest is very good (it is not possible in the current one, as far as I understand). In practical settings, this might be interesting to calibrate this parameter as well.

[mbaudin47]

The problem is that the current help page is not correct.

1. The plot was not updated according to the change in the model:

https://openturns.github.io/otfmi/master/_images/epid.png

Figure 1. The figure in the current doc.

The problems are:

• the X axis is not seconds, but 15 days,
• the S, I and R values in the Y axis are not consistent with the model: the value of $I$ should go up to approximately 300 at the maximum.

The plot should be more like:

Figure 2. The default model in https://shiny.bcgsc.ca/posepi1/.

2. The reference to the paper providing the calibrating data is not introduced in the help page:

Anonymous (1978). ‘Influenza in a boarding school’. In: British Medical Journal 587.1. url: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1603269/?page=2

This reference must be presented in the help page, because the value of the parameters in the model are calibrated so that it fits this particular data. This is not a generic SIR model, it is a model calibrated to the influenza propagation in this particular school in 1978.

Here is the figure from the paper.

Figure 3. The data from the (Anonymous, 1978) paper.

[mbaudin47]

@jschueller : Thank you for updating the plot! This fixes the first part of the issue:

• The plot was not updated according to the change in the model.
• The reference to the paper providing the calibrating data is not introduced in the help page.

But the second topic of the issue should also be covered in my opinion.