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:

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:

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)$ and the line 13 of the model is wrong. Hence 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 these settings that I 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.